blend_block.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
//blend_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.
//
//blend_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 blend_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 <cassert>
 
#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 BlendBlock
{
public:
    using Nd = std::shared_ptr<node::Node>;  
    using Var = std::shared_ptr<node::Variable>;  
    using OperandPtr = std::shared_ptr<const SystemT>;  

    BlendBlock() : precision_(DefaultPrecision()) {}

    BlendBlock(Var path_variable, std::vector<Nd> coefficients, std::vector<OperandPtr> operands)
        : path_variable_(std::move(path_variable)),
          coefficients_(std::move(coefficients)),
          operands_(std::move(operands)),
          precision_(DefaultPrecision())
    {
        assert(!operands_.empty() && "BlendBlock needs at least one operand");
        assert(coefficients_.size() == operands_.size() && "one coefficient per operand");
        derivative_coefficients_.reserve(coefficients_.size());
        for (auto const& c : coefficients_)
            derivative_coefficients_.push_back(c->Differentiate(path_variable_));
    }
 
    // Memory-isolating copy (ADR-0027): each 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.  The coefficient-system cache is reset (rebuilt
    // lazily per copy).  The coefficient nodes / path variable are shared --- they are never written
    // during evaluation.  This lets path tracking Clone a System into per-thread copies that share
    // no mutable state.
    BlendBlock(BlendBlock const& other)
        : path_variable_(other.path_variable_),
          coefficients_(other.coefficients_),
          derivative_coefficients_(other.derivative_coefficients_),
          precision_(other.precision_)
    {
        operands_.reserve(other.operands_.size());
        for (auto const& op : other.operands_)
            operands_.push_back(std::make_shared<const SystemT>(*op));
        // coefficient_system_ deliberately left null: rebuilt lazily per copy
    }

    BlendBlock& operator=(BlendBlock const& other)
    {
        if (this != &other)
        {
            path_variable_ = other.path_variable_;
            coefficients_ = other.coefficients_;
            derivative_coefficients_ = other.derivative_coefficients_;
            precision_ = other.precision_;
            operands_.clear();
            operands_.reserve(other.operands_.size());
            for (auto const& op : other.operands_)
                operands_.push_back(std::make_shared<const SystemT>(*op));
            coefficient_system_.reset();
        }
        return *this;
    }

    BlendBlock(BlendBlock&&) = default;
    BlendBlock& operator=(BlendBlock&&) = default;

    size_t NumFunctions() const
    {
        return operands_.empty() ? 0 : operands_.front()->NumNaturalFunctions();
    }

    std::vector<Nd> const& Coefficients() const { return coefficients_; }
    std::vector<OperandPtr> const& Operands() const { return operands_; }
    Var const& PathVariable() const { return path_variable_; }

    void Describe(std::ostream& out, size_t& row, VariableGroup const& /*vars*/, bool verbose) const
    {
        const size_t k = NumFunctions();
        out << "  ";
        describe_detail::PrintRowLabel(out, row, k);
        out << " = ";
        for (size_t i = 0; i < coefficients_.size(); ++i)
            out << (i ? " + " : "") << "(" << coefficients_[i] << ")*"
                << static_cast<char>('A' + (i < 26 ? i : 25));
        out << "   (blend of " << operands_.size() << " systems)\n";
        row += k;
        if (verbose)
            for (size_t i = 0; i < operands_.size(); ++i)
            {
                const char label = static_cast<char>('A' + (i < 26 ? i : 25));
                out << "    " << label << ":\n";
                auto f = operands_[i]->NaturalFunctionsAsNodes();
                for (size_t j = 0; j < f.size(); ++j)
                    out << "      " << label << "_" << j << " = " << f[j] << "\n";
            }
    }

    static int CombineDegrees(int a, int b)
    {
        if (a < 0 || b < 0)
            return -1;
        return std::max(a, b);
    }

    std::vector<int> Degrees() const
    {
        std::vector<int> d(NumFunctions(), 0);
        for (auto const& op : operands_)
        {
            auto od = op->Degrees();
            for (size_t i = 0; i < d.size() && i < od.size(); ++i) d[i] = CombineDegrees(d[i], od[i]);
        }
        return d;
    }

    std::vector<int> Degrees(VariableGroup const& vars) const
    {
        std::vector<int> d(NumFunctions(), 0);
        for (auto const& op : operands_)
        {
            auto od = op->Degrees(vars);
            for (size_t i = 0; i < d.size() && i < od.size(); ++i) d[i] = CombineDegrees(d[i], od[i]);
        }
        return d;
    }
 
    // A blend is the coupling homotopy, built from already-prepared (homogenized) operands; its
    // operands are shared_ptr<const System> and cannot be mutated, so Homogenize is a no-op.  It
    // is homogeneous/polynomial iff all its operands are.
    void Homogenize(VariableGroup const&, std::shared_ptr<node::Variable> const&) {}
    bool IsHomogeneous(VariableGroup const&) const
    {
        for (auto const& op : operands_) if (!op->IsHomogeneous()) return false;
        return true;
    }

    bool IsPolynomial(VariableGroup const&) const
    {
        for (auto const& op : operands_) if (!op->IsPolynomial()) return false;
        return true;
    }

    bool DependsOnPathVariable() const { return true; }

    bool HasConstantJacobian() const { return false; }

    void Differentiate() const {}

    void Reintern(node::ReinternMemo& memo)
    {
        if (path_variable_)
            path_variable_ = std::static_pointer_cast<node::Variable>(node::Reintern(path_variable_, memo));
        for (auto& c : coefficients_)
            c = node::Reintern(c, memo);
        derivative_coefficients_.clear();
        derivative_coefficients_.reserve(coefficients_.size());
        for (auto const& c : coefficients_)
            derivative_coefficients_.push_back(c->Differentiate(path_variable_));
        for (auto& op : operands_)
            std::const_pointer_cast<SystemT>(op)->ReinternNodes(memo);
        coefficient_system_.reset();  // rebuilt lazily against the reinterned coefficients
    }

    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;
        // operands are Systems: they self-align from the point they are evaluated at, so
        // nothing is pushed into them here (ADR-0057)
        // The coefficients are evaluated through the coefficient sub-system's SLP (which carries
        // its own precision; see EvalCoefficients), so the coefficient nodes and the shared path
        // variable are no longer evaluated during tracking.  Their precision is vestigial and
        // left untouched, keeping the shared node DAG read-only across threads (ADR-0027).
        precision_ = new_precision;
    }

    template <typename T>
    void EvalInPlace(Eigen::Ref<Vec<T>> result, Vec<T> const& vars, T const& path_value) const
    {
        SyncPrecision(vars);
        result.setZero();
        const Eigen::Index k = static_cast<Eigen::Index>(NumFunctions());
        const Vec<T> coeffs = EvalCoefficients<T>(path_value);
        for (size_t i = 0; i < operands_.size(); ++i)
        {
            const Vec<T> fi = operands_[i]->template Eval<T>(vars);
            result += coeffs(static_cast<Eigen::Index>(i)) * fi.head(k);
        }
    }

    template <typename T>
    void JacobianInPlace(Eigen::Ref<Mat<T>> J, Vec<T> const& vars, T const& path_value) const
    {
        SyncPrecision(vars);
        J.setZero();
        const Eigen::Index k = static_cast<Eigen::Index>(NumFunctions());
        const Vec<T> coeffs = EvalCoefficients<T>(path_value);
        for (size_t i = 0; i < operands_.size(); ++i)
        {
            const Mat<T> Ji = operands_[i]->template Jacobian<T>(vars);
            J += coeffs(static_cast<Eigen::Index>(i)) * Ji.topRows(k);
        }
    }

    template <typename T>
    void TimeDerivInPlace(Eigen::Ref<Vec<T>> result, Vec<T> const& vars, T const& path_value) const
    {
        SyncPrecision(vars);
        result.setZero();
        const Eigen::Index k = static_cast<Eigen::Index>(NumFunctions());
        const Vec<T> coeffs = EvalCoefficients<T>(path_value);
        const size_t n = operands_.size();
        for (size_t i = 0; i < n; ++i)
        {
            const Vec<T> fi = operands_[i]->template Eval<T>(vars);
            result += coeffs(static_cast<Eigen::Index>(n + i)) * fi.head(k);
        }
    }
 
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);
            }
        }
    }

    SystemT& EnsureCoefficientSystem() const
    {
        if (!coefficient_system_)
        {
            auto sys = std::make_shared<SystemT>();
            for (auto const& c : coefficients_)
                sys->AddFunction(c);
            for (auto const& c : derivative_coefficients_)
                sys->AddFunction(c);
            sys->AddVariableGroup(VariableGroup{path_variable_});
            coefficient_system_ = sys;
        }
        return *coefficient_system_;
    }

    template <typename T>
    Vec<T> EvalCoefficients(T const& path_value) const
    {
        // cs is a System: it self-aligns from t_point below, so nothing is pushed into it here
        auto& cs = EnsureCoefficientSystem();
        Vec<T> t_point(1);
        t_point(0) = path_value;
        return cs.template Eval<T>(t_point);
    }
 
    Var path_variable_;
    std::vector<Nd> coefficients_;
    std::vector<Nd> derivative_coefficients_;
    std::vector<OperandPtr> operands_;
    mutable unsigned precision_;

    mutable std::shared_ptr<SystemT> coefficient_system_;
 
    friend class boost::serialization::access;
 
    // Operands are shared_ptr<const SystemT>; serialize them as non-const pointers so we
    // never ask boost to deserialize into a const object (const-ness is a local concern).
    template <typename Archive>
    void save(Archive& ar, const unsigned /*version*/) const
    {
        ar & precision_;
        ar & path_variable_;
        ar & coefficients_;
        ar & derivative_coefficients_;
        std::vector<std::shared_ptr<SystemT>> ops;
        ops.reserve(operands_.size());
        for (auto const& o : operands_)
            ops.push_back(std::const_pointer_cast<SystemT>(o));
        ar & ops;
    }
 
    template <typename Archive>
    void load(Archive& ar, const unsigned /*version*/)
    {
        ar & precision_;
        ar & path_variable_;
        ar & coefficients_;
        ar & derivative_coefficients_;
        std::vector<std::shared_ptr<SystemT>> ops;
        ar & ops;
        operands_.assign(ops.begin(), ops.end());
    }
 
    BOOST_SERIALIZATION_SPLIT_MEMBER()
};
 
} // namespace blocks
} // namespace bertini