randomization_block.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
//randomization_block.hpp is free software: you can redistribute it and/or modify
//it under the terms of the GNU General Public License as published by
//the Free Software Foundation, either version 3 of the License, or
//(at your option) any later version.
//
//randomization_block.hpp is distributed in the hope that it will be useful,
//but WITHOUT ANY WARRANTY; without even the implied warranty of
//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//GNU General Public License for more details.
//
//You should have received a copy of the GNU General Public License
//along with randomization_block.hpp.  If not, see <http://www.gnu.org/licenses/>.
//
// Copyright(C) Bertini2 Development Team
//
// See <http://www.gnu.org/licenses/> for a copy of the license,
// as well as COPYING.  Bertini2 is provided with permitted
// additional terms in the b2/licenses/ directory.

 
#pragma once
 
#include <vector>
#include <memory>
#include <stdexcept>
 
#include <boost/serialization/vector.hpp>
#include <boost/serialization/shared_ptr.hpp>
#include <boost/serialization/split_member.hpp>
 
#include "bertini2/num_traits.hpp"
#include "bertini2/eigen_extensions.hpp"
#include "bertini2/function_tree.hpp"
#include "bertini2/function_tree/reintern.hpp"
#include "bertini2/system/blocks/describe.hpp"
 
namespace bertini {
namespace blocks {

template <typename SystemT>
class RandomizationBlock
{
public:
    using Var = std::shared_ptr<node::Variable>;  
    using OperandPtr = std::shared_ptr<SystemT>;  

    RandomizationBlock() : num_groups_(0), precision_(DefaultPrecision()) {}

    RandomizationBlock(OperandPtr operand,
                       Mat<complex_mp> coefficients,
                       std::vector<std::vector<int>> target_multidegrees,
                       std::vector<std::vector<int>> operand_multidegrees,
                       size_t num_groups)
        : operand_(std::move(operand)),
          coefficients_highest_precision_(std::move(coefficients)),
          target_multidegrees_(std::move(target_multidegrees)),
          operand_multidegrees_(std::move(operand_multidegrees)),
          num_groups_(num_groups),
          precision_(DefaultPrecision())
    {
        BuildWorking();
    }
 
    // Memory-isolating copy (ADR-0027): the operand System is deep-copied via the
    // System copy constructor, which shares its immutable node DAG and compiled SLP Program but
    // gives it its own per-thread evaluation Memory.  Everything else (the coefficient matrices,
    // multidegrees, working buffers) is value state, copied per copy.  This lets path tracking
    // Clone a System into per-thread copies that share no mutable state.
    RandomizationBlock(RandomizationBlock const& other)
        : operand_(other.operand_ ? std::make_shared<SystemT>(*other.operand_) : nullptr),
          coefficients_highest_precision_(other.coefficients_highest_precision_),
          coefficients_working_(other.coefficients_working_),
          target_multidegrees_(other.target_multidegrees_),
          operand_multidegrees_(other.operand_multidegrees_),
          num_groups_(other.num_groups_),
          homogenized_(other.homogenized_),
          hom_vars_(other.hom_vars_),
          hom_var_index_(other.hom_var_index_),
          precision_(other.precision_)
    {}

    RandomizationBlock& operator=(RandomizationBlock const& other)
    {
        if (this != &other)
        {
            operand_ = other.operand_ ? std::make_shared<SystemT>(*other.operand_) : nullptr;
            coefficients_highest_precision_ = other.coefficients_highest_precision_;
            coefficients_working_ = other.coefficients_working_;
            target_multidegrees_ = other.target_multidegrees_;
            operand_multidegrees_ = other.operand_multidegrees_;
            num_groups_ = other.num_groups_;
            homogenized_ = other.homogenized_;
            hom_vars_ = other.hom_vars_;
            hom_var_index_ = other.hom_var_index_;
            precision_ = other.precision_;
        }
        return *this;
    }

    RandomizationBlock(RandomizationBlock&&) = default;
    RandomizationBlock& operator=(RandomizationBlock&&) = default;

    size_t NumFunctions() const { return static_cast<size_t>(coefficients_highest_precision_.rows()); }

    std::vector<int> Degrees() const
    {
        std::vector<int> d(NumFunctions(), 0);
        for (size_t i = 0; i < target_multidegrees_.size(); ++i)
            for (int e : target_multidegrees_[i])
                d[i] += e;
        return d;
    }

    std::vector<int> Degrees(VariableGroup const& vars) const
    {
        const int g = GroupIndex(vars);
        if (g < 0)
            return Degrees();
        std::vector<int> d(NumFunctions(), 0);
        for (size_t i = 0; i < target_multidegrees_.size(); ++i)
            d[i] = target_multidegrees_[i][static_cast<size_t>(g)];
        return d;
    }

    bool IsPolynomial(VariableGroup const& /*vars*/) const { return operand_->IsPolynomial(); }
    bool IsHomogeneous(VariableGroup const& /*vars*/) const { return homogenized_ ? operand_->IsHomogeneous() : false; }

    void Homogenize(VariableGroup const& /*group*/, Var const& hom_var)
    {
        if (homogenized_)
            return;
        hom_vars_.push_back(hom_var);
        if (hom_vars_.size() == num_groups_)
        {
            operand_->Homogenize(VariableGroup(hom_vars_));
            homogenized_ = true;
            // the hom var of group g now sits at a fixed position of the (shared) variable ordering;
            // the block reads h_g from there during evaluation.
            hom_var_index_.assign(num_groups_, -1);
            auto const& ordering = operand_->Variables();
            for (size_t g = 0; g < num_groups_; ++g)
                for (size_t k = 0; k < ordering.size(); ++k)
                    if (ordering[k] == hom_vars_[g]) { hom_var_index_[g] = static_cast<int>(k); break; }
        }
    }

    void RestoreHomogenization(std::vector<Var> hom_vars)
    {
        if (hom_vars.size() != num_groups_)
            throw std::invalid_argument("RandomizationBlock: " + std::to_string(hom_vars.size())
                + " homogenizing variables for " + std::to_string(num_groups_) + " variable groups");
        hom_vars_ = std::move(hom_vars);
        homogenized_ = true;
        hom_var_index_.assign(num_groups_, -1);
        auto const& ordering = operand_->Variables();
        for (size_t g = 0; g < num_groups_; ++g)
            for (size_t k = 0; k < ordering.size(); ++k)
                if (ordering[k] == hom_vars_[g]) { hom_var_index_[g] = static_cast<int>(k); break; }
    }

    bool DependsOnPathVariable() const { return operand_->HavePathVariable(); }

    bool HasConstantJacobian() const { return false; }

    void Differentiate() const {}

    OperandPtr const& Operand() const { return operand_; }
    Mat<complex_mp> const& RandomizationMatrix() const { return coefficients_highest_precision_; }
 
    // Accessors for the function-tree expansion oracle (System::NaturalFunctionsAsNodes), which
    // rebuilds g_i = sum_j c_ij * node(f_j) * prod_g h_g^{(D_{i,g}-d_{j,g})}.
    bool IsHomogenized() const { return homogenized_; }
    std::vector<Var> const& HomVars() const { return hom_vars_; }
    std::vector<std::vector<int>> const& TargetMultidegrees() const { return target_multidegrees_; }
    std::vector<std::vector<int>> const& OperandMultidegrees() const { return operand_multidegrees_; }
    size_t NumGroups() const { return num_groups_; }

    void Reintern(node::ReinternMemo& memo)
    {
        for (auto& v : hom_vars_)
            v = std::static_pointer_cast<node::Variable>(node::Reintern(v, memo));
        operand_->ReinternNodes(memo);
    }

    void Describe(std::ostream& out, size_t& row, VariableGroup const& /*vars*/, bool verbose) const
    {
        const size_t n = NumFunctions();
        const size_t N = static_cast<size_t>(coefficients_highest_precision_.cols());
        out << "  ";
        describe_detail::PrintRowLabel(out, row, n);
        out << " = R . g   (R: " << n << "x" << N << " randomization matrix)\n";
        row += n;
 
        const size_t cap = describe_detail::kTerseRowCap;
 
        auto g = operand_->NaturalFunctionsAsNodes();
        const size_t gshown = (verbose || g.size() <= cap) ? g.size() : cap;
        for (size_t j = 0; j < gshown; ++j)
            out << "      g_" << j << " = " << g[j] << "\n";
        if (gshown < g.size())
            out << "      ... (" << (g.size() - gshown) << " more; describe(verbose=True) for all)\n";
 
        // the actual randomization matrix below -- 4 significant figures in terse (the default),
        // full precision in verbose.  Shown always (not gated behind verbose) so the numbers a
        // randomize() produced are visible at a glance.
        auto const& R   = coefficients_highest_precision_;
        const int   sig = describe_detail::CoeffSig(verbose);
        const Eigen::Index rcap   = static_cast<Eigen::Index>(cap);
        const Eigen::Index rshown = (verbose || R.rows() <= rcap) ? R.rows() : rcap;
        out << "    R =\n";
        for (Eigen::Index i = 0; i < rshown; ++i)
        {
            out << "      [ ";
            for (Eigen::Index j = 0; j < R.cols(); ++j)
            {
                if (j) out << ", ";
                describe_detail::PrintCoeff(out, R(i, j), sig);
            }
            out << " ]\n";
        }
        if (rshown < R.rows())
            out << "      ... (" << (R.rows() - rshown) << " more rows; describe(verbose=True) for all)\n";
    }

    unsigned Precision() const { return precision_; }

    void Precision(unsigned new_precision) const
    {
        // short-circuit when already materialized here.  Each holder of mp values keeps its
        // own "materialized at" tag; System deliberately keeps none and simply fans out on
        // every evaluation, which is cheap precisely because of this early return (ADR-0057).
        if (precision_==new_precision)
            return;
        if (new_precision > DoublePrecision())
        {
            auto& wm = std::get<Mat<complex_mp>>(coefficients_working_);
            for (Eigen::Index r = 0; r < wm.rows(); ++r)
                for (Eigen::Index c = 0; c < wm.cols(); ++c)
                {
                    wm(r, c).precision(new_precision);
                    if (new_precision > precision_)
                        wm(r, c) = coefficients_highest_precision_(r, c);
                }
        }
        // operand_ is a System: it self-aligns from the point it is evaluated at (ADR-0057)
        precision_ = new_precision;
    }

    template <typename T>
    void EvalInPlace(Eigen::Ref<Vec<T>> result, Vec<T> const& vars, T const& /*path_value*/) const
    {
        SyncPrecision(vars);
        const size_t N = static_cast<size_t>(coefficients_highest_precision_.cols());
        const Vec<T> f = operand_->template Eval<T>(vars).head(static_cast<Eigen::Index>(N));
        const auto& C = Working<T>();
 
        result.setZero();
        for (size_t i = 0; i < NumFunctions(); ++i)
        {
            T acc = Zero<T>();
            for (size_t j = 0; j < N; ++j)
            {
                const T& c = C(static_cast<Eigen::Index>(i), static_cast<Eigen::Index>(j));
                if (IsZero(c))
                    continue;
                acc += c * f(static_cast<Eigen::Index>(j)) * HPower<T>(i, j, vars);
            }
            result(static_cast<Eigen::Index>(i)) = acc;
        }
    }

    template <typename T>
    void JacobianInPlace(Eigen::Ref<Mat<T>> J, Vec<T> const& vars, T const& /*path_value*/) const
    {
        SyncPrecision(vars);
        const size_t N = static_cast<size_t>(coefficients_highest_precision_.cols());
        const Vec<T> f  = operand_->template Eval<T>(vars).head(static_cast<Eigen::Index>(N));
        const Mat<T> Jf = operand_->template Jacobian<T>(vars).topRows(static_cast<Eigen::Index>(N));
        const auto& C = Working<T>();
 
        J.setZero();
        for (size_t i = 0; i < NumFunctions(); ++i)
        {
            for (size_t j = 0; j < N; ++j)
            {
                const T& c = C(static_cast<Eigen::Index>(i), static_cast<Eigen::Index>(j));
                if (IsZero(c))
                    continue;
 
                const T W = HPower<T>(i, j, vars);                 // W_ij = prod_g h_g^{e_g}
                // operand-Jacobian term: c_ij * W_ij * (d f_j / d x), across all columns.
                J.row(static_cast<Eigen::Index>(i)) += (c * W) * Jf.row(static_cast<Eigen::Index>(j));
 
                // h-power terms: only the homogenizing-variable columns pick up f_j * dW/dh_{g'}.
                if (homogenized_)
                {
                    for (size_t g = 0; g < num_groups_; ++g)
                    {
                        const int e = Deficit(i, j, g);
                        if (e <= 0 || hom_var_index_[g] < 0)
                            continue;                              // e<0 means c_ij==0 (already skipped); e==0 -> no dependence
                        // dW/dh_g = e * h_g^{e-1} * prod_{g2 != g} h_{g2}^{e_{g2}}
                        const T dW = MixedPartial<T>(i, j, g, vars);
                        J(static_cast<Eigen::Index>(i), hom_var_index_[g]) += c * f(static_cast<Eigen::Index>(j)) * dW;
                    }
                }
            }
        }
    }

    template <typename T>
    void TimeDerivInPlace(Eigen::Ref<Vec<T>> result, Vec<T> const& /*vars*/, T const& /*path_value*/) const
    {
        result.setZero();
    }
 
private:
    template <typename T>
    void SyncPrecision(Vec<T> const& vars) const
    {
        if constexpr (!std::is_same<T, complex_dbl>::value)
        {
            if (vars.size() > 0)
            {
                const unsigned p = bertini::Precision(vars(0));
                if (p != precision_)
                    Precision(p);
            }
        }
    }
 
    template <typename T>
    const Mat<T>& Working() const { return std::get<Mat<T>>(coefficients_working_); }

    template <typename T>
    T Zero() const
    {
        T z(0);
        if constexpr (!std::is_same<T, complex_dbl>::value) z.precision(precision_);
        return z;
    }
 
    static bool IsZero(complex_dbl const& c) { return c == complex_dbl(0); }
    static bool IsZero(complex_mp const& c) { return c.real() == 0 && c.imag() == 0; }

    int Deficit(size_t i, size_t j, size_t g) const
    {
        return target_multidegrees_[i][g] - operand_multidegrees_[j][g];
    }

    template <typename T>
    T HPower(size_t i, size_t j, Vec<T> const& vars) const
    {
        T w(1);
        if constexpr (!std::is_same<T, complex_dbl>::value) w.precision(precision_);
        if (!homogenized_)
            return w;
        for (size_t g = 0; g < num_groups_; ++g)
        {
            const int e = Deficit(i, j, g);
            if (e == 0 || hom_var_index_[g] < 0)
                continue;
            w *= IntPow<T>(vars(hom_var_index_[g]), e);
        }
        return w;
    }

    template <typename T>
    T MixedPartial(size_t i, size_t j, size_t g, Vec<T> const& vars) const
    {
        const int eg = Deficit(i, j, g);
        T d(eg);
        if constexpr (!std::is_same<T, complex_dbl>::value) d.precision(precision_);
        d *= IntPow<T>(vars(hom_var_index_[g]), eg - 1);
        for (size_t g2 = 0; g2 < num_groups_; ++g2)
        {
            if (g2 == g)
                continue;
            const int e2 = Deficit(i, j, g2);
            if (e2 == 0 || hom_var_index_[g2] < 0)
                continue;
            d *= IntPow<T>(vars(hom_var_index_[g2]), e2);
        }
        return d;
    }

    template <typename T>
    static T IntPow(T const& base, int e)
    {
        T result(1);
        if constexpr (!std::is_same<T, complex_dbl>::value) result.precision(bertini::Precision(base));
        for (int k = 0; k < e; ++k)
            result *= base;
        return result;
    }

    int GroupIndex(VariableGroup const& vars) const
    {
        auto const& groups = operand_->VariableGroups();
        for (size_t g = 0; g < groups.size(); ++g)
            if (SameGroup(groups[g], vars))
                return static_cast<int>(g);
        return -1;
    }
    static bool SameGroup(VariableGroup const& a, VariableGroup const& b)
    {
        if (a.size() != b.size())
            return false;
        for (size_t k = 0; k < a.size(); ++k)
            if (a[k] != b[k])
                return false;
        return true;
    }
 
    void BuildWorking() const
    {
        const auto& M = coefficients_highest_precision_;
        auto& wd = std::get<Mat<complex_dbl>>(coefficients_working_);
        auto& wm = std::get<Mat<complex_mp>>(coefficients_working_);
        wd.resize(M.rows(), M.cols());
        wm.resize(M.rows(), M.cols());
        for (Eigen::Index r = 0; r < M.rows(); ++r)
            for (Eigen::Index c = 0; c < M.cols(); ++c)
            {
                wd(r, c) = complex_dbl(M(r, c));
                wm(r, c) = M(r, c);
            }
    }
 
    OperandPtr operand_;
    Mat<complex_mp> coefficients_highest_precision_;             
    mutable std::tuple<Mat<complex_dbl>, Mat<complex_mp>> coefficients_working_;
    std::vector<std::vector<int>> target_multidegrees_;            
    std::vector<std::vector<int>> operand_multidegrees_;           
    size_t num_groups_;                                            
    bool homogenized_ = false;
    std::vector<Var> hom_vars_;                                    
    std::vector<int> hom_var_index_;                               
    mutable unsigned precision_;
 
    friend class boost::serialization::access;
 
    // Operand held as shared_ptr<SystemT> (non-const, since Homogenize mutates the twin); split
    // like BlendBlock so serialization never deserializes into a const object.
    template <typename Archive>
    void save(Archive& ar, const unsigned /*version*/) const
    {
        ar & operand_;
        ar & coefficients_highest_precision_;
        ar & std::get<0>(coefficients_working_);
        ar & std::get<1>(coefficients_working_);
        ar & target_multidegrees_;
        ar & operand_multidegrees_;
        ar & num_groups_;
        ar & homogenized_;
        ar & hom_vars_;
        ar & hom_var_index_;
        ar & precision_;
    }
 
    template <typename Archive>
    void load(Archive& ar, const unsigned /*version*/)
    {
        ar & operand_;
        ar & coefficients_highest_precision_;
        ar & std::get<0>(coefficients_working_);
        ar & std::get<1>(coefficients_working_);
        ar & target_multidegrees_;
        ar & operand_multidegrees_;
        ar & num_groups_;
        ar & homogenized_;
        ar & hom_vars_;
        ar & hom_var_index_;
        ar & precision_;
    }
 
    BOOST_SERIALIZATION_SPLIT_MEMBER()
};
 
} // namespace blocks
} // namespace bertini