polynomial_block.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
//polynomial_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.
//
//polynomial_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 polynomial_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 <boost/serialization/access.hpp>
#include <boost/serialization/vector.hpp>
#include <boost/serialization/shared_ptr.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/straight_line_program.hpp"
 
namespace bertini {
namespace blocks {

class PolynomialBlock
{
public:
    using NE  = std::shared_ptr<node::NamedExpression>;  
    using Nd  = std::shared_ptr<node::Node>;  
    using Var = std::shared_ptr<node::Variable>;  

    PolynomialBlock() : precision_(DefaultPrecision()) {}
 
    // ---- construction (System forwards AddFunction / AddSubFunction here) ----
    void AddFunction(Nd const& f)    { functions_.push_back(f);    Invalidate(); }
    void AddConstant(NE const& f)    { constant_subfunctions_.push_back(f); Invalidate(); }

    void SetVariableOrdering(VariableGroup const& vars) const { variables_ = vars; Invalidate(); }
    void SetPathVariable(Var const& t) const { path_variable_ = t; Invalidate(); }
    void ClearPathVariable() const { path_variable_.reset(); Invalidate(); }
 
    // Mutable access for the owning System's construction-time manipulations (Homogenize walks
    // the trees in place; Reorder/Simplify reassign entries).  The System keeps the variable
    // groups / ordering; the block keeps the functions and their derivatives.
    std::vector<Nd>&       Functions()       { return functions_; }
    std::vector<Nd> const& Functions() const { return functions_; }
    std::vector<NE> const& ConstantSubfunctions() const { return constant_subfunctions_; }
    size_t NumConstants() const { return constant_subfunctions_.size(); }

    bool IsDifferentiated() const { return is_differentiated_; }
    void Invalidate() const { is_differentiated_ = false; }

    std::vector<int> Degrees() const
    {
        std::vector<int> d; d.reserve(functions_.size());
        for (auto const& f : functions_) d.push_back(f->Degree());
        return d;
    }

    std::vector<int> Degrees(VariableGroup const& vars) const
    {
        std::vector<int> d; d.reserve(functions_.size());
        for (auto const& f : functions_) d.push_back(f->Degree(vars));
        return d;
    }

    void Homogenize(VariableGroup const& group, Var const& hom_var)
    {
        for (auto& f : functions_)
            f = f->Homogenized(group, hom_var);
        Invalidate();
    }

    bool IsHomogeneous(VariableGroup const& vars) const
    {
        for (auto const& f : functions_) if (!f->IsHomogeneous(vars)) return false;
        return true;
    }

    bool IsPolynomial(VariableGroup const& vars) const
    {
        for (auto const& f : functions_) if (!f->IsPolynomial(vars)) return false;
        return true;
    }
 
    // ---- block contract: metadata ----
    size_t NumFunctions() const { return functions_.size(); }
    bool DependsOnPathVariable() const { return static_cast<bool>(path_variable_); }

    void Describe(std::ostream& out, size_t& row, VariableGroup const& /*vars*/, bool /*verbose*/) const
    {
        for (auto const& f : functions_)
            out << "  f_" << row++ << " = " << f << "\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;
        // The SLP (its per-thread Memory) is the sole evaluator and carries its own precision; the
        // function / derivative / variable nodes are no longer evaluated during tracking, so their
        // precision is vestigial and left untouched (keeps the shared node DAG read-only).
        slp_.precision(new_precision);
        precision_ = new_precision;
    }
 
    // ---- block contract: evaluation (value-in) ----
    // The compiled SLP is the sole evaluator; the function/derivative trees survive only as the
    // thing the SLP is compiled from (and that Simplify/Differentiate operate on).
    template <typename T>
    void EvalInPlace(Eigen::Ref<Vec<T>> result, Vec<T> const& vars, T const& path_value) const
    {
        EnsureDifferentiated();
        SetValues<T>(vars, path_value);
        slp_.template GetFuncValsInPlace<T>(result);
    }

    template <typename T>
    void JacobianInPlace(Eigen::Ref<Mat<T>> J, Vec<T> const& vars, T const& path_value) const
    {
        EnsureDifferentiated();
        SetValues<T>(vars, path_value);
        slp_.template GetJacobianInPlace<T>(J);
    }

    template <typename T>
    void TimeDerivInPlace(Eigen::Ref<Vec<T>> result, Vec<T> const& vars, T const& path_value) const
    {
        if (!path_variable_) { result.setZero(); return; }
        EnsureDifferentiated();
        SetValues<T>(vars, path_value);
        slp_.template GetTimeDerivInPlace<T>(result);
    }
 
    // ---- accessors the SLP compiler reads (mirror the System names) ----
    VariableGroup const& VariableOrdering() const { return variables_; }
    bool HavePathVariable() const { return static_cast<bool>(path_variable_); }
    Var GetPathVariable() const { return path_variable_; }
    size_t NumNaturalFunctions() const { return functions_.size(); }
    std::vector<Nd> const& GetNaturalFunctions() const { return functions_; }
    // The SLP is built from the explicit per-variable derivative trees (space/time derivatives).
    std::vector<Nd> const& GetSpaceDerivatives() const
    {
        if (space_derivatives_.empty())
            DifferentiateUsingDerivatives();
        return space_derivatives_;
    }

    std::vector<Nd> const& GetTimeDerivatives() const
    {
        if (path_variable_ && time_derivatives_.empty())
            DifferentiateUsingDerivatives();
        return time_derivatives_;
    }

    void Differentiate() const
    {
        if (is_differentiated_) return;
        DifferentiateUsingDerivatives();
        is_differentiated_ = true;  // set before Compile, which reads the deriv state
        slp_ = SLPCompiler().Compile(*this);
    }

    void SimplifyFunctions() const
    {
        // Intentionally inert on the top-level functions, preserving long-standing behavior:
        // when functions were Function-wrapped, Handle::Simplified() returned self (an opaque
        // boundary), so top-level simplification never happened.  Now that functions are bare
        // roots, actually simplifying them here would be a behavior change, deferred to its own step.
        Invalidate();
    }

    void SimplifyDerivatives() const
    {
        using bertini::Simplify;
        for (auto& n : space_derivatives_) n = Simplify(n);
        for (auto& n : time_derivatives_)  n = Simplify(n);
    }

    void Reintern(node::ReinternMemo& memo)
    {
        for (auto& f : functions_)
            f = node::Reintern(f, memo);
        for (auto& c : constant_subfunctions_)
            c = std::static_pointer_cast<node::NamedExpression>(node::Reintern(c, memo));
        for (auto& v : variables_)
            v = std::static_pointer_cast<node::Variable>(node::Reintern(v, memo));
        if (path_variable_)
            path_variable_ = std::static_pointer_cast<node::Variable>(node::Reintern(path_variable_, memo));
        space_derivatives_.clear();
        time_derivatives_.clear();
        Invalidate();
    }
 
private:
    void EnsureDifferentiated() const { if (!is_differentiated_) Differentiate(); }
 
    template <typename T>
    void SetValues(Vec<T> const& vars, T const& path_value) const
    {
        slp_.SetVariableValues(vars);
        if (path_variable_) slp_.template SetPathVariable<T>(path_value);
    }
 
    void DifferentiateUsingDerivatives() const
    {
        const size_t n = functions_.size();
        space_derivatives_.resize(n * variables_.size());
        for (size_t jj = 0; jj < variables_.size(); ++jj)
            for (size_t ii = 0; ii < n; ++ii)
                space_derivatives_[ii + jj * n] = functions_[ii]->Differentiate(variables_[jj]);
        if (path_variable_)
        {
            time_derivatives_.resize(n);
            for (size_t ii = 0; ii < n; ++ii)
                time_derivatives_[ii] = functions_[ii]->Differentiate(path_variable_);
        }
    }
 
    mutable std::vector<Nd> functions_;
    mutable std::vector<NE> constant_subfunctions_;
 
    mutable std::vector<Nd>  space_derivatives_;
    mutable std::vector<Nd>  time_derivatives_;
    mutable StraightLineProgram slp_;
 
    mutable VariableGroup variables_;   
    mutable Var path_variable_;         
 
    mutable bool is_differentiated_ = false;
    mutable unsigned precision_;
 
    friend class boost::serialization::access;
 
    template <typename Archive>
    void serialize(Archive& ar, const unsigned /*version*/)
    {
        ar & constant_subfunctions_;
        ar & functions_;
        ar & is_differentiated_;
        ar & space_derivatives_;
        ar & time_derivatives_;
        ar & slp_;
        ar & variables_;
        ar & path_variable_;
        ar & precision_;
    }
};
 
} // namespace blocks
} // namespace bertini