system/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) 2015 - 2021 by 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.
 
// individual authors of this file include:
// silviana amethyst, university of wisconsin eau claire

 
#ifndef BERTINI_SYSTEM_HPP
#define BERTINI_SYSTEM_HPP
 
#include <assert.h>
#include <vector>
 
 
#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/type_index.hpp>
 
#include "bertini2/mpfr_complex.hpp"
#include "bertini2/mpfr_extensions.hpp"
#include "bertini2/eigen_extensions.hpp"
 
 
#include "bertini2/function_tree.hpp"
#include "bertini2/system/patch.hpp"
 
#include "bertini2/system/straight_line_program.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 {
 
    
    enum class EvalMethod
    {
        FunctionTree, // using virtual methods and recursion
        SLP // using straight line programs
            // now!  20230714, Eindhoven, Netherlands
    };
 
    enum class DerivMethod
    {
        JacobianNode, // using Jacobian nodes, which are either 1 or 0 when evaluated based on the variable of differentiation
        Derivatives // classic differentiation, using more space in memory but not requiring a variable of differentation when evaluatiing
    };
 

    EvalMethod DefaultEvalMethod();
 
    DerivMethod DefaultDerivMethod();
 

    bool DefaultAutoSimplify();

    class System{
    
    public:
        // a few local using statements to reduce typing etc.
        using Fn = std::shared_ptr<node::Function>;
        using Var = std::shared_ptr<node::Variable>;
        using Nd = std::shared_ptr<node::Node>;
        using Jac = std::shared_ptr<node::Jacobian>;
        
        System() : is_differentiated_(false), have_path_variable_(false), have_ordering_(false), precision_(DefaultPrecision()), is_patched_(false)
        {}

        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);

        void precision(unsigned new_precision) const;

        unsigned precision() const
        {
            return precision_;
        }

        void Differentiate() const;
 
 
        
        void ResetFunctions() const
        {
            // TODO: it has the unfortunate side effect of resetting constant functions, too.
            switch (eval_method_){
            case EvalMethod::FunctionTree:
                for (const auto& iter : functions_) 
                    iter->Reset();
                break;
            case EvalMethod::SLP:
                // nothing
                break;
            }   
            
        }

        void ResetJacobian() const
        {
            switch (eval_method_)
            {
                case EvalMethod::FunctionTree:{
 
                    switch (deriv_method_){
                        case DerivMethod::JacobianNode:
                        {
                            for (const auto& iter : jacobian_) 
                                iter->Reset();
                            break;
                        }
                        case DerivMethod::Derivatives:
                        {
                            for (const auto& iter : space_derivatives_) 
                                iter->Reset();
                            break;
                        }
                    }
 
                    break;
                }
                case EvalMethod::SLP:
                {
                    // nothing to do, it's not a resetting kind of thing.
                    break;                  
                }
 
            }
        }
 
        void ResetTimeDerivatives() const
        {
            switch (eval_method_)
            {
                case EvalMethod::FunctionTree:{
                    
                    switch (deriv_method_){
                        case DerivMethod::JacobianNode:
                        {
                            for (const auto& iter : jacobian_) 
                                iter->Reset();
                            break;
                        }
                        case DerivMethod::Derivatives:
                        {
                            for (const auto& iter : time_derivatives_) 
                                iter->Reset();
                            break;
                        }
                    }
 
                    break;
                }
                case EvalMethod::SLP:
                {
                    // nothing to do, it's not a resetting kind of thing.
                    break;                  
                }
 
            }
        }

        void Reset() const
        {
            ResetFunctions();
            ResetJacobian();
            ResetTimeDerivatives();
        }
 
        template<typename T>
        void EvalInPlace(Vec<T> & function_values) const
        {
            
            if (function_values.size() != 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());
            }
 
            switch (eval_method_){
                case EvalMethod::FunctionTree:
                {
                    unsigned counter(0);
                    for (auto iter=functions_.begin(); iter!=functions_.end(); iter++, counter++) {
                        (*iter)->EvalInPlace<T>(function_values(counter));
                    }
                }
 
                case EvalMethod::SLP:
                    {
                        slp_.GetFuncValsInPlace<T>(function_values);
                    }
            }
 
 
            if (IsPatched())
                patch_.EvalInPlace(function_values,
                                    std::get<Vec<T> >(current_variable_values_)); // does a patch not have a caching mechanism?
                                    // .segment(NumNaturalFunctions(),NumTotalVariableGroups())
            
        }
        
        
        
        
        
        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()!=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());
            ResetFunctions();
            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()!=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());//TODO: remove this eval
            SetPathVariable(path_variable_value);
 
            ResetFunctions(); // todo, elimiante this.  i feel like setting the variables or path variable should be enough to set the flag/ take the action
 
            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() != NumTotalFunctions() || J.cols() != 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");
            }
            
            const auto& vars = Variables();
 
            if (!is_differentiated_)
                Differentiate();
 
            switch (eval_method_)
            {
                case EvalMethod::FunctionTree:
                {
                    switch (deriv_method_){
                        case DerivMethod::JacobianNode:{
                            for (int ii = 0; ii < NumNaturalFunctions(); ++ii)
                                for (int jj = 0; jj < NumVariables(); ++jj)
                                    jacobian_[ii]->EvalJInPlace<T>(J(ii,jj),vars[jj]);
                            break;
                        }
                        case DerivMethod::Derivatives:
                        {
                            for (int jj = 0; jj < NumVariables(); ++jj)
                                for (int ii = 0; ii < NumNaturalFunctions(); ++ii)
                                    space_derivatives_[ii+jj*NumNaturalFunctions()]->EvalInPlace<T>(J(ii,jj));
                            break;
                        }
                    }
                    break;
                } // function tree branch
 
                case EvalMethod::SLP:
                {
                    this->slp_.GetJacobianInPlace<T>(J); // the variable values should have been copied into place elsewhere.  that's not this function's responsibility.
                    break;                  
                }
            }
            
            if (IsPatched())
                patch_.JacobianInPlace(J,std::get<Vec<T> >(current_variable_values_));
            
        }
 
        
        
        
        
        
 

        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()!=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()!=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()!=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()); // TODO: remove this eval
            SetPathVariable(path_variable_value);
            ResetJacobian();
            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()!=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()); //TODO: remove this eval()
            SetPathVariable(path_variable_value);
            ResetTimeDerivatives();
            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()); //TODO: remove this eval()
            ResetTimeDerivatives();
            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() < 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();
 
            switch (eval_method_)
            {
                case EvalMethod::FunctionTree:{
                    switch (deriv_method_){
                        case DerivMethod::JacobianNode:
                        {
                            for (int ii = 0; ii < NumNaturalFunctions(); ++ii)
                                jacobian_[ii]->EvalJInPlace<T>(ds_dt(ii), path_variable_);
                            break;
                        }
                        case DerivMethod::Derivatives:
                        {
                            for (int ii = 0; ii < NumNaturalFunctions(); ++ii)
                                time_derivatives_[ii]->EvalInPlace<T>(ds_dt(ii));
                            break;
                        }
                    }
                break;
                } // function tree branch
 
                case EvalMethod::SLP:
                {
                    this->slp_.GetTimeDerivInPlace(ds_dt); // the variable values should have been copied into place elsewhere.  that's not this function's responsibility.
                    break;                  
                }
            }
 
            // the patch doesn't move with time.  derivatives 0.
            if (IsPatched())
                for (int ii = 0; ii < NumTotalVariableGroups(); ++ii)
                    ds_dt(ii+NumNaturalFunctions()) = T(0);
            
        }
 
        
        
        
        
        
        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();

        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;

        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()!= NumVariables())
                throw std::runtime_error("variable vector of different length from system-owned variables in SetVariables");
 
            const auto& vars = Variables();
 
            #ifndef BERTINI_DISABLE_PRECISION_CHECKS
                if (!std::is_same<T,dbl>::value && (Precision(new_values) != this->precision()))
                    throw std::runtime_error("precision of input point in SetVariables (" + std::to_string(Precision(new_values)) + ") must match the precision of the system (" + std::to_string(this->precision()) + ").");
 
                if (!std::is_same<T,dbl>::value && (vars[0]->node::NamedSymbol::precision() != this->precision()) )
                    throw std::runtime_error("internally, precision of variables (" + std::to_string(vars[0]->node::NamedSymbol::precision()) + ") in SetVariables must match the precision of the system (" + std::to_string(this->precision()) + ").");
            #endif
 
            if (!is_differentiated_)
                Differentiate();
 
 
            switch (eval_method_){
                case EvalMethod::FunctionTree:{
                    auto counter = 0;
 
                    for (auto iter=vars.begin(); iter!=vars.end(); iter++, counter++) {
                        (*iter)->set_current_value(new_values(counter));
                    }
 
                    std::get<Vec<T> >(current_variable_values_) = new_values;
                    break;
                }
                case EvalMethod::SLP:{
                    std::get<Vec<T> >(current_variable_values_) = new_values; // if this isn't here, then patch evaluation breaks.
                    slp_.SetVariableValues(new_values);
                    break;
                }
            } // switch
 
            
        }
 
 

        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");
 
            if (!is_differentiated_)
                Differentiate();
 
            switch (eval_method_){
                case EvalMethod::FunctionTree:{
                    path_variable_->set_current_value(new_value);
                    break;
                }
                case EvalMethod::SLP:{
                    path_variable_->set_current_value(new_value);
                    slp_.SetPathVariable(new_value);
                }
            }
        }
 
 
        template<typename T>
        void SetAndReset(Vec<T> const& new_space, T const& new_time) const
        {
            SetVariables(new_space);
            SetPathVariable(new_time);
 
            Reset();
        }
 
        template<typename T>
        void SetAndReset(Vec<T> const& new_space) const
        {
            SetVariables(new_space);
 
            Reset();
        }

        template<typename T>
        void SetImplicitParameters(Vec<T> new_values) const
        {
            if (new_values.size()!= implicit_parameters_.size())
                throw std::runtime_error("trying to set implicit parameter values, but there is a size mismatch");
 
            size_t counter = 0;
            for (auto iter=implicit_parameters_.begin(); iter!=implicit_parameters_.end(); iter++, counter++)
                (*iter)->set_current_value(new_values(counter));
 
        }
 
 
 
 
 
 
 

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

        void AddVariableGroup(VariableGroup const& v);
 

        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(Fn const& F);

        void AddParameters(std::vector<Fn> const& F);
 
 

        void AddSubfunction(Fn const& F);

        void AddSubfunctions(std::vector<Fn> const& F);
 
 

        void AddFunction(Fn const& F);

        void AddFunction(Nd const& F, std::string const& name = "unnamed_function");
 

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

        void AddConstant(Fn const& C);
 

        void AddConstants(std::vector<Fn> const& C);
 
 
 

        void AddPathVariable(Var const& v);
 

        bool HavePathVariable() const;
 
        auto& GetPathVariable() const{
            if (this->HavePathVariable())
                return this->path_variable_;
            else
                throw std::runtime_error("trying to get path variable for a system which doesn't have a path variable defined");
        }

        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()!=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);
            }
 
 
 

        auto Function(unsigned index) const
        {
            return functions_[index];
        }
 
        
        auto GetNaturalFunctions() const
        {
            return functions_;
        }
 
 

        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);
        }
        friend std::ostream& operator <<(std::ostream& out, const System & s);
 

        void ClearVariables();
 
 
        void CopyVariableStructure(System const& other);
        
        void DontAssumeUniformPrecision()
        {
            AssumeUniformPrecision(false);
        }
 
        void PleaseAssumeUniformPrecision()
        {
            AssumeUniformPrecision(true);
        }
 
        void AssumeUniformPrecision(bool val)
        {
            assume_uniform_precision_ = false;
        }

        auto IsAssumingUniformPrecision() const
        {
            return assume_uniform_precision_;
        }
 
        inline
        void PleaseAutoSimplify()
        {
            SetAutoSimplify(true);
        }
 
        inline
        void DontAutoSimplify()
        {
            SetAutoSimplify(false);
        }
 
        void SetAutoSimplify(bool val)
        {
            auto_simplify_ = val;
        }
 

        auto IsAutoSimplifying() const
        {
            return auto_simplify_;
        }

        void SimplifyFunctions();

        void SimplifyDerivatives() const;

        void Simplify();
 

        void SetEvalMethod(EvalMethod method)
        {
            eval_method_ = method;
        }

        EvalMethod  GetEvalMethod() const
        {
            return eval_method_;
        }
 
 
 

        void SetDerivMethod(DerivMethod method)
        {
            deriv_method_ = method;
        }

        DerivMethod  GetDerivMethod() const
        {
            return deriv_method_;
        }
 

        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:

        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 ungrouped_variable_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++);
                        break;
                    }
                    case VariableGroupType::Ungrouped:
                    {
                        x_dehomogenized(dehom_index++) = x(hom_index++);
                        ungrouped_variable_counter++;
                        break;
                    }
                    default:
                    {
                        throw std::runtime_error("unacceptable VariableGroupType in FIFOVariableOrdering");
                    }
                }
            }
 
            return x_dehomogenized;
        }
 
        void DifferentiateUsingDerivatives() const;
        void DifferentiateUsingJacobianNode() const;

        void ConstructOrdering() const;
 
 
        VariableGroup ungrouped_variables_; 
        std::vector< VariableGroup > variable_groups_; 
        std::vector< VariableGroup > hom_variable_groups_; 
 
        VariableGroup homogenizing_variables_; 
 
 
        bool have_path_variable_ = false; 
        Var path_variable_; 
        
        VariableGroup implicit_parameters_; 
        std::vector< Fn > explicit_parameters_; 
 
        std::vector< Fn > constant_subfunctions_; 
        std::vector< Fn > subfunctions_; 
        std::vector< Fn > functions_; 
        
        class Patch patch_; 
        bool is_patched_ = false;   
 
        mutable std::vector< Jac > jacobian_; 
 
        mutable std::vector< Nd > space_derivatives_; 
 
        mutable std::vector< Nd > time_derivatives_; 
 
        mutable bool is_differentiated_ = false; 
 
        mutable StraightLineProgram slp_; 
 
        std::vector< VariableGroupType > time_order_of_variable_groups_;
 
        mutable std::tuple< Vec<dbl>, Vec<mpfr_complex> > current_variable_values_;
 
        mutable VariableGroup variable_ordering_; 
        mutable bool have_ordering_ = false;
 
        mutable unsigned precision_; 
 
        bool assume_uniform_precision_ = false; 
 
        EvalMethod eval_method_ = DefaultEvalMethod(); 
        DerivMethod deriv_method_ = DefaultDerivMethod(); 
 
        bool auto_simplify_ = DefaultAutoSimplify();
 
 
 
 
 
 
        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 & have_path_variable_;
            ar & path_variable_;            
 
            ar & implicit_parameters_;
            ar & explicit_parameters_;
 
            ar & constant_subfunctions_;
            ar & subfunctions_;
            ar & functions_;
 
 
            ar & patch_;
            ar & is_patched_;
 
 
            ar & assume_uniform_precision_;
 
            ar & eval_method_;
            ar & deriv_method_;
 
            ar & auto_simplify_;
 
            // now for the cached / mutable things
            ar & precision_;
 
 
            // if (Archive::is_loading::value == true){
            //  is_differentiated_ = false;}
            // // else
            // // {
                ar & is_differentiated_;
                ar & jacobian_;
                ar & space_derivatives_;
                ar & time_derivatives_;
            // }
 
 
            ar & slp_; // does this need to be re-constructed after de-serialization?
 
            ar & time_order_of_variable_groups_;
 
            // 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<dbl>>(current_variable_values_).resize(NumVariables());
            //  std::get<Vec<mpfr_complex>>(current_variable_values_).resize(NumVariables());
            // }
            // // next two lines no matter what
            // ar & std::get<Vec<dbl>>(current_variable_values_);
            // ar & std::get<Vec<mpfr_complex>>(current_variable_values_);
 
        }
 
 
        EIGEN_MAKE_ALIGNED_OPERATOR_NEW
    };
 

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

    System Clone(System const& sys);

    void Simplify(System & sys);
}
 
 
 
 
 
 
 
 
 
#endif // for the ifndef include guards