ell_core.hpp

Go to the documentation of this file.

 
// -*- coding: utf-8 -*-
#pragma once
 
#include <limits>
#include <utility>
#include <valarray>
 
#include "ell_calc.hpp"
#include "ell_config.hpp"
#include "ell_matrix.hpp"
 
class EllCore {
    using Vec = std::valarray<double>;
 
    size_t _n;
    double _kappa;
    Matrix _mq;
    EllCalc _helper;
    double _tsq{};
    // Pre-allocated scratch buffers (Rust-style: eliminate per-call allocation)
    Vec _scratch;  // w = L^{-1}g (stable) / grad_t (non-stable)
    Vec _z;        // z = D^{-1}w (stable)
    Vec _v;        // q = L^{-T}z then v (stable rank-1 update)
 
  public:
    bool no_defer_trick = false;
 
  private:
    auto operator=(const EllCore& E) -> EllCore& = delete;
 
    EllCore(double kappa, Matrix&& mq, size_t ndim)
        : _n{ndim},
          _kappa{kappa},
          _mq{std::move(mq)},
          _helper{_n},
          _scratch(0.0, ndim),
          _z(0.0, ndim),
          _v(0.0, ndim) {}
 
  public:
    EllCore(const Vec& val, size_t ndim) : EllCore{1.0, Matrix(ndim), ndim} {
        this->_mq.diagonal() = val;
    }
 
    EllCore(const double alpha, const size_t ndim) : EllCore{alpha, Matrix(ndim), ndim} {
        this->_mq.identity();
    }
 
    EllCore(EllCore&& E) = default;
 
    EllCore& operator=(EllCore&&) = default;
 
    ~EllCore() = default;
 
    explicit EllCore(const EllCore& E) = default;
 
    auto copy() const -> EllCore { return EllCore(*this); }
 
    constexpr auto tsq() const -> double { return this->_tsq; }
 
    void set_use_parallel_cut(bool value) { this->_helper.use_parallel_cut = value; }
 
    template <typename T> auto update_bias_cut(Vec& grad, const T& beta) -> CutStatus {
        return this->_update_core(grad, beta, [this](const T& beta_l, const double tsq_l) {
            return this->_update_cut_bias_cut(beta_l, tsq_l);
        });
    }
 
    template <typename T> auto update_central_cut(Vec& grad, const T& beta) -> CutStatus {
        return this->_update_core(grad, beta, [this](const T& beta_l, const double tsq_l) {
            return this->_update_cut_central_cut(beta_l, tsq_l);
        });
    }
 
    template <typename T> auto update_q(Vec& grad, const T& beta) -> CutStatus {
        return this->_update_core(grad, beta, [this](const T& beta_l, const double tsq_l) {
            return this->_update_cut_q(beta_l, tsq_l);
        });
    }
 
    template <typename T> auto update_stable_bias_cut(Vec& grad, const T& beta) -> CutStatus {
        return this->_update_stable_core(grad, beta, [this](const T& beta_l, const double tsq_l) {
            return this->_update_cut_bias_cut(beta_l, tsq_l);
        });
    }
 
    template <typename T> auto update_stable_central_cut(Vec& grad, const T& beta) -> CutStatus {
        return this->_update_stable_core(grad, beta, [this](const T& beta_l, const double tsq_l) {
            return this->_update_cut_central_cut(beta_l, tsq_l);
        });
    }
 
    template <typename T> auto update_stable_q(Vec& grad, const T& beta) -> CutStatus {
        return this->_update_stable_core(grad, beta, [this](const T& beta_l, const double tsq_l) {
            return this->_update_cut_q(beta_l, tsq_l);
        });
    }
 
  private:
    template <typename T, typename Fn>
    auto _update_core(Vec& grad, const T& beta, Fn&& cut_strategy) -> CutStatus {
        // reuse _v as grad_t = M * g (pre-allocated, no allocation)
        for (size_t i = 0; i != this->_n; ++i) {
            auto s = 0.0;
            for (size_t j = 0; j != this->_n; ++j) {
                s += this->_mq(i, j) * grad[j];
            }
            this->_v[i] = s;
        }
 
        auto omega = 0.0;
        for (size_t i = 0; i != this->_n; ++i) {
            omega += this->_v[i] * grad[i];
        }
        this->_tsq = this->_kappa * omega;
 
        if (omega <= std::numeric_limits<double>::min()) {
            grad = this->_v;
            return CutStatus::NoEffect;
        }
 
        auto result = std::forward<Fn>(cut_strategy)(beta, this->_tsq);
        if (result.status != CutStatus::Success) {
            return result.status;
        }
 
        // n (n+1) / 2 + n
        const auto r = result.sigma / omega;
        for (size_t i = 0; i != this->_n; ++i) {
            const auto rQg = r * this->_v[i];
            for (size_t j = 0; j != i; ++j) {
                this->_mq(i, j) -= rQg * this->_v[j];
                this->_mq(j, i) = this->_mq(i, j);
            }
            this->_mq(i, i) -= rQg * this->_v[i];
        }
 
        this->_kappa *= result.delta;
 
        if (this->no_defer_trick) {
            this->_mq *= this->_kappa;
            this->_kappa = 1.0;
        }
 
        const auto rho_over_omega = result.rho / omega;
        for (size_t i = 0; i != this->_n; ++i) {
            grad[i] = this->_v[i] * rho_over_omega;
        }
        return result.status;
    }
 
    template <typename T, typename Fn>
    auto _update_stable_core(Vec& g, const T& beta, Fn&& cut_strategy) -> CutStatus {
        // _scratch = w = L^{-1}g (forward substitution)
        this->_scratch = g;
        for (size_t j = 0; j != this->_n - 1; ++j) {
            for (size_t i = j + 1; i != this->_n; ++i) {
                this->_mq(j, i) = this->_mq(i, j) * this->_scratch[j];
                this->_scratch[i] -= this->_mq(j, i);
            }
        }
 
        // _z = D^{-1} * w
        this->_z = this->_scratch;
        for (size_t i = 0; i != this->_n; ++i) {
            this->_z[i] *= this->_mq(i, i);
        }
 
        // omega = sum(w_i * z_i)
        auto omega = 0.0;
        for (size_t i = 0; i != this->_n; ++i) {
            omega += this->_z[i] * this->_scratch[i];
        }
 
        this->_tsq = this->_kappa * omega;
 
        if (omega <= std::numeric_limits<double>::min()) {
            return CutStatus::NoEffect;
        }
 
        auto result = std::forward<Fn>(cut_strategy)(beta, this->_tsq);
        if (result.status != CutStatus::Success) {
            return result.status;
        }
 
        // _v = grad_t = L^{-T} * z (back substitution) — preserved for final g
        this->_v = this->_z;
        for (auto i = this->_n - 1; i != 0; --i) {
            for (auto j = i; j != this->_n; ++j) {
                this->_v[i - 1] -= this->_mq(j, i - 1) * this->_v[j];
            }
        }
 
        // rank-one LDL^T update — reuse _scratch as working vector v
        // (_scratch (w) no longer needed; overwrite with gradient)
        const auto mu = result.sigma / (1.0 - result.sigma);
        auto oldt = omega / mu;
        this->_scratch = g;  // v = gradient
        for (size_t j = 0; j != this->_n; ++j) {
            const auto p = this->_scratch[j];
            const auto temp = this->_z[j];
            const auto newt = oldt + p * temp;
            const auto beta2 = temp / newt;
            this->_mq(j, j) *= oldt / newt;
            for (auto k = j + 1; k != this->_n; ++k) {
                this->_scratch[k] -= this->_mq(j, k);
                this->_mq(k, j) += beta2 * this->_scratch[k];
            }
            oldt = newt;
        }
        this->_kappa *= result.delta;
        // _v still holds grad_t = L^{-T}*z (preserved from back substitution)
        const auto rho_over_omega = result.rho / omega;
        for (size_t i = 0; i != this->_n; ++i) {
            g[i] = this->_v[i] * rho_over_omega;
        }
        return result.status;
    }
 
    auto _update_cut_bias_cut(double beta, double tsq) const -> CutResult {
        return this->_helper.calc_bias_cut(beta, tsq);
    }
 
    auto _update_cut_bias_cut(const std::valarray<double>& beta, double tsq) const -> CutResult {
        if (beta.size() < 2) {
            return this->_helper.calc_bias_cut(beta[0], tsq);
        }
        return this->_helper.calc_parallel_bias_cut(beta[0], beta[1], tsq);
    }
 
    auto _update_cut_central_cut(double /*unused*/, double tsq) const -> CutResult {
        return this->_helper.calc_central_cut(tsq);
    }
 
    auto _update_cut_central_cut(const std::valarray<double>& beta, double tsq) const -> CutResult {
        if (beta.size() < 2) {
            return this->_helper.calc_central_cut(tsq);
        }
        return this->_helper.calc_parallel_central_cut(beta[1], tsq);
    }
 
    auto _update_cut_q(double beta, double tsq) const -> CutResult {
        return this->_helper.calc_bias_cut_q(beta, tsq);
    }
 
    auto _update_cut_q(const std::valarray<double>& beta, double tsq) const -> CutResult {
        if (beta.size() < 2) {
            return this->_helper.calc_bias_cut_q(beta[0], tsq);
        }
        return this->_helper.calc_parallel_bias_cut_q(beta[0], beta[1], tsq);
    }
 
};  // } EllCore