Program Listing for File ell.hpp

↰ Return to documentation for file (ellcpp/ell.hpp)

// -*- coding: utf-8 -*-
#pragma once

#include <cmath>
#include <ellcpp/utility.hpp>
#include <tuple>
#include <xtensor/xarray.hpp>

// forward declaration
enum class CUTStatus;

class ell
{
  public:
    using Arr = xt::xarray<double, xt::layout_type::row_major>;
    // using params_t = std::tuple<double, double, double>;
    // using return_t = std::tuple<int, params_t>;

    bool use_parallel_cut = true;
    bool no_defer_trick = false;

  protected:
    double _mu {};
    double _rho {};
    double _sigma {};
    double _delta {};
    double _tsq {};

    const int _n;

    const double _nFloat;
    const double _nPlus1;
    const double _nMinus1;
    const double _halfN;
    const double _halfNplus1;
    const double _halfNminus1;
    const double _nSq;
    const double _c1;
    const double _c2;
    const double _c3;

    double _kappa;
    Arr _Q;
    Arr _xc;

    auto operator=(const ell& E) -> ell& = delete;

    template <typename V, typename U>
    ell(V&& kappa, Arr&& Q, U&& x) noexcept
        : _n {int(x.size())}
        , _nFloat {double(_n)}
        , _nPlus1 {_nFloat + 1.}
        , _nMinus1 {_nFloat - 1.}
        , _halfN {_nFloat / 2.}
        , _halfNplus1 {_nPlus1 / 2.}
        , _halfNminus1 {_nMinus1 / 2.}
        , _nSq {_nFloat * _nFloat}
        , _c1 {_nSq / (_nSq - 1)}
        , _c2 {2. / _nPlus1}
        , _c3 {_nFloat / _nPlus1}
        , _kappa {std::forward<V>(kappa)}
        , _Q {std::move(Q)}
        , _xc {std::forward<U>(x)}
    {
    }

  public:
    ell(const Arr& val, Arr x) noexcept
        : ell {1., xt::diag(val), std::move(x)}
    {
    }

    ell(const double& alpha, Arr x) noexcept
        : ell {alpha, xt::eye(x.size()), std::move(x)}
    {
    }

    ell(ell&& E) = default;

    ~ell() { }

    explicit ell(const ell& E) = default;

    [[nodiscard]] auto copy() const -> ell
    {
        return ell(*this);
    }

    [[nodiscard]] auto xc() const -> Arr
    {
        return _xc;
    }

    void set_xc(const Arr& xc)
    {
        _xc = xc;
    }

    template <typename T>
    auto update(const std::tuple<Arr, T>& cut) -> std::tuple<CUTStatus, double>;

  protected:
    auto _update_cut(const double& beta) -> CUTStatus
    {
        return this->_calc_dc(beta);
    }

    CUTStatus _update_cut(const Arr& beta)
    { // parallel cut
        if (beta.shape()[0] < 2)
        {
            return this->_calc_dc(beta[0]);
        }
        return this->_calc_ll_core(beta[0], beta[1]);
    }

    auto _calc_ll_core(const double& b0, const double& b1) -> CUTStatus;

    auto _calc_ll_cc(const double& b1, const double& b1sq) -> void;

    auto _calc_dc(const double& beta) noexcept -> CUTStatus;

    auto _calc_cc(const double& tau) noexcept -> void;
}; // } ell