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) 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_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 <stdexcept>
#include <string>
#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<complex_mp> >& coefficients_mpfr = std::get<std::vector<Vec<complex_mp> > >(this->coefficients_working_);
            std::vector<Vec<complex_dbl> >& coefficients_dbl = std::get<std::vector<Vec<complex_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) = complex_dbl(coefficients_highest_precision_[ii](jj));
                    coefficients_mpfr[ii](jj) = complex_mp(coefficients_highest_precision_[ii](jj));
 
                    assert(coefficients_highest_precision_[ii](jj) == other.coefficients_highest_precision_[ii](jj));
                }
            }
        }

        Patch& operator=(Patch const& other) = default;
 

        Patch(std::vector<unsigned> const& sizes) : coefficients_highest_precision_(sizes.size()), variable_group_sizes_(sizes), precision_(DefaultPrecision())
        {
            using bertini::Precision;
            using bertini::multiprecision::RandomComplex;
 
            std::vector<Vec<complex_mp> >& coefficients_mpfr = std::get<std::vector<Vec<complex_mp> > >(coefficients_working_);
            std::vector<Vec<complex_dbl> >& coefficients_dbl = std::get<std::vector<Vec<complex_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)
                {
                    // bounded-modulus draw (away from 0 and infinity), matching the linear-product /
                    // mhom / binomial start-system coefficients -- a heavy-tailed coefficient here scales
                    // the patch equation badly and feeds the same near-t=0 conditioning trouble.
                    multiprecision::RandomComplexBoundedModulusAssign(coefficients_highest_precision_[ii](jj), MaxPrecisionAllowed());
                }
 
                coefficients_dbl[ii].resize(sizes[ii]);
                for (unsigned jj=0; jj<sizes[ii]; ++jj)
                {
                    coefficients_dbl[ii](jj) = complex_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_);
            }
        }
 

        Patch(std::vector<unsigned> const& sizes, std::vector<Vec<complex_mp>> const& coefficients)
            : variable_group_sizes_(sizes), precision_(DefaultPrecision())
        {
            if (coefficients.size() != sizes.size())
                throw std::invalid_argument("Patch: " + std::to_string(coefficients.size())
                    + " coefficient vectors for " + std::to_string(sizes.size()) + " variable groups");
 
            std::vector<Vec<complex_mp> >& coefficients_mpfr = std::get<std::vector<Vec<complex_mp> > >(coefficients_working_);
            std::vector<Vec<complex_dbl> >& coefficients_dbl = std::get<std::vector<Vec<complex_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)
            {
                if (static_cast<size_t>(coefficients[ii].size()) != sizes[ii])
                    throw std::invalid_argument("Patch: group " + std::to_string(ii) + " has "
                        + std::to_string(coefficients[ii].size()) + " coefficients for "
                        + std::to_string(sizes[ii]) + " variables");
 
                coefficients_highest_precision_[ii].resize(sizes[ii]);
                coefficients_dbl[ii].resize(sizes[ii]);
                coefficients_mpfr[ii].resize(sizes[ii]);
                for (unsigned jj=0; jj<sizes[ii]; ++jj)
                {
                    // set the precision first so the copy is exact at the source's precision
                    coefficients_highest_precision_[ii](jj).precision(coefficients[ii](jj).precision());
                    coefficients_highest_precision_[ii](jj) = coefficients[ii](jj);
                    coefficients_dbl[ii](jj) = complex_dbl(coefficients_highest_precision_[ii](jj));
                    coefficients_mpfr[ii](jj) = coefficients_highest_precision_[ii](jj);
                    coefficients_mpfr[ii](jj).precision(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<complex_mp> >& coefficients_mpfr = std::get<std::vector<Vec<complex_mp> > >(p.coefficients_working_);
            std::vector<Vec<complex_dbl> >& coefficients_dbl = std::get<std::vector<Vec<complex_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)
                    // real bounded-modulus draw: away from 0 and infinity (same recipe as the complex
                    // patch / start systems) but kept REAL, so a real patch keeps a real path real.
                    multiprecision::RandomRealBoundedModulusAssign(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) = complex_dbl(p.coefficients_highest_precision_[ii](jj));
            }
 
            return p;
        }
 

        unsigned Precision() const
        {
            return precision_;
        }

        template <typename T>
        void SyncPrecision(Vec<T> const& x) const
        {
            if constexpr (!std::is_same<T, complex_dbl>::value)
            {
                if (x.size() > 0)
                {
                    const unsigned p = bertini::Precision(x(0));
                    if (p != precision_)
                        Precision(p);
                }
            }
        }

        void Precision(unsigned new_precision) const
        {
            if (precision_==new_precision)
                return;
 
            using bertini::Precision;
            std::vector<Vec<complex_mp> >& coefficients_mpfr = std::get<std::vector<Vec<complex_mp> > >(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");
 
            SyncPrecision(x);
 
            #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());
            }
            // Formerly an assert that the input's precision MATCHED the patch's -- disabled long
            // ago because it could not be honoured.  A patch is an evaluable type, so it now
            // self-aligns instead of demanding to be aligned for: ADR-0057.
//                );
            #endif
 
            // unpack from the tuple of working coefficients
            const std::vector<Vec<T> >& coefficients = std::get<std::vector<Vec<T> > >(coefficients_working_);
 
            unsigned offset(static_cast<unsigned>(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");
 
            SyncPrecision(x);   // the jacobian is constant in x, but it is written FROM the
                                // working coefficients, so they must be at x's precision
 
 
            #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(static_cast<unsigned>(jacobian.rows() - NumVariableGroups())); // by precondition this number is at least 0.  the precondition is ensured by the public wrapper
 
            // A patch row is sparse -- one block of coefficients per variable group, zero elsewhere.
            // Zero the rows the patch owns before writing those coefficients, so the patch FULLY
            // defines its own rows and the caller need not pre-zero the matrix.  The block-composed
            // Jacobian path allocates J uninitialized and assigns only the block (function) rows;
            // without this the patch rows' off-coefficient entries are read uninitialized -- benign
            // (zeroed pages) on Linux/macOS, but garbage on Windows, where a degree-2 homotopy's
            // Jacobian "evaluated" to ~1e252 and wrecked the AMP condition-number estimate.
            jacobian.block(offset, 0, NumVariableGroups(), jacobian.cols()).setZero();
 
            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
        {
            // Self-align to the point, exactly as EvalInPlace and the Jacobian do -- ADR-0057.
            // Without this, rescaling read whatever precision the working coefficients happened to
            // be left at by the last evaluation, and the rescaled point came back REPORTING the
            // requested precision while carrying only the older one's correct digits.  Nothing
            // downstream could see it: the assert that used to stand here demanded the caller align
            // the precisions beforehand, and asserts are compiled out of a release build, so the
            // one check that could have caught it was absent exactly where it mattered (b2#461).
            SyncPrecision(x);
 
            #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.");
            #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 static_cast<unsigned>(variable_group_sizes_.size());
        }
 

        unsigned NumVariables() const
        {
            unsigned num_vars(0);
            for (auto v : variable_group_sizes_)
                num_vars += v;
            return num_vars;
        }

        std::vector< Vec< complex_mp > > const& Coefficients() const
        {
            return coefficients_highest_precision_;
        }

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

        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< complex_mp > > coefficients_highest_precision_; 
 
        mutable std::tuple< std::vector< Vec< complex_mp > >, std::vector< Vec< complex_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