rplus.hpp

Go to the documentation of this file.

/*
 * Copyright © 2025-2026 Akreon LLC
 * SPDX-License-Identifier: Apache-2.0
 */
 
#pragma once
 
#include "cone.hpp"
#include "akreon/linalg/operator/quadratic_rep.hpp"
 
namespace conic {
 
template <linalg::VectorContainer Vec>
struct RplusCirc {
  using T = linalg::vector_element_t<Vec>;
  static constexpr std::size_t N = linalg::vector_dim_v<Vec>;
 
  constexpr auto operator()(const Vec& x, const Vec& y) const {
    Vec res{};
    for (std::size_t i = 0; i < N; ++i) res(i) = x(i) * y(i);
    return res;
  }
};
 
template <linalg::VectorContainer Slack>
struct NonnegOrthantCone
    : public Cone<Slack, NonnegOrthantCone<Slack>, ConeType::SYMMETRIC> {
  using ConeBase = Cone<Slack, NonnegOrthantCone<Slack>, ConeType::SYMMETRIC>;
 
  using ConeBase::eps;
  using ConeBase::M;
  using typename ConeBase::HessianMat;
  using typename ConeBase::T;
 
  constexpr Slack x0_impl() const {
    Slack s{};
    for (std::size_t i = 0; i < M; ++i) s(i) = T(1);
    return s;
  }
 
  constexpr Slack y0_impl() const { return x0_impl(); }
 
  constexpr bool in_cone_impl(const Slack& x) const {
    for (std::size_t i = 0; i < M; ++i)
      if (!(x(i) > T(0))) return false;
    return true;
  }
 
  constexpr bool in_dual_cone_impl(const Slack& y) const {
    return in_cone_impl(y);
  }
 
  constexpr T alpha_max_impl(const Slack& x,
                             const linalg::Vector<T, M>& p) const {
    constexpr T eta{0.99};
    T a{1};
    for (std::size_t i = 0; i < M; ++i) {
      if (p(i) < T(0)) {
        const T cand = -x(i) / p(i);
        a = (cand < a) ? cand : a;
      }
    }
    return std::min(T(1), eta * a);
  }
 
  constexpr T alpha_max_dual_impl(const Slack& y,
                                  const linalg::Vector<T, M>& p) const {
    return alpha_max_impl(y, p);
  }
 
  constexpr auto H_impl(const Slack& x, const Slack& y) const {
    const auto [w, event] = w_impl(x, y);
    if (linalg::is_fatal_event(event))
      return std::tuple(
          linalg::QuadraticRep<HessianMat, RplusCirc<Slack>>(e_impl()), event);
    Slack w_inv{};
    bool safeguarded = false;
    for (std::size_t i = 0; i < M; ++i) {
      const T wi_safe = std::max(w(i), eps);
      safeguarded = safeguarded || (wi_safe != w(i));
      w_inv(i) = T(1) / wi_safe;
    }
    const auto out_event =
        linalg::is_fatal_event(event) ? event
        : (event != linalg::NumericalEvent::NONE || safeguarded)
            ? linalg::NumericalEvent::NUMERICAL_SAFEGUARD
            : linalg::NumericalEvent::NONE;
    return std::tuple(linalg::QuadraticRep<HessianMat, RplusCirc<Slack>>(w_inv),
                      out_event);
  }
 
  constexpr auto grad_impl(const Slack& x) const {
    Slack dx{};
    bool safeguarded = false;
    for (std::size_t i = 0; i < M; ++i) {
      const T xi_safe = std::max(x(i), eps);
      safeguarded = safeguarded || (xi_safe != x(i));
      dx(i) = T(-1) / xi_safe;
    }
    return std::tuple(dx, safeguarded
                              ? linalg::NumericalEvent::NUMERICAL_SAFEGUARD
                              : linalg::NumericalEvent::NONE);
  }
 
  constexpr T nu_impl() const { return static_cast<T>(M); }
 
  constexpr auto e_impl() const {
    linalg::vector_traits<Slack> vec{};
    return vec.fill(T(1));
  }
 
  constexpr auto w_impl(const Slack& x, const Slack& y) const {
    Slack w{};
    bool safeguarded{false};
    for (std::size_t i = 0; i < M; ++i) {
      const T x_safe = std::max(x(i), eps);
      const T y_safe = std::max(y(i), eps);
      safeguarded = safeguarded || (x_safe != x(i)) || (y_safe != y(i));
      w(i) = linalg::constexpr_sqrt_safe(x_safe / y_safe);
    }
    return std::tuple(w, safeguarded
                             ? linalg::NumericalEvent::NUMERICAL_SAFEGUARD
                             : linalg::NumericalEvent::NONE);
  }
 
  constexpr auto circ_impl(const Slack& x, const Slack& y) const {
    return RplusCirc<Slack>()(x, y);
  }
 
  template <typename Scaling>
  constexpr auto scaled_affine_correction_impl(const Slack& x, const Scaling&,
                                               const Slack& px,
                                               const Slack& py) const {
    Slack correction{};
    bool safeguarded = false;
    for (std::size_t i = 0; i < M; ++i) {
      const T xi_safe = std::max(x(i), eps);
      safeguarded = safeguarded || (xi_safe != x(i));
      correction(i) = px(i) * py(i) * (T(1) / xi_safe);
    }
    return std::tuple(correction,
                      safeguarded ? linalg::NumericalEvent::NUMERICAL_SAFEGUARD
                                  : linalg::NumericalEvent::NONE);
  }
 
  constexpr auto trace_impl(const Slack& x, const Slack& y) const {
    const auto x_circ_y = circ_impl(x, y);
    T res{};
    for (std::size_t i = 0; i < M; ++i) res += x_circ_y(i);
    return res;
  }
 
  constexpr auto inv_impl(const Slack& x) const {
    Slack res{};
    bool safeguarded = false;
    for (std::size_t i = 0; i < M; ++i) {
      const T xi_safe = std::max(x(i), eps);
      safeguarded = safeguarded || (xi_safe != x(i));
      res(i) = T(1) / xi_safe;
    }
    return std::tuple(res, safeguarded
                               ? linalg::NumericalEvent::NUMERICAL_SAFEGUARD
                               : linalg::NumericalEvent::NONE);
  }
 
  constexpr auto gamma_F_impl(const Slack& x, const Slack& y) const {
    T res{};
    bool safeguarded = false;
    for (std::size_t i = 0; i < M; ++i) {
      const auto xy_safe = std::max(x(i) * y(i), eps);
      safeguarded = safeguarded || (xy_safe != x(i) * y(i));
      res += linalg::constexpr_ln_safe(xy_safe, eps);
    }
    const auto xTy_safe =
        std::max(linalg::inner_product(x, y), eps) / nu_impl();
    safeguarded =
        safeguarded || (xTy_safe * nu_impl() != linalg::inner_product(x, y));
    res = nu_impl() * linalg::constexpr_ln_safe(xTy_safe, eps) - res;
    return std::tuple(res, safeguarded
                               ? linalg::NumericalEvent::NUMERICAL_SAFEGUARD
                               : linalg::NumericalEvent::NONE);
  }
};
 
} // namespace conic