patch.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
//include/bertini2/system/patch.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.
//
//include/bertini2/system/patch.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 include/bertini2/system/patch.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_PATCH_HPP
#define BERTINI_PATCH_HPP
 
#include "bertini2/mpfr_complex.hpp"
 
#include "bertini2/mpfr_extensions.hpp"
#include "bertini2/num_traits.hpp"
#include "bertini2/eigen_extensions.hpp"
 
#include <vector>
 
#include <Eigen/Dense>
 
#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>
 
#include <boost/serialization/vector.hpp>
 
namespace bertini {
 
 
 

    class Patch
    {
 
    public:

        Patch() : precision_(DefaultPrecision())
        {}
 

        Patch(Patch const& other)
        {   
            variable_group_sizes_ = other.variable_group_sizes_;
            precision_ = DefaultPrecision();
 
            // a little shorthand unpacking the tuple
            std::vector<Vec<mpfr_complex> >& coefficients_mpfr = std::get<std::vector<Vec<mpfr_complex> > >(this->coefficients_working_);
            std::vector<Vec<dbl> >& coefficients_dbl = std::get<std::vector<Vec<dbl> > >(this->coefficients_working_);
 
            coefficients_highest_precision_.resize(other.NumVariableGroups());
            coefficients_mpfr.resize(variable_group_sizes_.size());
            coefficients_dbl.resize(variable_group_sizes_.size());
 
            for (unsigned ii(0); ii<other.NumVariableGroups(); ++ii)
            {
                auto curr_size = variable_group_sizes_[ii];
 
                coefficients_highest_precision_[ii].resize(curr_size);
                coefficients_dbl[ii].resize(curr_size);
                coefficients_mpfr[ii].resize(curr_size);
 
                for (unsigned jj(0); jj<curr_size; jj++)
                {
                    coefficients_highest_precision_[ii](jj).precision(other.coefficients_highest_precision_[ii](jj).precision());
 
                    coefficients_highest_precision_[ii](jj) = other.coefficients_highest_precision_[ii](jj);
 
                    coefficients_dbl[ii](jj) = dbl(coefficients_highest_precision_[ii](jj));
                    coefficients_mpfr[ii](jj) = mpfr_complex(coefficients_highest_precision_[ii](jj));
 
                    assert(coefficients_highest_precision_[ii](jj) == other.coefficients_highest_precision_[ii](jj));
                }
            }
        }
 

        Patch(std::vector<unsigned> const& sizes) : variable_group_sizes_(sizes), coefficients_highest_precision_(sizes.size()), precision_(DefaultPrecision())
        {
            using bertini::Precision;
            using bertini::multiprecision::RandomComplex;
 
            std::vector<Vec<mpfr_complex> >& coefficients_mpfr = std::get<std::vector<Vec<mpfr_complex> > >(coefficients_working_);
            std::vector<Vec<dbl> >& coefficients_dbl = std::get<std::vector<Vec<dbl> > >(coefficients_working_);
 
            coefficients_highest_precision_.resize(sizes.size());
            coefficients_dbl.resize(sizes.size());
            coefficients_mpfr.resize(sizes.size());
 
            for (size_t ii=0; ii<sizes.size(); ++ii)
            {
                // this produces coefficients at maximum precision.
                coefficients_highest_precision_[ii].resize(sizes[ii]);
                for (unsigned jj=0; jj<sizes[ii]; ++jj)
                {
                    multiprecision::RandomComplexAssign(coefficients_highest_precision_[ii](jj), MaxPrecisionAllowed());
                }
 
                coefficients_dbl[ii].resize(sizes[ii]);
                for (unsigned jj=0; jj<sizes[ii]; ++jj)
                {
                    coefficients_dbl[ii](jj) = dbl(coefficients_highest_precision_[ii](jj));
                }
 
                // assignment preserves precision of source.  
                // https://github.com/boostorg/multiprecision/issues/75
                coefficients_mpfr[ii].resize(sizes[ii]);
                for (unsigned jj=0; jj<sizes[ii]; ++jj)
                {
                    coefficients_mpfr[ii](jj) = coefficients_highest_precision_[ii](jj);
                    coefficients_mpfr[ii](jj).precision(precision_);
                }
 
                assert(Precision(coefficients_mpfr[ii](0))==precision_);
            }
        }
 

        static Patch Random(std::vector<unsigned> const& sizes)
        {
            return Patch(sizes);
        }
 

        static Patch RandomReal(std::vector<unsigned> const& sizes)
        {   
            using bertini::multiprecision::RandomReal;
            using bertini::Precision;
 
            Patch p;
 
            p.variable_group_sizes_ = sizes;
 
            
 
            std::vector<Vec<mpfr_complex> >& coefficients_mpfr = std::get<std::vector<Vec<mpfr_complex> > >(p.coefficients_working_);
            std::vector<Vec<dbl> >& coefficients_dbl = std::get<std::vector<Vec<dbl> > >(p.coefficients_working_);
 
            p.coefficients_highest_precision_.resize(sizes.size());
            coefficients_mpfr.resize(sizes.size());
            coefficients_dbl.resize(sizes.size());
            
            for (size_t ii=0; ii<sizes.size(); ++ii)
            {
                p.coefficients_highest_precision_[ii].resize(sizes[ii]);
                for (unsigned jj=0; jj<sizes[ii]; ++jj)
                    multiprecision::RandomRealAssign(p.coefficients_highest_precision_[ii](jj), MaxPrecisionAllowed());
 
                coefficients_mpfr[ii] = p.coefficients_highest_precision_[ii]; 
                Precision(coefficients_mpfr[ii],DefaultPrecision());
                assert(Precision(coefficients_mpfr[ii](0))==DefaultPrecision());
 
                coefficients_dbl[ii].resize(sizes[ii]);
                for (unsigned jj=0; jj<sizes[ii]; ++jj)
                    coefficients_dbl[ii](jj) = dbl(p.coefficients_highest_precision_[ii](jj));
            }
 
            return p;
        }
 

        unsigned Precision() const
        {
            return precision_;
        }

        void Precision(unsigned new_precision) const
        {
            if (precision_==new_precision)
                return;
 
            using bertini::Precision;
            std::vector<Vec<mpfr_complex> >& coefficients_mpfr = std::get<std::vector<Vec<mpfr_complex> > >(coefficients_working_);
 
            for (unsigned ii = 0; ii < NumVariableGroups(); ++ii)
            {
                for (unsigned jj=0; jj<variable_group_sizes_[ii]; ++jj)
                {
                    if (new_precision>precision_)
                    {
                        coefficients_mpfr[ii](jj) = coefficients_highest_precision_[ii](jj);
                        coefficients_mpfr[ii](jj).precision(new_precision);
                    }
                    else
                        coefficients_mpfr[ii](jj).precision(new_precision);
 
                    assert(Precision(coefficients_mpfr[ii](jj))==new_precision);
                }
            }
            
            precision_ = new_precision;
        }
 
 
        
 
        

        template<typename Derived, typename T>
        void EvalInPlace(Eigen::MatrixBase<Derived> & function_values, Vec<T> const& x) const
        {
            static_assert(std::is_same<typename Derived::Scalar,T>::value,"scalar types must match");
 
            #ifndef BERTINI_DISABLE_ASSERTS
            if (! (function_values.size()>=NumVariableGroups()) )
            {
                std::stringstream ss;
                ss << "container for function values must be of length at least as long as the number of variable groups.  the input vector into which to write is of length " << function_values.size();
                throw std::runtime_error(ss.str());
            }
//          assert((bertini::Precision(x(0))==DoublePrecision() || bertini::Precision(x(0)) == Precision())
//                  && "precision of input vector must match current working precision of patch during evaluation"
//                );
            #endif
 
            // unpack from the tuple of working coefficients
            const std::vector<Vec<T> >& coefficients = std::get<std::vector<Vec<T> > >(coefficients_working_);
 
            unsigned offset(function_values.size() - NumVariableGroups()); // by precondition this number is at least 0.  the precondition is ensured by the public wrapper
            unsigned counter(0);
            for (unsigned ii = 0; ii < NumVariableGroups(); ++ii)
            {
                T& value = function_values(ii+offset);
                value = T(-1);
                for (unsigned jj=0; jj<variable_group_sizes_[ii]; ++jj)
                {   
                    value += x(counter)*coefficients[ii](jj);
                    counter++;
                }
            }
        }

        template<typename T>
        Vec<T> Eval(Vec<T> const& x) const
        {
            Vec<T> function_values = Vec<T>::Zero(NumVariableGroups());
            EvalInPlace(function_values, x);
            return function_values;
        }

        template<typename Derived, typename T>
        void JacobianInPlace(Eigen::MatrixBase<Derived> & jacobian, Vec<T> const& x) const
        {
            static_assert(std::is_same<typename Derived::Scalar,T>::value,"scalar types must match");
 
 
            #ifndef BERTINI_DISABLE_ASSERTS
            assert(jacobian.rows()>=NumVariableGroups() && "input jacobian must have at least as many rows as variable groups");
            assert(jacobian.cols()==NumVariables() && "input jacobian must have as many columns as the patch has variables");
            assert(
                   (bertini::Precision(x(0))==DoublePrecision() || bertini::Precision(x(0)) == Precision())  
                        && "precision of input vector must match current working precision of patch during evaluation"
                   );
            #endif
            
            const std::vector<Vec<T> >& coefficients = std::get<std::vector<Vec<T> > >(coefficients_working_);
 
            unsigned offset(jacobian.rows() - NumVariableGroups()); // by precondition this number is at least 0.  the precondition is ensured by the public wrapper
            unsigned counter(0);
            for (unsigned ii = 0; ii < NumVariableGroups(); ++ii)
                for (unsigned jj=0; jj<variable_group_sizes_[ii]; ++jj)
                    jacobian(ii+offset,counter++) = coefficients[ii](jj);
        }

        template<typename T>
        Mat<T> Jacobian(Vec<T> const& x) const
        {
            Mat<T> jacobian = Mat<T>::Zero(NumVariableGroups(), NumVariables());
            JacobianInPlace(jacobian, x);
            return jacobian;
        }
 
        
        template<typename T>
        void RescalePointToFitInPlace(Vec<T> & x) const
        {
            #ifndef BERTINI_DISABLE_ASSERTS
                assert(x.size() == NumVariables() && "input point for rescaling to fit a patch must have same length as total number of variables being patched, in all variable groups.");
                assert((bertini::Precision(x(0))==DoublePrecision() || bertini::Precision(x(0)) == Precision())
                        && "precision of input vector must match current working precision of patch during rescaling"
                       );
            #endif
 
            const std::vector<Vec<T> >& coefficients = std::get<std::vector<Vec<T> > >(coefficients_working_);
 
            unsigned starting_index_counter(0);
            for (unsigned ii=0; ii<NumVariableGroups(); ii++)
            {
                auto subvec = x.segment(starting_index_counter,variable_group_sizes_[ii]);
                subvec /= subvec.transpose() * coefficients[ii];
                starting_index_counter += variable_group_sizes_[ii];
            }
        }

        template <typename T>
        Vec<T> RescalePoint(Vec<T> const& x) const
        {
            Vec<T> x_rescaled = x;
            RescalePointToFitInPlace(x_rescaled);
            return x_rescaled;
        }

        unsigned NumVariableGroups() const
        {
            return variable_group_sizes_.size();
        }
 

        unsigned NumVariables() const
        {
            unsigned num_vars(0);
            for (auto v : variable_group_sizes_)
                num_vars += v;
            return num_vars;
        }
 
 
        friend std::ostream& operator<<(std::ostream & out, Patch const& p)
        {
            out << p.NumVariableGroups() << " variable groups being patched\n";
            unsigned counter(0);
            for (auto& c : p.coefficients_highest_precision_)
            {
                out << "patch " << counter++ << " has " << c.size() << " coefficients:\n";
                out << c << "\n";
            }
            out << "current patch precision: " << p.Precision() << "\n";
            return out;
        }

        bool operator==(Patch const& rhs) const
        {
            if (NumVariableGroups()!=rhs.NumVariableGroups())
                return false;
            if (NumVariables()!=rhs.NumVariables())
                return false;
            for (unsigned ii=0; ii<NumVariableGroups(); ii++)
                if (variable_group_sizes_[ii]!=rhs.variable_group_sizes_[ii])
                    return false;
            for (unsigned ii=0; ii<NumVariableGroups(); ii++)
                if (coefficients_highest_precision_[ii]!=rhs.coefficients_highest_precision_[ii])
                    return false;
 
 
            return true;
        }

        bool operator!=(Patch const& rhs) const
        {
            return !(*this==rhs);
        }
 
    private:

        //
        //    Data members
        //
 
        std::vector< Vec< mpfr_complex > > coefficients_highest_precision_; 
 
        mutable std::tuple< std::vector< Vec< mpfr_complex > >, std::vector< Vec< dbl > > > coefficients_working_; 
 
        std::vector<unsigned> variable_group_sizes_; 
 
        mutable unsigned precision_; 
 
        // add serialization support through boost.
 
        friend class boost::serialization::access;
 
        template <typename Archive>
        void serialize(Archive& ar, const unsigned version) {
            ar & precision_;
 
            ar & coefficients_highest_precision_;
 
            ar & std::get<0>(coefficients_working_);
            ar & std::get<1>(coefficients_working_);
            ar & variable_group_sizes_;
            
        }
 
    };
 
 
}
 
 
#endif  // re: include guards