cone.hpp

Go to the documentation of this file.

/*
 * Copyright © 2025-2026 Akreon LLC
 * SPDX-License-Identifier: Apache-2.0
 */
 
#pragma once
 
#include <climits>
#include <optional>
 
#include "akreon/convex/cone_traits.hpp"
#include "akreon/linalg/operator/hessian.hpp"
#include "akreon/linalg/operator/ldl.hpp"
 
namespace conic {
 
template <linalg::VectorContainer Slack, typename Derived, ConeType CType>
struct Cone {
  using T = linalg::vector_element_t<Slack>;
 
  static constexpr std::size_t M = linalg::vector_dim_v<Slack>;
 
  static constexpr T eps = 1e-8;
 
  using HessianMat = linalg::Matrix<T, M, M>;
 
  constexpr Slack x0() const { return self().x0_impl(); }
 
  constexpr Slack y0() const { return self().y0_impl(); }
 
  constexpr bool in_cone(const Slack& x) const {
    return self().in_cone_impl(x);
  }
 
  constexpr bool in_dual_cone(const Slack& y) const {
    return self().in_dual_cone_impl(y);
  }
 
  constexpr auto H(const Slack& x, const Slack& y) const {
    return self().H_impl(x, y);
  }
 
  constexpr T alpha_max(const Slack& x, const linalg::Vector<T, M>& p) const {
    return self().alpha_max_impl(x, p);
  }
 
  constexpr T alpha_max_dual(const Slack& y,
                             const linalg::Vector<T, M>& p) const {
    return self().alpha_max_dual_impl(y, p);
  }
 
  constexpr T nu() const { return self().nu_impl(); }
 
  constexpr auto grad(const Slack& x) const { return self().grad_impl(x); }
 
  constexpr auto e() const {
    static_assert(ConeType::SYMMETRIC == CType,
                  "Nonsymmetric cone is not equipped with Jordan Algebras");
    return self().e_impl();
  }
 
  constexpr auto w(const Slack& x, const Slack& y) const {
    static_assert(ConeType::SYMMETRIC == CType,
                  "Nonsymmetric cone is not equipped with Jordan Algebras");
    return self().w_impl(x, y);
  }
 
  constexpr auto circ(const Slack& x, const Slack& y) const {
    static_assert(ConeType::SYMMETRIC == CType,
                  "Nonsymmetric cone is not equipped with Jordan Algebras");
    return self().circ_impl(x, y);
  }
 
  constexpr auto trace(const Slack& x, const Slack& y) const {
    static_assert(ConeType::SYMMETRIC == CType,
                  "Nonsymmetric cone is not equipped with Jordan Algebras");
    return self().trace_impl(x, y);
  }
 
  constexpr auto inv(const Slack& x) const {
    static_assert(ConeType::SYMMETRIC == CType,
                  "Nonsymmetric cone is not equipped with Jordan Algebras");
    return self().inv_impl(x);
  }
 
  constexpr T trace_weight() const {
    const auto identity = e();
    return trace(identity, identity) /
           linalg::inner_product(identity, identity);
  }
 
  template <typename Scaling>
  constexpr auto scaled_affine_correction(const Slack& x, const Scaling& P,
                                          const Slack& px,
                                          const Slack& py) const {
    static_assert(ConeType::SYMMETRIC == CType,
                  "NT correction requires a symmetric cone");
    return self().scaled_affine_correction_impl(x, P, px, py);
  }
 
  template <typename Scaling>
  constexpr auto scaled_affine_correction_impl(const Slack& x, const Scaling& P,
                                               const Slack& px,
                                               const Slack& py) const {
    const auto [root, root_event] = self().sqrt_impl(P.x);
    if (linalg::is_fatal_event(root_event))
      return std::tuple(Slack{}, root_event);
    const auto [inv_root, inv_event] = inv(root);
    if (linalg::is_fatal_event(inv_event))
      return std::tuple(Slack{}, inv_event);
    const auto W_inv = Scaling(root);
    const auto W = Scaling(inv_root);
    const auto lambda = W_inv.apply(x);
    HessianMat jordan{};
    for (std::size_t j = 0; j < M; ++j) {
      Slack basis{};
      basis(j) = T(1);
      const auto column = circ(lambda, basis);
      for (std::size_t i = 0; i < M; ++i) jordan(i, j) = column(i);
    }
    const auto [correction, event] =
        linalg::LDL(jordan).solve(circ(W_inv.apply(px), W.apply(py)));
    if (linalg::is_fatal_event(event)) return std::tuple(Slack{}, event);
    return std::tuple(W_inv.apply(correction), event);
  }
 
  constexpr auto gamma_F(const Slack& x, const Slack& y) const {
    static_assert(ConeType::SYMMETRIC == CType,
                  "Nonsymmetric cone is not equipped with Jordan Algebras");
    return self().gamma_F_impl(x, y);
  }
 
  constexpr auto self() const { return static_cast<const Derived&>(*this); }
};
 
} // namespace conic