system.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
//system.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.
//
//system.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 system.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.

 
#ifndef BERTINI_SYSTEM_HPP
#define BERTINI_SYSTEM_HPP
 
#include <assert.h>
#include <optional>
#include <vector>
#include <set>
#include <string>
 
 
#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>
#include <boost/serialization/export.hpp>
#include <boost/serialization/shared_ptr.hpp>
#include <boost/serialization/vector.hpp>
#include <boost/serialization/deque.hpp>
#include <boost/serialization/std_variant.hpp>
#include <boost/type_index.hpp>
 
#include "bertini2/mpfr_complex.hpp"
#include "bertini2/mpfr_extensions.hpp"
#include "bertini2/eigen_extensions.hpp"
 
 
#include "bertini2/detail/sha256.hpp"
#include "bertini2/function_tree.hpp"
#include "bertini2/function_tree/reintern.hpp"
#include "bertini2/function_tree/canonical_decoding.hpp"
#include "bertini2/system/patch.hpp"
 
#include "bertini2/system/straight_line_program.hpp"
#include "bertini2/system/blocks/block.hpp"
 
#include <boost/archive/binary_oarchive.hpp>
#include <boost/archive/binary_iarchive.hpp>
#include <boost/iostreams/stream_buffer.hpp>
#include <boost/iostreams/stream.hpp>
#include <boost/iostreams/device/back_inserter.hpp>
 
namespace bertini {
 
    // (included above) so the evaluation blocks can see them.
 
    namespace node { struct EncodingContext; }  // function_tree/canonical_encoding.hpp (ADR-0042)
 
    class Slice;  // system/slice.hpp -- a thin wrapper over a LinearFormsBlock; System::Slices()
                  // recovers them from a slice-derived system without binding the block variant.

    struct NodeMatrix
    {
        std::vector<std::shared_ptr<node::Node>> entries; 
        size_t rows = 0; 
        size_t cols = 0; 
    };

    NodeMatrix SymbolicJacobian(std::vector<std::shared_ptr<node::Node>> const& functions,
                                VariableGroup const& variables);

    class System{
 
    public:
        // a few local using statements to reduce typing etc.
        using NE = std::shared_ptr<node::NamedExpression>;  
        using Var = std::shared_ptr<node::Variable>;  
        using Nd = std::shared_ptr<node::Node>;  

        System() : have_path_variable_(false), is_patched_(false), is_differentiated_(false), have_ordering_(false)
        {}

        explicit
        System(std::vector<Nd> const& functions);

        explicit
        System(std::string const& input);

        System(System const& other);

        System(System && other) : System()
        {
            swap(*this, other);
        }

        System& operator=(const System & other);

        System& operator=(System && other) = default;

        friend void swap(System & a, System & b);
 
        // There is deliberately NO way to set a precision on a System, and nothing to fan out.
        // Every evaluable type -- each block, the patch, the SLP -- self-aligns to the precision
        // of the point it is handed, under SyncPrecision.  See ADR-0057.
 
        // There is deliberately NO precision() accessor.  A System carries no precision of its
        // own: evaluation happens at the precision of the point it is handed, and the objects
        // that actually hold multiprecision values -- each block's working coefficients, the
        // patch, the SLP's memory -- each keep their own "materialized at" tag and short-circuit
        // when already there.  A System-level copy was a duplicate of those, and exposing it
        // invited callers to read it, compute a target, and hand it back, which is exactly the
        // hand-alignment ADR-0057 removes.

        void Differentiate() const;

        NodeMatrix SymbolicJacobian(bool usercoordinates = true) const;
 

        template<typename T>
        void EvalInPlace(Vec<T> & function_values) const
        {
 
            if (function_values.size() != static_cast<Eigen::Index>(NumTotalFunctions()))
            {
                std::stringstream ss;
                ss << "trying to evaluate system in-place, but number length of vector into which to write the values (" << function_values.size() << ") doesn't match number of system user-defined functions plus patches ( " << NumNaturalFunctions() << "+" << NumPatches() << ") = " << NumTotalFunctions() << ").  Use System.NumTotalFunctions() to make the container for in-place evaluation";
                throw std::runtime_error(ss.str());
            }
 
            if (!is_differentiated_)
                Differentiate();   // syncs + (for the polynomial block) compiles the SLP
 
            EvalBlocksInPlace<T>(function_values);
            if (IsPatched())
                patch_.EvalInPlace(function_values,
                                   std::get<Vec<T> >(current_variable_values_));
            CoerceBlockOutputPrecision(function_values);
        }
 
 
 
 

        template<typename T>
        Vec<T> Eval() const
        {
            Vec<T> function_values(NumTotalFunctions()); // create vector with correct number of entries.
            EvalInPlace(function_values);
 
            return function_values;
        }
 
 
 

        template<typename T, typename Derived>
        void EvalInPlace(Vec<T>& function_values, const Eigen::MatrixBase<Derived>& variable_values) const
        {
            static_assert(std::is_same<typename Derived::Scalar,T>::value,"scalar types must match");
 
            if (variable_values.size()!=static_cast<Eigen::Index>(NumVariables()))
            {
                std::stringstream ss;
                ss << "trying to evaluate system, but number of input variables (" << variable_values.size() << ") doesn't match number of system variables (" << NumVariables() << ").";
                throw std::runtime_error(ss.str());
            }
            if (have_path_variable_)
                throw std::runtime_error("not using a time value for evaluation of system, but path variable IS defined.");
 
            SetVariables(variable_values.eval());
            EvalInPlace(function_values);
        }
 
 
 

        template<typename Derived>
        typename Derived::PlainObject Eval(const Eigen::MatrixBase<Derived>& variable_values) const
        {
            typedef typename Derived::Scalar T;
 
            Vec<T> function_values(NumTotalFunctions()); // create vector with correct number of entries.
            EvalInPlace(function_values, variable_values);
            return function_values;
 
        }

        template<typename T>
        Vec<T> Eval(const Vec<T> & variable_values) const
        {
            Vec<T> function_values(NumTotalFunctions()); // create vector with correct number of entries.
            EvalInPlace(function_values, variable_values);
            return function_values;
        }
 
 
 
 

        template<typename Derived, typename T>
        void EvalInPlace(Vec<T> & function_values, const Eigen::MatrixBase<Derived>& variable_values, const T & path_variable_value) const
        {
            static_assert(std::is_same<typename Derived::Scalar, T>::value, "scalar types must be the same");
 
            if (variable_values.size()!=static_cast<Eigen::Index>(NumVariables()))
                throw std::runtime_error("trying to evaluate system, but number of variables doesn't match.");
            if (!have_path_variable_)
                throw std::runtime_error("trying to use a time value for evaluation of system, but no path variable defined.");
 
            SetVariables(variable_values.eval());
            SetPathVariable(path_variable_value);
 
 
            EvalInPlace(function_values);
        }
 
 
 
 
 
 
 
 

        template<typename Derived, typename T>
        Vec<T> Eval(const Eigen::MatrixBase<Derived>& variable_values, const T & path_variable_value) const
        {
            Vec<T> function_values(NumTotalFunctions()); // create vector with correct number of entries.
            EvalInPlace(function_values, variable_values, path_variable_value);
            return function_values;
        }
 

        template<typename T>
        Vec<T> Eval(const Vec<T> & variable_values, const T & path_variable_value) const
        {
            Vec<T> function_values(NumTotalFunctions()); // create vector with correct number of entries.
            EvalInPlace(function_values, variable_values, path_variable_value);
            return function_values;
        }
 
 

        template <typename T>
        void JacobianInPlace(Mat<T> & J) const
        {
 
 
            if(J.rows() != static_cast<Eigen::Index>(NumTotalFunctions()) || J.cols() != static_cast<Eigen::Index>(NumVariables()))
            {
                throw std::runtime_error("trying to evaluate jacobian of system in place, but input J doesn't have right number of columns or rows");
            }
 
            if (!is_differentiated_)
                Differentiate();
 
            JacobianBlocksInPlace<T>(J);
            if (IsPatched())
                patch_.JacobianInPlace(J, std::get<Vec<T> >(current_variable_values_));
            CoerceBlockOutputPrecision(J);
        }
 
 
 
 
 
 
 

        template<typename T>
        Mat<T> Jacobian() const
        {
 
            Mat<T> J(NumTotalFunctions(), NumVariables());
            JacobianInPlace(J);
 
            return J;
        }
 
 
 

        template<typename T>
        void JacobianInPlace(Mat<T> & J, const Vec<T> &  variable_values) const
        {
 
            if (variable_values.size()!=static_cast<Eigen::Index>(NumVariables()))
                throw std::runtime_error("trying to evaluate jacobian, but number of variables doesn't match.");
 
            if (HavePathVariable())
                throw std::runtime_error("not using a time value for computation of jacobian, but a path variable is defined.");
 
            SetAndReset(variable_values);
 
            JacobianInPlace(J);
        }
 
 
 
 
 

        template<typename T>
        Mat<T> Jacobian(const Vec<T> & variable_values) const
        {
            if (variable_values.size()!=static_cast<Eigen::Index>(NumVariables()))
                throw std::runtime_error("trying to evaluate jacobian, but number of variables doesn't match.");
 
            if (HavePathVariable())
                throw std::runtime_error("not using a time value for computation of jacobian, but a path variable is defined.");
 
            Mat<T> J(NumTotalFunctions(), NumVariables());
            JacobianInPlace(J,variable_values);
            return J;
        }
 
 
 

        template<typename Derived, typename T>
        void JacobianInPlace(Mat<T> & J, const Eigen::MatrixBase<Derived> & variable_values, const T & path_variable_value) const
        {
            static_assert(std::is_same<typename Derived::Scalar, T>::value, "scalar types must be the same");
 
            if (variable_values.size()!=static_cast<Eigen::Index>(NumVariables()))
                throw std::runtime_error("trying to evaluate jacobian, but number of variables doesn't match.");
 
            if (!HavePathVariable())
                throw std::runtime_error("trying to use a time value for computation of jacobian, but no path variable defined.");
 
            SetVariables(variable_values.eval());
            SetPathVariable(path_variable_value);
            JacobianInPlace(J);
        }
 
 
 
 

        template<typename Derived, typename T>
        Mat<T> Jacobian(const Eigen::MatrixBase<Derived> & variable_values, const T & path_variable_value) const
        {
            static_assert(std::is_same<typename Derived::Scalar, T>::value, "scalar types must be the same");
 
            if (variable_values.size()!=NumVariables())
                throw std::runtime_error("trying to evaluate jacobian, but number of variables doesn't match.");
 
            if (!HavePathVariable())
                throw std::runtime_error("trying to use a time value for computation of jacobian, but no path variable defined.");
 
            Mat<T> J(NumTotalFunctions(), NumVariables());
            JacobianInPlace(J,variable_values, path_variable_value);
            return J;
        }
 

        template<typename T>
        Mat<T> Jacobian(const Vec<T> & variable_values, const T & path_variable_value) const
        {
            if (variable_values.size()!=static_cast<Eigen::Index>(NumVariables()))
                throw std::runtime_error("trying to evaluate jacobian, but number of variables doesn't match.");
 
            if (!HavePathVariable())
                throw std::runtime_error("using a time value for computation of jacobian, but no path variable is defined.");
 
            Mat<T> J(NumTotalFunctions(), NumVariables());
            JacobianInPlace(J,variable_values,path_variable_value);
            return J;
        }
 

        template<typename Derived, typename T>
        void TimeDerivativeInPlace(Vec<T> & ds_dt,
                            const Eigen::MatrixBase<Derived> & variable_values,
                            const T & path_variable_value) const
        {
            static_assert(std::is_same<typename Derived::Scalar, T>::value, "scalar types must be the same");
 
            SetVariables(variable_values.eval());
            SetPathVariable(path_variable_value);
            TimeDerivativeInPlace(ds_dt);
        }
 
 
 
 
 

        template<typename Derived, typename T>
        Vec<T> TimeDerivative(const Eigen::MatrixBase<Derived> & variable_values, const T & path_variable_value) const
        {
            static_assert(std::is_same<typename Derived::Scalar, T>::value, "scalar types must be the same");
 
            Vec<T> ds_dt(NumTotalFunctions());
            TimeDerivativeInPlace(ds_dt, variable_values, path_variable_value);
            return ds_dt;
        }
 
 
 
 
 
 
 

        template<typename Derived, typename T>
        void TimeDerivativeInPlace(Vec<T> & ds_dt,
                            const Eigen::MatrixBase<Derived> & variable_values) const
        {
            static_assert(std::is_same<typename Derived::Scalar, T>::value, "scalar types must be the same");
 
            SetVariables(variable_values.eval());
            TimeDerivativeInPlace(ds_dt);
        }
 
 
 
 
 

        template<typename Derived, typename T>
        Vec<T> TimeDerivative(const Eigen::MatrixBase<Derived> & variable_values) const
        {
            static_assert(std::is_same<typename Derived::Scalar, T>::value, "scalar types must be the same");
 
            Vec<T> ds_dt(NumTotalFunctions());
            TimeDerivativeInPlace(ds_dt, variable_values);
            return ds_dt;
        }
 

        template<typename T>
        void TimeDerivativeInPlace(Vec<T> & ds_dt) const
        {
 
            if(ds_dt.size() < static_cast<Eigen::Index>(NumNaturalFunctions()))
            {
                std::stringstream ss;
                ss << "trying to evaluate system in place, but number of input functions (" << ds_dt.size() << ") doesn't match number of system functions (" << NumNaturalFunctions() << ").";
                throw std::runtime_error(ss.str());
            }
 
            if (!HavePathVariable())
                throw std::runtime_error("computing time derivative of system with no path variable defined");
 
            if (!is_differentiated_)
                Differentiate();
 
            TimeDerivBlocksInPlace<T>(ds_dt);
            // the patch doesn't move with time.  derivatives 0.
            if (IsPatched())
                for (size_t ii = 0; ii < NumTotalVariableGroups(); ++ii)
                    ds_dt(static_cast<Eigen::Index>(ii + NumNaturalFunctions())) = T(0);
            CoerceBlockOutputPrecision(ds_dt);
        }
 
 
 
 
 

        template<typename T>
        Vec<T> TimeDerivative() const
        {
            if (!HavePathVariable())
                throw std::runtime_error("computing time derivative of system with no path variable defined");
 
            Vec<T> ds_dt(NumTotalFunctions());
            TimeDerivativeInPlace(ds_dt);
            return ds_dt;
        }

        void Homogenize();

        void Homogenize(VariableGroup const& provided_hom_vars);

        bool IsHomogeneous() const;

        bool IsPolynomial() const;
 

        //
        //  Nummers --   functions which get the numbers of things.
        //
 

        size_t NumNaturalFunctions() const;

        size_t NumPatches() const
        {
            return patch_.NumVariableGroups();
        }

        size_t NumTotalFunctions() const;

        size_t NumVariables() const;

        size_t NumNaturalVariables() const;

        size_t NumHomVariables() const;

        VariableGroup const& HomogenizingVariables() const { return homogenizing_variables_; }

        VariableGroup const& UngroupedVariables() const { return ungrouped_variables_; }

        size_t NumTotalVariableGroups() const;

         size_t NumVariableGroups() const;
 

         size_t NumUngroupedVariables() const;

         size_t NumHomVariableGroups() const;
 

        size_t NumConstants() const;

        size_t NumParameters() const;
 

        size_t NumImplicitParameters() const;
 
 
 
 

        //
        // Setters -- templated.
        //

        template<typename T>
        void SetVariables(const Vec<T> & new_values) const
        {
            if (new_values.size()!= static_cast<Eigen::Index>(NumVariables()))
            {
                throw std::runtime_error("variable vector of different length from system-owned variables in SetVariables");
            }
 
            // ALIGN to the incoming point rather than refusing it (#377).  A System's
            // precision is not an invariant the caller must maintain by hand: the caller's
            // intent is always just "evaluate f at x", and x carries the precision that
            // evaluation should happen at.  Aligning here pushes the new precision down to
            // every block (each block re-tags its SLP's memory and rebuilds its constants
            // from their exact recipes), to the stored path value, and to the patch, so the
            // whole System is consistent for the evaluation that follows.
            //
            // Refusing was not merely unergonomic, it was inescapable in one direction: a
            // block's SLP takes its memory precision from the ambient DefaultPrecision() at
            // lazy compile time while the System keeps whatever it was told, and once those
            // diverge, System::precision(n) cannot repair it -- both it and
            // StraightLineProgram::precision(n) short-circuit when handed the value they
            // already hold.  See the SLP-side note in SetVariableValues.
            //
            // Deliberately NOT behind BERTINI_DISABLE_PRECISION_CHECKS: this is required
            // behaviour, not a debug assertion.  A system with no variables (a constant) has
            // an empty point and so carries no precision to read.
            // NO fan-out.  Every evaluable type self-aligns from the point it is handed --
            // each block, the patch and the SLP all do it under SyncPrecision -- so there is
            // nothing to push down from here and nothing to keep in sync (ADR-0057).
 
            // Blocks are value-in: the polynomial block feeds this stored vector into its SLP, the
            // structured blocks compute on it directly, and the patch reads it too (see
            // EvalBlocksInPlace).  The shared Variable nodes are no longer written during evaluation
            // (ADR-0027), so the node DAG stays read-only across threads.
            std::get<Vec<T> >(current_variable_values_) = new_values;
        }
 
 

        template<typename T>
        void SetPathVariable(T const& new_value) const
        {
            if (!have_path_variable_)
                throw std::runtime_error("trying to set the value of the path variable, but one is not defined for this system");
 
            // Store the path value in the System's own per-thread buffer, NOT the shared node.
            // Blocks are value-in and receive the path value as an argument (see EvalBlocksInPlace /
            // CurrentPathValue), so the path-variable node is never read during evaluation (ADR-0027).
            std::get<T>(current_path_value_) = new_value;
        }
 

        template<typename T>
        void SetAndReset(Vec<T> const& new_space, T const& new_time) const
        {
            SetVariables(new_space);
            SetPathVariable(new_time);
        }

        template<typename T>
        void SetAndReset(Vec<T> const& new_space) const
        {
            SetVariables(new_space);
        }
 
 
 
 
 
 
 
 

        //
        //  Adders  --   functions which add things to the system.
        //
 

        void AddVariableGroup(VariableGroup const& v);
 

        void SetVariableGroups(std::vector<VariableGroup> const& groups);
 

        void AddHomVariableGroup(VariableGroup const& v);
 
 

        void AddUngroupedVariable(Var const& v);
 

        void AddUngroupedVariables(VariableGroup const& v);
 

        void AddImplicitParameter(Var const& v);
 

        void AddImplicitParameters(VariableGroup const& v);
 
 
 

        void AddParameter(NE const& F);
 
 

        void AddFunction(Nd const& F);
 

        void AddFunctions(std::vector<Nd> const& F);

        void AddBlock(Block b) { ThrowIfSealed("AddBlock"); blocks_.push_back(std::move(b)); }

        std::vector<Block> const& Blocks() const { return blocks_; }

        std::vector<Slice> Slices() const;

        void ClearFunctions()
        {
            ThrowIfSealed("ClearFunctions");
            for (auto it = blocks_.begin(); it != blocks_.end(); )
            {
                if (std::holds_alternative<blocks::PolynomialBlock>(*it))
                    it = blocks_.erase(it);
                else
                    ++it;
            }
            InvalidateDifferentiation();
        }

        bool HasBlocks() const { return !blocks_.empty(); }

        bool HasStructuredBlocks() const
        {
            for (auto const& b : blocks_)
                if (!std::holds_alternative<blocks::PolynomialBlock>(b))
                    return true;
            return false;
        }

        void ClearBlocks() { ThrowIfSealed("ClearBlocks"); blocks_.clear(); }
 

        //
        //  content identity (ADR-0042)
        //

        std::string CanonicalEncodingText() const;

        static System FromCanonicalEncoding(std::string const& text);

        detail::Digest256 ContentDigest() const;

        std::size_t Hash() const;

        bool IsSame(System const& other) const;

        void Seal();

        bool IsSealed() const;

        void ReinternNodes(node::ReinternMemo& memo);

        System ExpandToFunctionTree() const;

        std::vector<Nd> NaturalFunctionsAsNodes() const;
 

        System Randomize() const;

        System Randomize(int codimension) const;

        System Randomize(Mat<complex_mp> const& R) const;

        Mat<complex_mp> RandomizationMatrix() const;
 
 
 
 
 
 

        void AddConstant(NE const& C);
 
 
 

        void AddPathVariable(Var const& v);
 

        bool HavePathVariable() const;

        const Var& GetPathVariable() const;

        std::set<std::string> VariableNameSet() const;

        VariableGroup VariableOrdering() const;
 

        const VariableGroup& Variables() const;

        VariableGroup const& AffineVariableGroup(size_t index) const
        {
            return variable_groups_[index];
        }

        std::vector<unsigned> VariableGroupSizes() const
        {
            return VariableGroupSizesFIFO();
        }

        template<typename T>
        Vec<T> DehomogenizePoint(Vec<T> const& x) const
            {
 
                if (x.size()!=static_cast<Eigen::Index>(NumVariables())){
                    std::stringstream message;
                    message << "dehomogenizing point with incorrect number of coordinates. input has ";
                    message << x.size();
                    message << " but system expects ";
                    message << NumVariables();
                    throw std::runtime_error(message.str());
                }
 
                if (!have_ordering_)
                    ConstructOrdering();
 
                return DehomogenizePointFIFO(x);
            }
 

        template<typename T>
        auto InfinityNormOfDehomogenized(Vec<T> const& x) const
        {
            auto const d = DehomogenizePoint(x);
            if (!HaveAuxiliaryCoordinates())
                return d.template lpNorm<Eigen::Infinity>();
 
            // the value type, not decltype(abs(...)): a multiprecision abs() is an expression
            // template, and a variable of that type cannot be reassigned in the loop
            using std::abs;
            using RealT = typename NumTraits<T>::Real;
            RealT biggest{0};
            for (Eigen::Index k{0}; k < d.size(); ++k)
                if (!CoordinateIsAuxiliary(static_cast<unsigned>(k)))
                {
                    RealT const a = abs(d(k));
                    if (a > biggest) biggest = a;
                }
            return biggest;
        }
 

        template<typename T, typename R>
        bool IsFinite(Vec<T> const& x, R const& threshold) const
        {
            return InfinityNormOfDehomogenized(x) <= threshold;
        }
 

        template<typename T, typename R>
        bool IsReal(Vec<T> const& x, R const& threshold) const
        {
            using std::abs; using std::imag;
            using RealT = typename NumTraits<T>::Real;
            auto const d = DehomogenizePoint(x);
            RealT biggest{0};
            for (Eigen::Index k{0}; k < d.size(); ++k)
                if (!CoordinateIsAuxiliary(static_cast<unsigned>(k)))
                {
                    RealT const a = abs(imag(d(k)));
                    if (a > biggest) biggest = a;
                }
            return biggest <= threshold;
        }
 

        void SetAuxiliaryVariableGroups(std::vector<unsigned> groups);
 

        void SetAuxiliaryCoordinates(std::vector<unsigned> coordinates);
 

        std::vector<unsigned> const& AuxiliaryVariableGroups() const
        {
            return auxiliary_variable_groups_;
        }

        std::vector<unsigned> const& AuxiliaryCoordinates() const
        {
            return auxiliary_coordinates_;
        }

        bool HaveAuxiliaryCoordinates() const
        {
            return !auxiliary_variable_groups_.empty() || !auxiliary_coordinates_.empty();
        }

        bool CoordinateIsAuxiliary(unsigned index) const;

        size_t NumAuxiliaryCoordinates() const;
 

        template<typename T>
        Vec<T> HomogenizePoint(Vec<T> const& x) const
            {
 
                if (x.size()!=static_cast<Eigen::Index>(NumNaturalVariables())){
                    std::stringstream message;
                    message << "homogenizing point with incorrect number of coordinates. input has ";
                    message << x.size();
                    message << " but system expects ";
                    message << NumNaturalVariables();
                    throw std::runtime_error(message.str());
                }
 
                if (!have_ordering_)
                    ConstructOrdering();
 
                auto x_homogenized = HomogenizePointFIFO(x);
 
                if (IsPatched())
                    RescalePointToFitPatchInPlace(x_homogenized);
 
                return x_homogenized;
            }
 
 
 

        auto Function(unsigned index) const
        {
            auto fns = GetFunctions();
            if (index >= fns.size())
                throw std::out_of_range("System::Function index out of range");
            return fns[index];
        }
 

        std::vector<Nd> GetNaturalFunctions() const
        {
            return NaturalFunctionsAsNodes();
        }
 

        std::vector<Nd> GetFunctions() const
        {
            auto fns = NaturalFunctionsAsNodes();
            auto patches = PatchFunctionsAsNodes();
            fns.insert(fns.end(), patches.begin(), patches.end());
            return fns;
        }
 

        std::vector<Nd> PatchFunctionsAsNodes() const;
 

        void CopyFunctions(System const& other)
        {
            AddFunctions(other.GetNaturalFunctions());
        }
 

        unsigned FIFOIndexOfGroup(VariableGroup const& g) const
        {
            unsigned affine_counter = 0, hom_counter = 0, fifo_i = 0;
            for (auto const& iter : time_order_of_variable_groups_)
            {
                bool match = false;
                switch (iter)
                {
                    case VariableGroupType::Affine:      match = (variable_groups_[affine_counter++] == g); break;
                    case VariableGroupType::Homogeneous: match = (hom_variable_groups_[hom_counter++] == g); break;
                    case VariableGroupType::Ungrouped:   break;
                    default: throw std::runtime_error("unacceptable VariableGroupType in FIFOIndexOfGroup");
                }
                if (match)
                    return fifo_i;
                ++fifo_i;
            }
            throw std::runtime_error("FIFOIndexOfGroup: the given variable group is not a group of this system");
        }
 

        template<typename T>
        Vec<T> CoordinatesOfGroup(Vec<T> const& user_point, unsigned group_fifo_index) const
        {
            auto const span = SpanOfGroup(group_fifo_index);
            if (span.first + span.second > static_cast<unsigned>(user_point.size()))
                throw std::runtime_error("CoordinatesOfGroup: point is shorter than the system's variable structure implies");
            return user_point.segment(span.first, span.second);
        }
 

        std::pair<unsigned, unsigned> SpanOfGroup(unsigned group_fifo_index) const
        {
            unsigned affine_counter = 0, hom_counter = 0, fifo_i = 0, offset = 0;
            for (auto const& iter : time_order_of_variable_groups_)
            {
                unsigned group_size;
                switch (iter)
                {
                    case VariableGroupType::Affine:      group_size = static_cast<unsigned>(variable_groups_[affine_counter++].size()); break;
                    case VariableGroupType::Homogeneous: group_size = static_cast<unsigned>(hom_variable_groups_[hom_counter++].size()); break;
                    case VariableGroupType::Ungrouped:   group_size = 1; break;
                    default: throw std::runtime_error("unacceptable VariableGroupType in SpanOfGroup");
                }
                if (fifo_i == group_fifo_index)
                    return {offset, group_size};
                offset += group_size;
                ++fifo_i;
            }
            throw std::out_of_range("SpanOfGroup: group index out of range");
        }
 
 

        auto VariableGroups() const
        {
            return variable_groups_;
        }

        auto HomVariableGroups() const
        {
            return hom_variable_groups_;
        }
 

        const Var& PathVariable() const
        {
            return path_variable_;
        }
 

        std::vector< Nd > GetSpaceDerivatives() const;

        std::vector< Nd > GetTimeDerivatives() const;
 

        //
        //  Functions involving coefficients of the system
        //
 

        template <typename NumT>
        typename Eigen::NumTraits<NumT>::Real CoefficientBound(unsigned num_evaluations = 1) const;
 
 

        int DegreeBound() const;

         std::vector<int> Degrees() const;

         std::vector<int> Degrees(VariableGroup const& vars) const;

        void ReorderFunctionsByDegreeDecreasing();
 

        void ReorderFunctionsByDegreeIncreasing();
 
 
 
 

        //
        //  Functions regarding patches.
        //
 
 
 
 

        void AutoPatch()
        {
            AutoPatchFIFO();
        }

        void CopyPatches(System const& other);
 

        Patch GetPatch() const
        {
            return patch_;
        }

        bool IsPatched() const
        {
            return is_patched_;
        }
 

        template <typename T>
        Vec<T> RescalePointToFitPatch(Vec<T> const& x) const
        {
            return patch_.RescalePoint(x);
        }
 

        template<typename T>
        void RescalePointToFitPatchInPlace(Vec<T> & x) const
        {
            patch_.RescalePointToFitInPlace(x);
        }

        void Describe(std::ostream& out, bool verbose = false) const;

        friend std::ostream& operator <<(std::ostream& out, const System & s);
 

        void ClearVariables();
 
 
 
 
 

        void CopyVariableStructure(System const& other);
 
        // The Please/Dont AssumeUniformPrecision family was removed: the setter had
        // ignored its argument (always storing false) for ages, so the early-out it was
        // meant to enable was dead code, and skipping the propagation is unsound anyway.
        // A System no longer carries a precision at all, and there is nothing to fan out
        // (ADR-0057): every evaluable type self-aligns to the point it is handed.

        void SimplifyFunctions();

        void SimplifyDerivatives() const;

        void Simplify();
 
 
 
 

        System& operator+=(System const& rhs);

        friend const System operator+(System lhs, System const& rhs);

        System& operator*=(Nd const& N);

        friend const System operator*(System s, Nd const&  N);

        friend const System operator*(Nd const&  N, System const& s);
    private:

        System AssembleRandomized(std::shared_ptr<System> operand, Mat<complex_mp> coefficients) const;

        Mat<complex_mp> BuildRandomizationMatrix(std::shared_ptr<System> operand, std::size_t codimension) const;

        std::vector<unsigned> VariableGroupSizesFIFO() const;
 

        void AutoPatchFIFO();
 
 

        template<typename T>
        Vec<T> DehomogenizePointFIFO(Vec<T> const& x) const
        {
            #ifndef BERTINI_DISABLE_ASSERTS
            assert(homogenizing_variables_.size()==0 || homogenizing_variables_.size()==NumVariableGroups() && "must have either 0 homogenizing variables, or the number of homogenizing variables must match the number of affine variable groups.");
            #endif
 
            bool is_homogenized = homogenizing_variables_.size()!=0;
            Vec<T> x_dehomogenized(NumNaturalVariables());
 
            unsigned affine_group_counter = 0;
            unsigned hom_group_counter = 0;
 
            unsigned hom_index = 0; // index into x, the point we are dehomogenizing
            unsigned dehom_index = 0; // index into x_dehomogenized, the point we are computing
 
            for (auto& iter : time_order_of_variable_groups_)
            {
                switch (iter){
                    case VariableGroupType::Affine:
                    {
                        if (is_homogenized)
                        {
                            auto h = x(hom_index++);
                            for (unsigned ii = 0; ii < variable_groups_[affine_group_counter].size(); ++ii)
                                x_dehomogenized(dehom_index++) = x(hom_index++) / h;
                            affine_group_counter++;
                        }
                        else
                        {
                            for (unsigned ii = 0; ii < variable_groups_[affine_group_counter].size(); ++ii)
                                x_dehomogenized(dehom_index++) = x(hom_index++);
                        }
                        break;
                    }
                    case VariableGroupType::Homogeneous:
                    {
                        for (unsigned ii = 0; ii < hom_variable_groups_[hom_group_counter].size(); ++ii)
                            x_dehomogenized(dehom_index++) = x(hom_index++);
                        hom_group_counter++; // was missing; mattered only for multiple hom groups of differing sizes
                        break;
                    }
                    case VariableGroupType::Ungrouped:
                    {
                        x_dehomogenized(dehom_index++) = x(hom_index++);
                        break;
                    }
                    default:
                    {
                        throw std::runtime_error("unacceptable VariableGroupType in FIFOVariableOrdering");
                    }
                }
            }
 
            return x_dehomogenized;
        }
 

        template<typename T>
        Vec<T> HomogenizePointFIFO(Vec<T> const& x) const
        {
            #ifndef BERTINI_DISABLE_ASSERTS
            assert(homogenizing_variables_.size()==0 || homogenizing_variables_.size()==NumVariableGroups() && "must have either 0 homogenizing variables, or the number of homogenizing variables must match the number of affine variable groups.");
            #endif
 
            bool is_homogenized = homogenizing_variables_.size()!=0;
            if (!is_homogenized)
                return x;
 
            Vec<T> x_homogenized(NumVariables());
 
            unsigned affine_group_counter = 0;
            unsigned hom_group_counter = 0;
 
            unsigned dehom_index = 0; // index into x, the user-coordinates point
            unsigned hom_index = 0; // index into x_homogenized, the point we are computing
 
            for (auto& iter : time_order_of_variable_groups_)
            {
                switch (iter){
                    case VariableGroupType::Affine:
                    {
                        x_homogenized(hom_index++) = T(1);
                        for (unsigned ii = 0; ii < variable_groups_[affine_group_counter].size(); ++ii)
                            x_homogenized(hom_index++) = x(dehom_index++);
                        affine_group_counter++;
                        break;
                    }
                    case VariableGroupType::Homogeneous:
                    {
                        for (unsigned ii = 0; ii < hom_variable_groups_[hom_group_counter].size(); ++ii)
                            x_homogenized(hom_index++) = x(dehom_index++);
                        hom_group_counter++;
                        break;
                    }
                    case VariableGroupType::Ungrouped:
                    {
                        x_homogenized(hom_index++) = x(dehom_index++);
                        break;
                    }
                    default:
                    {
                        throw std::runtime_error("unacceptable VariableGroupType in HomogenizePointFIFO");
                    }
                }
            }
 
            return x_homogenized;
        }
 
        // --- the System's polynomial block (its functions live here after the fold) ---

        blocks::PolynomialBlock& PolyBlock()
        {
            for (auto& b : blocks_)
                if (auto* p = std::get_if<blocks::PolynomialBlock>(&b))
                    return *p;
            blocks_.emplace_back(blocks::PolynomialBlock{});
            return std::get<blocks::PolynomialBlock>(blocks_.back());
        }

        blocks::PolynomialBlock const* PolyBlockPtr() const
        {
            for (auto const& b : blocks_)
                if (auto* p = std::get_if<blocks::PolynomialBlock>(&b))
                    return p;
            return nullptr;
        }

        std::vector<Nd> const& PolyFunctions() const
        {
            static const std::vector<Nd> none;
            auto* p = PolyBlockPtr();
            return p ? p->Functions() : none;
        }

        void InvalidateDifferentiation() const
        {
            is_differentiated_ = false;
            if (auto* p = PolyBlockPtr())
                p->Invalidate();
        }

        void SyncPolyBlock() const
        {
            if (auto* p = PolyBlockPtr())
            {
                p->SetVariableOrdering(Variables());
                if (have_path_variable_) p->SetPathVariable(path_variable_);
                else                     p->ClearPathVariable();
            }
        }

        void ConstructOrdering() const;

        template <typename T>
        T CurrentPathValue() const
        {
            if (have_path_variable_)
                return std::get<T>(current_path_value_);
            return T(0);
        }

        template <typename Derived>
        void CoerceBlockOutputPrecision(Eigen::MatrixBase<Derived>& result) const
        {
            using Scalar = typename Derived::Scalar;
            if constexpr (!std::is_same<Scalar, complex_dbl>::value)
            {
                using bertini::Precision;
                // the working precision IS the staged point's -- derived, not stored (ADR-0057)
                auto const& staged = std::get<Vec<complex_mp> >(current_variable_values_);
                Precision(result, staged.size() ? Precision(staged) : DefaultPrecision());
            }
        }

        template <typename T>
        void EvalBlocksInPlace(Vec<T>& function_values) const
        {
            const auto& vars = std::get<Vec<T> >(current_variable_values_);
            const T t = CurrentPathValue<T>();
            Eigen::Index row = 0;
            for (auto const& blk : blocks_)
                std::visit([&](auto const& b){
                    const Eigen::Index n = static_cast<Eigen::Index>(b.NumFunctions());
                    b.template EvalInPlace<T>(function_values.segment(row, n), vars, t);
                    row += n;
                }, blk);
        }

        template <typename T>
        void JacobianBlocksInPlace(Mat<T>& J) const
        {
            const auto& vars = std::get<Vec<T> >(current_variable_values_);
            const T t = CurrentPathValue<T>();
            Eigen::Index row = 0;
            for (auto const& blk : blocks_)
                std::visit([&](auto const& b){
                    const Eigen::Index n = static_cast<Eigen::Index>(b.NumFunctions());
                    Mat<T> jb(n, J.cols());                 // blocks write into a contiguous target
                    b.template JacobianInPlace<T>(jb, vars, t);
                    J.block(row, 0, n, J.cols()) = jb;
                    row += n;
                }, blk);
        }

        template <typename T>
        void TimeDerivBlocksInPlace(Vec<T>& ds_dt) const
        {
            const auto& vars = std::get<Vec<T> >(current_variable_values_);
            const T t = CurrentPathValue<T>();
            Eigen::Index row = 0;
            for (auto const& blk : blocks_)
                std::visit([&](auto const& b){
                    const Eigen::Index n = static_cast<Eigen::Index>(b.NumFunctions());
                    b.template TimeDerivInPlace<T>(ds_dt.segment(row, n), vars, t);
                    row += n;
                }, blk);
        }
 
 
        // Append this system's canonical encoding to `out`, sharing `ctx` so operand systems
        // (randomization / blend) and cross-function node sharing back-reference deterministically
        // within one encoding unit.  The public CanonicalEncodingText() wraps this with a fresh
        // context; the recursion for operand systems passes the same one.  (ADR-0042)
        void EncodeCanonical(std::ostream& out, node::EncodingContext& ctx) const;
 
        // Read this system's canonical encoding from `cursor` (the inverse of EncodeCanonical,
        // section by section), sharing `ctx` so back-references into operand systems resolve
        // the way the encoder numbered them.  Replaces this system's content wholesale.
        void DecodeCanonicalFrom(node::DecodingCursor& cursor, node::DecodingContext& ctx);
 
        // Guard for structural mutators: throws std::logic_error naming `operation` if the
        // system is sealed (ADR-0042).  First line of every structural mutator.
        void ThrowIfSealed(char const* operation) const;
 
        bool is_sealed_ = false; 
        mutable std::optional<detail::Digest256> sealed_digest_; 
 
        VariableGroup ungrouped_variables_; 
        std::vector< VariableGroup > variable_groups_; 
        std::vector< VariableGroup > hom_variable_groups_; 
 
        VariableGroup homogenizing_variables_; 
 
        std::vector<unsigned> auxiliary_variable_groups_; 
        std::vector<unsigned> auxiliary_coordinates_; 
 
        std::vector<Nd> pre_homogenization_functions_; 
 
 
        bool have_path_variable_ = false; 
        Var path_variable_; 
 
        VariableGroup implicit_parameters_; 
        std::vector< NE > explicit_parameters_; 
 
        // The polynomial path -- functions_, subfunctions_, constant_subfunctions_, their
        // derivatives, the SLP, and eval_method_ -- has been folded into a
        // blocks::PolynomialBlock held in blocks_ (see PolyBlock()/PolyBlockPtr()).  The System
        // is now a thin orchestrator over blocks + variable groups + patch.
 
        class Patch patch_; 
        bool is_patched_ = false;   
 
        mutable bool is_differentiated_ = false; 
 
        std::vector<Block> blocks_; 
 
        std::vector< VariableGroupType > time_order_of_variable_groups_;
 
        mutable std::tuple< Vec<complex_dbl>, Vec<complex_mp> > current_variable_values_;
        mutable std::tuple< complex_dbl, complex_mp > current_path_value_{}; 
 
        mutable VariableGroup variable_ordering_; 
        mutable bool have_ordering_ = false;
 
 
 
 
 
 
 
 
        friend class boost::serialization::access;
 
        /*definition of serialize function*/
        template<class Archive>
        void serialize(Archive & ar, const unsigned int /*version*/){
            ar & ungrouped_variables_;
            ar & variable_groups_;
            ar & hom_variable_groups_;
 
            ar & homogenizing_variables_;
 
            ar & auxiliary_variable_groups_;
            ar & auxiliary_coordinates_;
 
            ar & have_path_variable_;
            ar & path_variable_;
 
            ar & implicit_parameters_;
            ar & explicit_parameters_;
 
            ar & patch_;
            ar & is_patched_;
 
            // The polynomial path (functions / subfunctions / derivatives / SLP / eval+deriv
            // methods) now lives inside the PolynomialBlock, which is archived as part of blocks_.
            ar & blocks_;
 
 
            // now for the cached / mutable things
            // precision_ is gone (ADR-0057): a System has no precision of its own, and
            // transient evaluation state is not archived (see current_variable_values_ below)
 
            ar & is_differentiated_;
 
            ar & time_order_of_variable_groups_;
 
            ar & pre_homogenization_functions_;
 
            // seal flag round-trips; the memoized digest is recomputed on demand after load
            // (guards against archives written under an older encoding version) (ADR-0042)
            ar & is_sealed_;
 
            // if (Archive::is_loading::value == true){
                // have_ordering_ = false;}
            // else
            // {
                ar & have_ordering_;
                ar & variable_ordering_;
            // }
 
 
            // if (Archive::is_loading::value == true){
            //  std::get<Vec<complex_dbl>>(current_variable_values_).resize(NumVariables());
            //  std::get<Vec<complex_mp>>(current_variable_values_).resize(NumVariables());
            // }
            // // next two lines no matter what
            // ar & std::get<Vec<complex_dbl>>(current_variable_values_);
            // ar & std::get<Vec<complex_mp>>(current_variable_values_);
 
            // A loaded system joins the live node universe (ADR-0068).  An archive builds nodes
            // through their default constructors, never through Make, so what comes out is
            // content-equal to the original and shares none of its nodes: its x is not THE x,
            // and it could not be combined with the system it was copied from.  Every load
            // goes through here -- pickle, copy, deepcopy, the MPI broadcast, the records
            // loader -- so the repair is made here rather than asked of each of them.
            if (Archive::is_loading::value)
            {
                node::ReinternMemo memo;
                ReinternNodes(memo);
            }
        }
 
 
    public:
        EIGEN_MAKE_ALIGNED_OPERATOR_NEW
    };
 

    System Concatenate(System sys1, System const& sys2);
 
 

    std::string UniquePathVariableName(System const& target, std::string base = "t");

    System MakeHomotopy(System const& target, System const& start,
                        std::string const& path_variable_name = "",
                        std::shared_ptr<node::Node> const& gamma = nullptr);

    void CheckVariableStructuresMatch(System const& a, System const& b,
                                      std::string const& operation,
                                      std::string const& a_name, std::string const& b_name);

    void AdoptVariableStructure(System& into, System const& from);

    System MakeMovingHomotopy(System const& fixed, System const& start_moving, System const& end_moving,
                              std::string const& path_variable_name = "",
                              std::shared_ptr<node::Node> const& gamma = nullptr);

    System Clone(System const& sys);

    void Simplify(System & sys);

    std::shared_ptr<const System> InternSystem(std::shared_ptr<System> const& candidate);

    std::shared_ptr<const System> LoadSystemUnified(std::istream& in);

    inline constexpr char SystemEncodingVersion[] = "b2sysenc/2";

    // Explicit instantiation declarations for the two concrete numeric types.
    // Definitions live in core/src/system/system.cpp.
    // Suppresses re-instantiation of the heavy Eval/Jacobian/Set template bodies
    // (with their eval_method_ switch trees) in every including TU.
 
    extern template void System::EvalInPlace<complex_dbl>(Vec<complex_dbl>&) const;
    extern template void System::EvalInPlace<complex_mp>(Vec<complex_mp>&) const;
 
    extern template Vec<complex_dbl> System::Eval<complex_dbl>() const;
    extern template Vec<complex_mp> System::Eval<complex_mp>() const;
 
    extern template void System::JacobianInPlace<complex_dbl>(Mat<complex_dbl>&) const;
    extern template void System::JacobianInPlace<complex_mp>(Mat<complex_mp>&) const;
 
    extern template Mat<complex_dbl> System::Jacobian<complex_dbl>() const;
    extern template Mat<complex_mp> System::Jacobian<complex_mp>() const;
 
    extern template Mat<complex_dbl> System::Jacobian<complex_dbl>(const Vec<complex_dbl>&) const;
    extern template Mat<complex_mp> System::Jacobian<complex_mp>(const Vec<complex_mp>&) const;
 
    extern template void System::JacobianInPlace<complex_dbl>(Mat<complex_dbl>&, const Vec<complex_dbl>&) const;
    extern template void System::JacobianInPlace<complex_mp>(Mat<complex_mp>&, const Vec<complex_mp>&) const;
 
    extern template void System::TimeDerivativeInPlace<complex_dbl>(Vec<complex_dbl>&) const;
    extern template void System::TimeDerivativeInPlace<complex_mp>(Vec<complex_mp>&) const;
 
    extern template Vec<complex_dbl> System::TimeDerivative<complex_dbl>() const;
    extern template Vec<complex_mp> System::TimeDerivative<complex_mp>() const;
 
    extern template void System::SetVariables<complex_dbl>(const Vec<complex_dbl>&) const;
    extern template void System::SetVariables<complex_mp>(const Vec<complex_mp>&) const;
 
    extern template void System::SetPathVariable<complex_dbl>(complex_dbl const&) const;
    extern template void System::SetPathVariable<complex_mp>(complex_mp const&) const;
 
    extern template void System::SetAndReset<complex_dbl>(Vec<complex_dbl> const&, complex_dbl const&) const;
    extern template void System::SetAndReset<complex_mp>(Vec<complex_mp> const&, complex_mp const&) const;
 
    extern template void System::SetAndReset<complex_dbl>(Vec<complex_dbl> const&) const;
    extern template void System::SetAndReset<complex_mp>(Vec<complex_mp> const&) const;
 
}
 
 
 
 
 
 
 
 
 
#endif // for the ifndef include guards