amp_tracker.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
//amp_tracker.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.
//
//amp_tracker.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 amp_tracker.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_AMP_TRACKER_HPP
#define BERTINI_AMP_TRACKER_HPP
 
#include <assert.h>
 
#pragma once 
 
#include "bertini2/trackers/base_tracker.hpp"
 
 
namespace bertini{
 
    namespace tracking{
 
        
        using std::max;
        using std::min;
        using std::pow;
 
        using bertini::max;

        inline
        mpfr_float MinTimeForCurrentPrecision(unsigned precision, mpfr_float const& time_to_go, int safety_digits = 3)
        {
            mpfr_float t = pow( mpfr_float(10), safety_digits-long(precision)) * time_to_go;
            if (precision==DoublePrecision() && t<1e-150)
                return mpfr_float("1e-150");
            else
                return t;
        }

        inline
        mpfr_float MinStepSizeForPrecision(unsigned precision, mpfr_float const& time_to_go, int safety_digits = 3)
        {
            return MinTimeForCurrentPrecision(precision, time_to_go, safety_digits);
        }
 

        inline
        double ArithmeticCost(unsigned precision)
        {
            if (precision==DoublePrecision())
                return 1;
            else
                return 10.35 + 0.13 * precision;
        }
 
 
 

        inline
        mpfr_float StepsizeSatisfyingCriterionB(unsigned precision,
                                        unsigned digits_B,
                                        unsigned num_newton_iterations,
                                        unsigned predictor_order = 0)
        {
            return pow(mpfr_float(10), -( digits_B - precision )*num_newton_iterations/(predictor_order+1));
        }
        

        inline
        unsigned MinDigitsForLogOfStepsize(mpfr_float const& log_of_stepsize, mpfr_float const& time_to_go, unsigned safety_digits = 3)
        {
            return (ceil(log_of_stepsize) + ceil(log10(abs(time_to_go))) + safety_digits).convert_to<unsigned>();
        }

        inline
        unsigned MinDigitsForStepsizeInterval(mpfr_float const& min_stepsize, mpfr_float const& max_stepsize, mpfr_float const& time_to_go)
        {
            return max(MinDigitsForLogOfStepsize(-log10(min_stepsize),time_to_go),   
                       MinDigitsForLogOfStepsize(-log10(max_stepsize),time_to_go));
        }
        

        template<typename RealT>
        void MinimizeTrackingCost(unsigned & new_precision, RealT & new_stepsize, 
                          unsigned min_precision, RealT const& min_stepsize,
                          unsigned max_precision, RealT const& max_stepsize,
                          unsigned digits_B,
                          unsigned num_newton_iterations,
                          unsigned predictor_order = 0)
        {
            double min_cost = Eigen::NumTraits<double>::highest();
            new_precision = MaxPrecisionAllowed()+1; // initialize to an impossible value.
            new_stepsize = min_stepsize; // initialize to minimum permitted step size.
 
            auto minimizer_routine = 
                [&min_cost, &new_stepsize, &new_precision, &digits_B, num_newton_iterations, predictor_order, max_stepsize](unsigned p)
                {
                    RealT candidate_stepsize = min(StepsizeSatisfyingCriterionB(p, digits_B, num_newton_iterations, predictor_order),
                                                  max_stepsize);
                    using std::abs;
                    double current_cost = ArithmeticCost(p) / abs(double(candidate_stepsize));
 
                    if (current_cost < min_cost)
                    {
                        min_cost = current_cost;
                        new_stepsize = candidate_stepsize;
                        new_precision = p;
                    }
                };
 
            unsigned lowest_mp_precision_to_test = min_precision;
 
            if (min_precision<=DoublePrecision())
                minimizer_routine(DoublePrecision());           
 
 
            if (lowest_mp_precision_to_test < LowestMultiplePrecision())
                lowest_mp_precision_to_test = LowestMultiplePrecision();
            else
                lowest_mp_precision_to_test = (lowest_mp_precision_to_test/PrecisionIncrement()) * PrecisionIncrement(); // use integer arithmetic to round.
 
 
            if (max_precision < LowestMultiplePrecision())
                max_precision = LowestMultiplePrecision();
            else
                max_precision = (max_precision/PrecisionIncrement()) * PrecisionIncrement(); // use integer arithmetic to round.
 
            for (unsigned p = lowest_mp_precision_to_test; p <= max_precision; p+=PrecisionIncrement())
                minimizer_routine(p);
        }
 
 
 
 
 
 
 

        class AMPTracker : public Tracker<AMPTracker>
        {
            friend class Tracker<AMPTracker>;
        public:
            
            typedef Tracker<AMPTracker> Base;
            typedef typename TrackerTraits<AMPTracker>::EventEmitterType EmitterType;
 
            enum UpsampleRefinementOption
            {
               upsample_refine_off  = 0,
               upsample_refine_on   = 1
            };
        

            AMPTracker(class System const& sys) : Tracker(sys), current_precision_(DefaultPrecision())
            {   
                Set<PrecConf>(AMPConfigFrom(sys));
            }
            

            void PrecisionSetup(AdaptiveMultiplePrecisionConfig const& AMP_config)
            {
                Set<PrecConf>(AMP_config);
            }
 
 
            const unsigned GetCurrentPrecision() const
            {
                return current_precision_;
            }
 
            
            void PrecisionPreservation(bool should_preseve_precision)
            {
                preserve_precision_ = should_preseve_precision;
            }
 
            
            virtual ~AMPTracker() = default;
 
 
            Vec<mpfr_complex> CurrentPoint() const override
            {
                if (this->CurrentPrecision()==DoublePrecision())
                {
                    const auto& curr_vector = std::get<Vec<dbl>>(this->current_space_);
                    Vec<mpfr_complex> returnme(NumVariables());
                    for (unsigned ii = 0; ii < NumVariables(); ++ii)
                    {
                        returnme(ii) = mpfr_complex(curr_vector(ii));
                    }
                    return returnme;
                }
                else
                    return std::get<Vec<mpfr_complex>>(this->current_space_);
            }
 
            
 
        private:

            SuccessCode TrackerLoopInitialization(mpfr_complex const& start_time,
                                           mpfr_complex const& end_time,
                                           Vec<mpfr_complex> const& start_point) const override
            {
                #ifndef BERTINI_DISABLE_ASSERTS
                assert(
                        (!preserve_precision_ 
                         || 
                         -log10(tracking_tolerance_) <= start_point(0).precision())
                         && "when tracking a path, either preservation of precision must be turned off (so precision can be higher at the end of tracking), or the initial precision must be high enough to support the resulting points to the desired tolerance"
                         );
                #endif
 
                NotifyObservers(Initializing<AMPTracker,mpfr_complex>(*this,start_time, end_time, start_point));
 
                initial_precision_ = Precision(start_point(0));
                DefaultPrecision(initial_precision_);
                // set up the master current time and the current step size
                
                current_time_ = start_time;
                current_time_.precision(initial_precision_);
 
                
                endtime_highest_precision_ = end_time;
                endtime_highest_precision_.precision(initial_precision_);
 
                endtime_ = end_time;
                endtime_.precision(initial_precision_);
 
                current_stepsize_.precision(initial_precision_);
                if (reinitialize_stepsize_)
                {
                    mpfr_float segment_length = abs(start_time-end_time)/Get<Stepping>().min_num_steps;
                    SetStepSize(min(NumTraits<mpfr_float>::FromRational(Get<Stepping>().initial_step_size, current_precision_),segment_length));
                }
 
                // populate the current space value with the start point, in appropriate precision
                if (initial_precision_==DoublePrecision())
                    MultipleToDouble(start_point);
                else
                    MultipleToMultiple(initial_precision_, start_point);
                
                ChangePrecision<upsample_refine_off>(initial_precision_);
                
                ResetCounters();
 
                ChangePrecision<upsample_refine_off>(start_point(0).precision());
 
 
                auto initial_refinement_code = InitialRefinement();
 
                #ifndef BERTINI_DISABLE_ASSERTS
                if (initial_precision_==DoublePrecision()){
                    PrecisionSanityCheck<dbl>();
                }
                else{
                    PrecisionSanityCheck<mpfr_complex>();
                }
                #endif
 
                return initial_refinement_code;
            }
 
 
 
 
 
            void ResetCounters() const override
            {
                Tracker::ResetCountersBase();
                num_precision_decreases_ = 0;
                num_successful_steps_since_stepsize_increase_ = 0;
                num_successful_steps_since_precision_decrease_ = 0;
                // initialize to the frequency so guaranteed to compute it the first try    
                num_steps_since_last_condition_number_computation_ = this->Get<Stepping>().frequency_of_CN_estimation;
            }

            SuccessCode InitialRefinement() const
            {
                SuccessCode initial_refinement_code = RefineStoredPoint();
                if (initial_refinement_code!=SuccessCode::Success)
                {
                    do {
                        if (current_precision_ > Get<PrecConf>().maximum_precision)
                        {
                            NotifyObservers(SingularStartPoint<EmitterType>(*this));
                            return SuccessCode::SingularStartPoint;
                        }
 
                        if (current_precision_==DoublePrecision())
                            initial_refinement_code = ChangePrecision<upsample_refine_on>(LowestMultiplePrecision());
                        else
                            initial_refinement_code = ChangePrecision<upsample_refine_on>(current_precision_+PrecisionIncrement());
                    }
                    while (initial_refinement_code!=SuccessCode::Success);
                }
                return SuccessCode::Success;
            }       
 
 
 

            SuccessCode PreIterationCheck() const override
            {
                if (num_successful_steps_taken_ >= Get<Stepping>().max_num_steps)
                    return SuccessCode::MaxNumStepsTaken;
                if (current_stepsize_ < Get<Stepping>().min_step_size)
                    return SuccessCode::MinStepSizeReached;
                if (current_precision_ > Get<PrecConf>().maximum_precision)
                    return SuccessCode::MaxPrecisionReached;
 
                return SuccessCode::Success;
            }
 
 
 
            
            
 
 
            void PostTrackCleanup() const override
            {
                if (preserve_precision_)
                    ChangePrecision(initial_precision_);
                NotifyObservers(TrackingEnded<EmitterType>(*this));
            }

            void CopyFinalSolution(Vec<mpfr_complex> & solution_at_endtime) const override
            {
 
                // the current precision is the precision of the output solution point.
                if (current_precision_==DoublePrecision())
                {
                    unsigned num_vars = GetSystem().NumVariables();
                    solution_at_endtime.resize(num_vars);
                    for (unsigned ii=0; ii<num_vars; ii++)
                        solution_at_endtime(ii) = mpfr_complex(std::get<Vec<dbl> >(current_space_)(ii));
                }
                else
                {
                    unsigned num_vars = GetSystem().NumVariables();
                    solution_at_endtime.resize(num_vars);
                    for (unsigned ii=0; ii<num_vars; ii++)
                    {
                        solution_at_endtime(ii) = std::get<Vec<mpfr_complex> >(current_space_)(ii);
                        solution_at_endtime(ii).precision(current_precision_);
                    }
                }
            }
 
 
 

            SuccessCode TrackerIteration() const override
            {
                if (current_precision_==DoublePrecision())
                    return TrackerIteration<dbl>();
                else
                    return TrackerIteration<mpfr_complex>();
            }
 

            template <typename ComplexType>
            SuccessCode TrackerIteration() const // not an override, because it is templated
            {   
 
                using RealType = typename Eigen::NumTraits<ComplexType>::Real;
 
                #ifndef BERTINI_DISABLE_ASSERTS
                assert(PrecisionSanityCheck<ComplexType>() && "precision sanity check failed.  some internal variable is not in correct precision");
                #endif
 
                NotifyObservers(NewStep<EmitterType>(*this));
 
                Vec<ComplexType>& predicted_space = std::get<Vec<ComplexType> >(temporary_space_); // this will be populated in the Predict step
                Vec<ComplexType>& current_space = std::get<Vec<ComplexType> >(current_space_); // the thing we ultimately wish to update
                ComplexType current_time = ComplexType(current_time_);
                ComplexType delta_t = ComplexType(delta_t_);
 
                #ifndef BERTINI_DISABLE_ASSERTS
                PrecisionSanityCheck<ComplexType>();
                #endif
 
                SuccessCode predictor_code = Predict<ComplexType, RealType>(predicted_space, current_space, current_time, delta_t);
                if (predictor_code==SuccessCode::MatrixSolveFailureFirstPartOfPrediction)
                {
                    NotifyObservers(FirstStepPredictorMatrixSolveFailure<EmitterType>(*this));
                    InitialMatrixSolveError();
                    return predictor_code;
                }
                else if (predictor_code==SuccessCode::MatrixSolveFailure)
                {
                    NotifyObservers(PredictorMatrixSolveFailure<EmitterType>(*this));
                    NewtonConvergenceError();// decrease stepsize, and adjust precision as necessary
                    return predictor_code;
                }   
                else if (predictor_code==SuccessCode::HigherPrecisionNecessary)
                {   
                    NotifyObservers(PredictorHigherPrecisionNecessary<EmitterType>(*this));
                    AMPCriterionError<ComplexType>();
                    return predictor_code;
                }
 
 
                NotifyObservers(SuccessfulPredict<AMPTracker, ComplexType>(*this, predicted_space));
 
                Vec<ComplexType>& tentative_next_space = std::get<Vec<ComplexType> >(tentative_space_); // this will be populated in the Correct step
 
                ComplexType tentative_next_time = current_time + delta_t;
 
                SuccessCode corrector_code = Correct<ComplexType, RealType>(tentative_next_space,
                                                     predicted_space,
                                                     tentative_next_time);
 
                #ifndef BERTINI_DISABLE_ASSERTS
                assert(PrecisionSanityCheck<ComplexType>() && "precision sanity check failed.  some internal variable is not in correct precision");
                #endif
 
                if (corrector_code==SuccessCode::MatrixSolveFailure || corrector_code==SuccessCode::FailedToConverge)
                {
                    NotifyObservers(CorrectorMatrixSolveFailure<EmitterType>(*this));
                    NewtonConvergenceError();
                    return corrector_code;
                }
                else if (corrector_code == SuccessCode::HigherPrecisionNecessary)
                {
                    NotifyObservers(CorrectorHigherPrecisionNecessary<EmitterType>(*this));
                    AMPCriterionError<ComplexType>();
                    return corrector_code;
                }
                else if (corrector_code == SuccessCode::GoingToInfinity)
                {
                    // there is no corrective action possible...
                    return corrector_code;
                }
 
                NotifyObservers(SuccessfulCorrect<AMPTracker, ComplexType>(*this, tentative_next_space));
 
                // copy the tentative vector into the current space vector;
                current_space = tentative_next_space;
                return AdjustAMPStepSuccess<ComplexType>();
            }
 

            SuccessCode CheckGoingToInfinity() const override
            {
                if (current_precision_ == DoublePrecision())
                    return Base::CheckGoingToInfinity<dbl>();
                else
                    return Base::CheckGoingToInfinity<mpfr_complex>();
            }

            template <UpsampleRefinementOption refine_if_necessary = upsample_refine_off>
            SuccessCode UpdatePrecisionAndStepsize() const
            {
                SetStepSize(next_stepsize_);
                return ChangePrecision<refine_if_necessary>(next_precision_);
            }
 
 
 
 
 
 

            template <typename ComplexType>
            SuccessCode AdjustAMPStepSuccess() const
            {
                // TODO: think about why we consider reducing the stepsize?  this is despite documentation stating that it can only increase
                mpfr_float min_stepsize = current_stepsize_ * NumTraits<mpfr_float>::FromRational(Get<Stepping>().step_size_fail_factor, current_precision_);
                mpfr_float max_stepsize = min( current_stepsize_ * NumTraits<mpfr_float>::FromRational(Get<Stepping>().step_size_success_factor, current_precision_),  NumTraits<mpfr_float>::FromRational(Get<Stepping>().max_step_size, current_precision_));
 
 
                unsigned min_precision = MinRequiredPrecision_BCTol<ComplexType>();
                unsigned max_precision = max(min_precision,current_precision_);
 
                if (num_successful_steps_since_stepsize_increase_ < Get<Stepping>().consecutive_successful_steps_before_stepsize_increase)
                    max_stepsize = current_stepsize_; // disallow stepsize changing 
 
 
                if ( (num_successful_steps_since_precision_decrease_ < Get<PrecConf>().consecutive_successful_steps_before_precision_decrease)
                    ||
                    (num_precision_decreases_ >= Get<PrecConf>().max_num_precision_decreases))
                    min_precision = max(min_precision, current_precision_); // disallow precision changing 
 
 
                MinimizeTrackingCost(next_precision_, next_stepsize_, 
                            min_precision, min_stepsize,
                            max_precision, max_stepsize,
                            DigitsB<ComplexType>(),
                            Get<NewtonConfig>().max_num_newton_iterations,
                            predictor_order_);
 
 
                if ( (next_stepsize_ > current_stepsize_) || (next_precision_ < current_precision_) )
                    num_successful_steps_since_stepsize_increase_ = 0;
                else
                    ++num_successful_steps_since_stepsize_increase_;
                
 
                if (next_precision_ < current_precision_)
                { 
                    ++num_precision_decreases_;
                    num_successful_steps_since_precision_decrease_ = 0;
                }
                else
                    ++num_successful_steps_since_precision_decrease_;
 
                return UpdatePrecisionAndStepsize();
            }
 
 
 
            void InitialMatrixSolveError() const
            {
                next_stepsize_ = current_stepsize_;
 
                if (current_precision_==DoublePrecision())
                    next_precision_ = LowestMultiplePrecision();
                else
                    next_precision_ = current_precision_+(1+num_consecutive_failed_steps_) * PrecisionIncrement();
 
                UpdatePrecisionAndStepsize();
            }

            void NewtonConvergenceError() const
            {
                next_precision_ = current_precision_;
 
                next_stepsize_ = mpfr_float(Get<Stepping>().step_size_fail_factor, CurrentPrecision())*current_stepsize_;
                while (next_stepsize_ < MinStepSizeForPrecision(next_precision_, abs(current_time_ - endtime_)))
                {
                    if (next_precision_==DoublePrecision())
                        next_precision_=LowestMultiplePrecision();
                    else
                        next_precision_+=PrecisionIncrement();
                }
 
                UpdatePrecisionAndStepsize();
            }
 
 
 
 
 
 

            template<typename ComplexType>
            void AMPCriterionError() const
            {   
                using RealT = typename Eigen::NumTraits<ComplexType>::Real;
 
                unsigned min_next_precision; // sure, i could use a trigraph here, but it'd be terrible
                if (current_precision_==DoublePrecision())
                    min_next_precision = LowestMultiplePrecision(); // precision increases
                else
                    min_next_precision = current_precision_ + (1+num_consecutive_failed_steps_)*PrecisionIncrement(); // precision increases
 
 
                mpfr_float min_stepsize = MinStepSizeForPrecision(current_precision_, abs(current_time_ - endtime_));
                mpfr_float max_stepsize = current_stepsize_ * NumTraits<mpfr_float>::FromRational(Get<Stepping>().step_size_fail_factor,current_precision_);  // Stepsize decreases.
 
                if (min_stepsize > max_stepsize)
                {
                    // stepsizes are incompatible, must increase precision
                    next_precision_ = min_next_precision;
                    // decrease stepsize somewhat less than the fail factor
                    next_stepsize_ = max(current_stepsize_ * (1+NumTraits<mpfr_float>::FromRational(Get<Stepping>().step_size_fail_factor,current_precision_))/2, min_stepsize);
                }
                else
                {
                    unsigned digits_B = DigitsB<ComplexType>();
 
                    unsigned min_precision = max(min_next_precision,
                                                 digits_B,
                                                 DigitsC<ComplexType>(),
                                                 MinDigitsForStepsizeInterval(min_stepsize, max_stepsize, abs(current_time_ - endtime_)),
                                                 digits_final_
                                                 );
                    
                        unsigned a = ceil(digits_B - (predictor_order_+1)* -log10(max_stepsize)/Get<NewtonConfig>().max_num_newton_iterations).convert_to<unsigned>();
 
                    unsigned max_precision = max(min_precision, a);
 
                    MinimizeTrackingCost(next_precision_, next_stepsize_, 
                            min_precision, min_stepsize,
                            Get<PrecConf>().maximum_precision, max_stepsize,
                            digits_B,
                            Get<NewtonConfig>().max_num_newton_iterations,
                            predictor_order_);
                }
 
                UpdatePrecisionAndStepsize();
            }
 
 
 
 

            template<typename ComplexType>
            NumErrorT B_RHS() const
            {   
                return max(amp::CriterionBRHS(this->norm_J_, 
                                              this->norm_J_inverse_, 
                                              Get<NewtonConfig>().max_num_newton_iterations, 
                                              tracking_tolerance_, 
                                              this->size_proportion_, 
                                              Get<PrecConf>()), NumErrorT(0));
            }
 
 
 

            template<typename ComplexType>
            unsigned DigitsB() const
            {   
                return unsigned(B_RHS<ComplexType>());
            }
 
 
 
            

            template<typename ComplexType>
            NumErrorT C_RHS() const
            {   
                return max(amp::CriterionCRHS(this->norm_J_inverse_, 
                                              NumErrorT(std::get<Vec<ComplexType> > (current_space_).norm()), 
                                              tracking_tolerance_, 
                                              Get<PrecConf>()), NumErrorT(0));
            }
 
 

            template<typename ComplexType>
            unsigned DigitsC() const
            {   
                return unsigned(C_RHS<ComplexType>());
            }
 
 
 
 

            template <typename ComplexType>
            unsigned MinRequiredPrecision_BCTol() const
            {
                return max(DigitsB<ComplexType>(), 
                           DigitsC<ComplexType>(), 
                           digits_tracking_tolerance_, 
                           DoublePrecision()); 
            }
 
 
 
 
 
 

            //
            //  overrides for counter adjustment after a TrackerIteration()
            //

            void OnStepSuccess() const override
            {
                Tracker::IncrementBaseCountersSuccess();
                NotifyObservers(SuccessfulStep<EmitterType>(*this));
            }

            void OnStepFail() const override
            {
                Tracker::IncrementBaseCountersFail();
                num_successful_steps_since_precision_decrease_ = 0;
                num_successful_steps_since_stepsize_increase_ = 0;
                NotifyObservers(FailedStep<EmitterType>(*this));
            }
 
 
 
            void OnInfiniteTruncation() const override
            {
                NotifyObservers(InfinitePathTruncation<EmitterType>(*this));
            }
 

            //
            //       Predict and Correct functions
            //
 
 
            
 

            template<typename ComplexType, typename RealType>
            SuccessCode Predict(Vec<ComplexType> & predicted_space, 
                                const Vec<ComplexType>& current_space,
                                ComplexType const& current_time, ComplexType const& delta_t) const
            {
                
                                
                
                
                static_assert(std::is_same< typename Eigen::NumTraits<RealType>::Real, 
                                        typename Eigen::NumTraits<ComplexType>::Real>::value,
                                        "underlying complex type and the type for comparisons must match");
 
 
                if (predictor_->HasErrorEstimate())
                    return predictor_->Predict(predicted_space,
                                    this->error_estimate_,
                                    this->size_proportion_,
                                    this->norm_J_,
                                    this->norm_J_inverse_,
                                    tracked_system_,
                                    current_space, current_time, 
                                    delta_t,
                                    this->condition_number_estimate_,
                                    num_steps_since_last_condition_number_computation_, 
                                    Get<Stepping>().frequency_of_CN_estimation, 
                                    tracking_tolerance_,
                                    Get<PrecConf>());
                else
                    return predictor_->Predict(predicted_space,
                                    this->size_proportion_,
                                    this->norm_J_,
                                    this->norm_J_inverse_,
                                    tracked_system_,
                                    current_space, current_time, 
                                    delta_t,
                                    this->condition_number_estimate_,
                                    num_steps_since_last_condition_number_computation_, 
                                    Get<Stepping>().frequency_of_CN_estimation, 
                                    tracking_tolerance_,
                                    Get<PrecConf>());
            }
 
 

            template<typename ComplexType, typename RealType>
            SuccessCode Correct(Vec<ComplexType> & corrected_space, 
                                Vec<ComplexType> const& current_space, 
                                ComplexType const& current_time) const
            {
                static_assert(std::is_same< typename Eigen::NumTraits<RealType>::Real, 
                                        typename Eigen::NumTraits<ComplexType>::Real>::value,
                                        "underlying complex type and the type for comparisons must match");
 
 
 
 
                return corrector_->Correct(corrected_space,
                                    this->norm_delta_z_,
                                    this->norm_J_,
                                    this->norm_J_inverse_,
                                    this->condition_number_estimate_,
                                    tracked_system_,
                                    current_space,
                                    current_time,
                                    tracking_tolerance_,
                                    Get<NewtonConfig>().min_num_newton_iterations,
                                    Get<NewtonConfig>().max_num_newton_iterations,
                                    Get<PrecConf>());
            }
 
 

            SuccessCode RefineStoredPoint() const
            {
                SuccessCode code;
                if (current_precision_==DoublePrecision())
                {
                    code = RefineImpl<dbl>(std::get<Vec<dbl> >(temporary_space_),std::get<Vec<dbl> >(current_space_), dbl(current_time_));
                    if (code == SuccessCode::Success)
                        std::get<Vec<dbl> >(current_space_) = std::get<Vec<dbl> >(temporary_space_);
                }
                else
                {
                    code = RefineImpl<mpfr_complex>(std::get<Vec<mpfr_complex> >(temporary_space_),std::get<Vec<mpfr_complex> >(current_space_), current_time_);
                    if (code == SuccessCode::Success)
                        std::get<Vec<mpfr_complex> >(current_space_) = std::get<Vec<mpfr_complex> >(temporary_space_);
                }
                return code;
            }
 
 

            template <typename ComplexType>
            SuccessCode RefineImpl(Vec<ComplexType> & new_space,
                                Vec<ComplexType> const& start_point, ComplexType const& current_time) const
            {
                using RealType = typename Eigen::NumTraits<ComplexType>::Real;
                
                static_assert(std::is_same< typename Eigen::NumTraits<RealType>::Real, 
                                        typename Eigen::NumTraits<ComplexType>::Real>::value,
                                        "underlying complex type and the type for comparisons must match");
 
                auto target_precision = Precision(current_time);
                assert(Precision(start_point)==target_precision);
                ChangePrecision(target_precision);
                Precision(new_space,target_precision);
                
 
 
 
                return corrector_->Correct(new_space,
                                           this->norm_delta_z_,
                                           this->norm_J_,
                                           this->norm_J_inverse_,
                                           this->condition_number_estimate_,
                                           tracked_system_,
                                           start_point,
                                           current_time,
                                           tracking_tolerance_,
                                           Get<NewtonConfig>().min_num_newton_iterations,
                                           Get<NewtonConfig>().max_num_newton_iterations,
                                           Get<PrecConf>());
            }
 
 
 
 
 

            template <typename ComplexType>
            SuccessCode RefineImpl(Vec<ComplexType> & new_space,
                                Vec<ComplexType> const& start_point, ComplexType const& current_time,
                                NumErrorT const& tolerance, unsigned max_iterations) const
            {
                auto target_precision = Precision(current_time);
                assert(Precision(start_point)==target_precision);
                ChangePrecision(target_precision);
                Precision(new_space,target_precision);
 
                return corrector_->Correct(new_space,
                                        this->norm_delta_z_,
                                        this->norm_J_,
                                        this->norm_J_inverse_,
                                        this->condition_number_estimate_,
                                        tracked_system_,
                                        start_point,
                                        current_time, 
                                        tolerance,
                                        1,
                                        max_iterations,
                                        Get<PrecConf>());
            }
 
 
 
 
 
 
 
 
 
 
 
 
 

            //
            //  Functions for converting between precision types
            //
 
        public:
            template <UpsampleRefinementOption refine_if_necessary = upsample_refine_off>
            SuccessCode ChangePrecision(unsigned new_precision) const
            {
                if (new_precision==current_precision_) // no op
                    return SuccessCode::Success;
 
                NotifyObservers(PrecisionChanged<EmitterType>(*this,current_precision_,new_precision));
                
 
                bool upsampling_needed = new_precision > current_precision_;
                // reset the counter for estimating the condition number.  
                num_steps_since_last_condition_number_computation_ = this->Get<Stepping>().frequency_of_CN_estimation;
 
                if (new_precision==DoublePrecision() && current_precision_>DoublePrecision())
                {
                    // convert from multiple precision to double precision
                    MultipleToDouble();
                }
                else if(new_precision > DoublePrecision() && current_precision_ == DoublePrecision())
                {
                    // convert from double to multiple precision
                    DoubleToMultiple(new_precision);
                    #ifndef BERTINI_DISABLE_ASSERTS
                    assert(PrecisionSanityCheck<mpfr_complex>() && "precision sanity check failed.  some internal variable is not in correct precision");
                    #endif
                }
                else
                {
                    MultipleToMultiple(new_precision);
                    #ifndef BERTINI_DISABLE_ASSERTS
                    assert(PrecisionSanityCheck<mpfr_complex>() && "precision sanity check failed.  some internal variable is not in correct precision");
                    #endif
                }
 
                if (refine_if_necessary && upsampling_needed)
                    return RefineStoredPoint();
                else
                    return SuccessCode::Success;
            }
 
 
 
 
 
            private:
            
            void DoubleToDouble(Vec<dbl> const& source_point) const
            {   
                #ifndef BERTINI_DISABLE_ASSERTS
                assert(source_point.size() == GetSystem().NumVariables() && "source point for converting to multiple precision is not the same size as the number of variables in the system being solved.");
                #endif
 
                current_precision_ = DoublePrecision();
                DefaultPrecision(DoublePrecision());
 
                GetSystem().precision(16);
 
                std::get<Vec<dbl> >(current_space_) = source_point;
            }

            void DoubleToDouble() const
            {
                DoubleToDouble(std::get<Vec<dbl> >(current_space_));
            }
            
 

            void MultipleToDouble(Vec<mpfr_complex> const& source_point) const
            {   
                #ifndef BERTINI_DISABLE_ASSERTS
                assert(source_point.size() == GetSystem().NumVariables() && "source point for converting to multiple precision is not the same size as the number of variables in the system being solved.");
                #endif
 
                AdjustCurrentPrecision(DoublePrecision());
 
                GetSystem().precision(DoublePrecision()); 
 
                // copy the current space in.
                if (std::get<Vec<dbl> >(current_space_).size()!=source_point.size())
                    std::get<Vec<dbl> >(current_space_).resize(source_point.size());
 
                for (unsigned ii=0; ii<source_point.size(); ii++)
                    std::get<Vec<dbl> >(current_space_)(ii) = dbl(source_point(ii));
 
                endtime_.precision(DoublePrecision()); // i question this one  2021-04-12
            }

            void MultipleToDouble() const
            {
                MultipleToDouble(std::get<Vec<mpfr_complex> >(current_space_));
            }
 
 

            void DoubleToMultiple(unsigned new_precision, Vec<dbl> const& source_point) const
            {   
                #ifndef BERTINI_DISABLE_ASSERTS
                assert(source_point.size() == GetSystem().NumVariables() && "source point for converting to multiple precision is not the same size as the number of variables in the system being solved.");
                assert(new_precision > DoublePrecision() && "must convert to precision higher than DoublePrecision when converting to multiple precision");
                #endif
                
                AdjustCurrentPrecision(new_precision);
                CopyToCurrentSpace(source_point);
                AdjustInternalsPrecision(new_precision);
                AdjustTemporariesPrecision(new_precision);
 
                #ifndef BERTINI_DISABLE_ASSERTS
                PrecisionSanityCheck<mpfr_complex>();
                #endif
            }
 
 

            void DoubleToMultiple(unsigned new_precision) const
            {
                DoubleToMultiple( new_precision, std::get<Vec<dbl> >(current_space_));
            }
 
 
 
 

            void MultipleToMultiple(unsigned new_precision, Vec<mpfr_complex> const& source_point) const
            {   
                #ifndef BERTINI_DISABLE_ASSERTS
                assert(source_point.size() == GetSystem().NumVariables() && "source point for converting to multiple precision is not the same size as the number of variables in the system being solved.");
                assert(new_precision > DoublePrecision() && "must convert to precision higher than DoublePrecision when converting to multiple precision");
                #endif
 
                AdjustCurrentPrecision(new_precision);
                CopyToCurrentSpace(source_point);
                AdjustInternalsPrecision(new_precision);
                AdjustTemporariesPrecision(new_precision);
 
                #ifndef BERTINI_DISABLE_ASSERTS
                PrecisionSanityCheck<mpfr_complex>();
                #endif
            }
 

            void MultipleToMultiple(unsigned new_precision) const
            {
                MultipleToMultiple( new_precision, std::get<Vec<mpfr_complex> >(current_space_));
            }
 
 
 
 
 
            void AdjustCurrentPrecision(unsigned new_precision) const
            {
                previous_precision_ = current_precision_;
                current_precision_ = new_precision;
                DefaultPrecision(new_precision);
            }
 
            
            void CopyToCurrentSpace(Vec<dbl> const& source_point) const
            {
                auto& space = std::get<Vec<mpfr_complex> >(current_space_);
                if (space.size()!=source_point.size())
                    space.resize(source_point.size());
                for (unsigned ii=0; ii<source_point.size(); ii++)
                    space(ii) = mpfr_complex(source_point(ii));
            }
 
            void CopyToCurrentSpace(Vec<mpfr_complex> const& source_point) const
            {
                auto& space = std::get<Vec<mpfr_complex> >(current_space_);
                if (space.size()!=source_point.size())
                    space.resize(source_point.size());
                space = source_point;
            }
 
            void AdjustInternalsPrecision(unsigned new_precision) const
            {
                GetSystem().precision(new_precision);
                predictor_->ChangePrecision(new_precision);
                corrector_->ChangePrecision(new_precision);
 
                endtime_ = endtime_highest_precision_;
 
                endtime_.precision(new_precision);
                current_stepsize_.precision(new_precision);
                delta_t_.precision(new_precision);
                current_time_.precision(new_precision);
 
                Precision(std::get<Vec<mpfr_complex> >(current_space_),new_precision);
            }

            void AdjustTemporariesPrecision(unsigned new_precision) const
            {
                const auto num_vars = GetSystem().NumVariables();
 
                //  the current_space value is adjusted in the appropriate ChangePrecision function
                std::get<Vec<mpfr_complex> >(tentative_space_).resize(num_vars);
                Precision(std::get<Vec<mpfr_complex> >(tentative_space_), new_precision);
 
                std::get<Vec<mpfr_complex> >(temporary_space_).resize(num_vars);
                Precision(std::get<Vec<mpfr_complex> >(temporary_space_), new_precision);
            }
 
 
 
            
            template<typename ComplexType>
            bool PrecisionSanityCheck() const
            {   
 
                if constexpr (std::is_same<ComplexType, dbl>::value){
                    return true;
                }
 
                if constexpr (std::is_same<ComplexType, mpfr_complex>::value){
                    assert(DefaultPrecision()==current_precision_ && "current precision differs from the default precision");
                    assert(GetSystem().precision() == current_precision_ && "tracked system is out of precision");
                    
                    assert(std::get<Vec<mpfr_complex> >(current_space_)(0).precision() == current_precision_ && "current space out of precision");
                    assert(std::get<Vec<mpfr_complex> >(tentative_space_)(0).precision() == current_precision_ && "tentative space out of precision");
                    assert(std::get<Vec<mpfr_complex> >(temporary_space_)(0).precision() == current_precision_ && "temporary space out of precision");
                    assert(Precision(current_stepsize_) == current_precision_ && "current_stepsize_ out of precision");
                    assert(Precision(delta_t_) == current_precision_ && "delta_t_ out of precision");
                    assert(Precision(endtime_) == current_precision_ && "endtime_ out of precision");
                    assert(Precision(current_time_) == current_precision_ && "current_time_ out of precision");
                    assert(current_precision_ <= MaxPrecisionAllowed() && "current_precision_ exceeds max precision");
                    assert(predictor_->precision() == current_precision_ && "predictor_ out of precision");
                    return  true;
                }
 
            }
 
 
 

            //
 
 
 
 

            // state variables
            bool preserve_precision_ = false; 
 
            mutable unsigned previous_precision_; 
            mutable unsigned current_precision_; 
            mutable unsigned next_precision_; 
            mutable unsigned num_precision_decreases_; 
            mutable unsigned initial_precision_; 
            mutable unsigned num_successful_steps_since_precision_decrease_; 
 
            mutable mpfr_complex endtime_highest_precision_;
 
        public:
 
            unsigned CurrentPrecision() const override
            {
                return current_precision_;
            }
        }; // re: class Tracker
 
    } // namespace tracking
} // namespace bertini
 
 
#endif