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) 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_AMP_TRACKER_HPP
#define BERTINI_AMP_TRACKER_HPP
 
#include <assert.h>
#include <optional>
 
#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
        real_mp MinTimeForCurrentPrecision(unsigned precision, real_mp const& time_to_go, int safety_digits = 3)
        {
            real_mp t = pow( real_mp(10), safety_digits-long(precision)) * time_to_go;
            if (precision==DoublePrecision() && t<1e-150)
                return real_mp("1e-150");
            else
                return t;
        }

        inline
        real_mp MinStepSizeForPrecision(unsigned precision, real_mp 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 101.47 + 1.59 * precision;
        }
 
 
 

        template<typename RealT>
        inline RealT StepsizeSatisfyingCriterionB(unsigned precision,
                                        unsigned digits_B,
                                        unsigned num_newton_iterations,
                                        unsigned predictor_order = 0)
        {
            return pow(RealT(10), -(static_cast<int>(digits_B) - static_cast<int>(precision))
                                   * static_cast<int>(num_newton_iterations) / (predictor_order+1.0));
        }
 

        inline
        unsigned MinDigitsForLogOfStepsize(real_mp const& log_of_stepsize, real_mp 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(real_mp const& min_stepsize, real_mp const& max_stepsize, real_mp const& time_to_go)
        {
            return max(MinDigitsForLogOfStepsize(-log10(min_stepsize),time_to_go),
                       MinDigitsForLogOfStepsize(-log10(max_stepsize),time_to_go));
        }

        inline
        unsigned ExtraDigitsBeforePrecisionDecrease(unsigned current_precision)
        {
            if (current_precision <= DoublePrecision())
                return 0;
            return (current_precision / PrecisionIncrement()) * 2;
        }

        inline
        void ApplyPrecisionDecreaseMargin(int & digits_required, unsigned current_digits, unsigned extra)
        {
            if (digits_required < static_cast<int>(current_digits))
            {
                digits_required += static_cast<int>(extra);
                if (digits_required > static_cast<int>(current_digits))
                    digits_required = static_cast<int>(current_digits);
            }
        }
 

        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();
 
            unsigned minimizing_precision = 0; // initialize to an impossible value.
 
            // a few casts so that we can work in double precision, because doing this in multiprec sucks
            double min_stepsize_lowprec = static_cast<double>(min_stepsize);
            double max_stepsize_lowprec = static_cast<double>(max_stepsize);
 
            auto minimizer_routine =
                [&min_cost, &minimizing_precision, digits_B, num_newton_iterations, predictor_order, min_stepsize_lowprec, max_stepsize_lowprec](unsigned p)
                {
                    double criterion_b_stepsize = StepsizeSatisfyingCriterionB<double>(p, digits_B, num_newton_iterations, predictor_order);
                    if (criterion_b_stepsize < min_stepsize_lowprec)
                        return; // precision too low to satisfy min_stepsize constraint
                    double candidate_stepsize = min(criterion_b_stepsize, max_stepsize_lowprec);
                    using std::abs;
                    double current_cost = ArithmeticCost(p) / abs(candidate_stepsize);
 
                    if (current_cost < min_cost)
                    {
                        min_cost = current_cost;
                        minimizing_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);
 
            if (minimizing_precision==0){
                throw std::runtime_error("MinimizeTrackingCost failed to find a suitable stepsize and precision");
            }
 
            new_precision = minimizing_precision; // copy the computed value
            // next, because the above computed the new stepsize in double precision, which may be lowprec, we compute in full precision
            new_stepsize = max(
                               min(
                                    StepsizeSatisfyingCriterionB<RealT>(new_precision, digits_B, num_newton_iterations, predictor_order),
                                    max_stepsize
                                ),
                               min_stepsize
                               );
        }
 
 
 
 
 
 
 

        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())
            {
                // Deriving the criteria from the system needs a degree bound, and a system that is
                // not polynomial has none.  Rather than refuse to exist, keep the default config
                // -- which is internally consistent, if meaningless for such a system -- so the
                // caller can still reach PrecisionSetup and install bounds they can defend.  That
                // is the escape hatch the refusal in SetAMPConfigFrom names; without this guard it
                // would be unreachable, because the refusal would fire here first.
                if (sys.IsPolynomial())
                    Set<PrecConf>(AMPConfigFrom(sys));
            }
 

            void PrecisionSetup(AdaptiveMultiplePrecisionConfig const& AMP_config)
            {
                Set<PrecConf>(AMP_config);
            }
 

            unsigned GetCurrentPrecision() const
            {
                return current_precision_;
            }
 

            void PrecisionPreservation(bool should_preseve_precision)
            {
                preserve_precision_ = should_preseve_precision;
            }

            void SetStartPrecision(std::optional<unsigned> p = std::nullopt)
            {
                override_start_precision_ = p;
            }
 
 
            virtual ~AMPTracker() = default;
 

            Vec<complex_mp> CurrentPoint() const override
            {
                if (this->CurrentPrecision()==DoublePrecision())
                {
                    const auto& curr_vector = std::get<Vec<complex_dbl>>(this->current_space_);
                    Vec<complex_mp> returnme(NumVariables());
                    for (unsigned ii = 0; ii < NumVariables(); ++ii)
                    {
                        returnme(ii) = complex_mp(curr_vector(ii));
                    }
                    return returnme;
                }
                else
                    return std::get<Vec<complex_mp>>(this->current_space_);
            }
 
 
 
        private:

            SuccessCode TrackerLoopInitialization(complex_mp const& start_time,
                                           complex_mp const& end_time,
                                           Vec<complex_mp> const& start_point) const override
            {
                initial_precision_ = override_start_precision_.value_or(Precision(start_point(0)));
 
                #ifndef BERTINI_DISABLE_ASSERTS
                assert(
                        (!preserve_precision_
                         ||
                         -log10(this->TrackingTolerance()) <= initial_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,complex_mp>(*this,start_time, end_time, start_point));
                SetThreadPrecision(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_);
 
                // The first step is built at THIS track's start precision, and the precision is set
                // after the assignment, never before it.  current_precision_ still holds wherever the
                // previous track ended, and an assignment takes its source's precision
                // (preserve_related_precision), so building it at current_precision_ rounded the
                // exact initial_step_size at the previous track's precision.  A track that starts in
                // double then rounded that a second time, to a double, and the last bit of its first
                // step -- and so its whole step sequence -- depended on which path the tracker had
                // tracked before: thread scheduling made the threaded solve irreproducible, and path
                // order made a serial solve depend on it (#378).
                if (reinitialize_stepsize_)
                {
                    real_mp segment_length = abs(start_time-end_time)/Get<Stepping>().min_num_steps;
                    SetStepSize(min(real_mp(Get<Stepping>().initial_step_size, initial_precision_),segment_length));
                }
                current_stepsize_.precision(initial_precision_);
 
                // 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();
 
                auto initial_refinement_code = InitialRefinement();
 
                #ifndef BERTINI_DISABLE_ASSERTS
                if (initial_precision_==DoublePrecision()){
                    PrecisionSanityCheck<complex_dbl>();
                }
                else{
                    PrecisionSanityCheck<complex_mp>();
                }
                #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
            {
                // The budget counts EVERY step, not only the successful ones -- see the
                // same check in FixedPrecisionTracker and issue #410.  This is the guard
                // that bounds a path whose steps fail without ever advancing.
                if (NumTotalStepsTaken() >= 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 CopySolution(Vec<complex_mp> & solution_at_endtime) const override
            {
 
                // the current precision is the precision of the output solution point.
                if (current_precision_==DoublePrecision())
                {
                    unsigned num_vars = static_cast<unsigned>(GetSystem().NumVariables());
                    solution_at_endtime.resize(num_vars);
                    for (unsigned ii=0; ii<num_vars; ii++)
                        solution_at_endtime(ii) = complex_mp(std::get<Vec<complex_dbl> >(current_space_)(ii));
                }
                else
                {
                    unsigned num_vars = static_cast<unsigned>(GetSystem().NumVariables());
                    solution_at_endtime.resize(num_vars);
                    for (unsigned ii=0; ii<num_vars; ii++)
                    {
                        solution_at_endtime(ii) = std::get<Vec<complex_mp> >(current_space_)(ii);
                        solution_at_endtime(ii).precision(current_precision_);
                    }
                }
            }
 
 
 

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

            template <typename ComplexT>
            SuccessCode TrackerIteration() const // not an override, because it is templated
            {
 
                using RealT = typename Eigen::NumTraits<ComplexT>::Real;
 
                #ifndef BERTINI_DISABLE_ASSERTS
                assert(PrecisionSanityCheck<ComplexT>() && "precision sanity check failed.  some internal variable is not in correct precision");
                #endif
 
                NotifyObservers(NewStep<EmitterType>(*this));
 
                Vec<ComplexT>& predicted_space = std::get<Vec<ComplexT> >(temporary_space_); // this will be populated in the Predict step
                Vec<ComplexT>& current_space = std::get<Vec<ComplexT> >(current_space_); // the thing we ultimately wish to update
                ComplexT current_time = ComplexT(current_time_);
                ComplexT delta_t = ComplexT(delta_t_);
 
                #ifndef BERTINI_DISABLE_ASSERTS
                PrecisionSanityCheck<ComplexT>();
                #endif
 
                SuccessCode predictor_code = Predict<ComplexT, RealT>(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));
                    auto adj_code = AMPCriterionError<ComplexT>();
                    if (adj_code != SuccessCode::Success)
                        return adj_code;
                    return predictor_code;
                }
 
 
                NotifyObservers(SuccessfulPredict<AMPTracker, ComplexT>(*this, predicted_space));
 
                Vec<ComplexT>& tentative_next_space = std::get<Vec<ComplexT> >(tentative_space_); // this will be populated in the Correct step
 
                ComplexT tentative_next_time = current_time + delta_t;
 
                SuccessCode corrector_code = Correct<ComplexT, RealT>(tentative_next_space,
                                                     predicted_space,
                                                     tentative_next_time);
 
                #ifndef BERTINI_DISABLE_ASSERTS
                assert(PrecisionSanityCheck<ComplexT>() && "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));
                    auto adj_code = AMPCriterionError<ComplexT>();
                    if (adj_code != SuccessCode::Success)
                        return adj_code;
                    return corrector_code;
                }
                else if (corrector_code == SuccessCode::GoingToInfinity)
                {
                    // there is no corrective action possible...
                    return corrector_code;
                }
 
                NotifyObservers(SuccessfulCorrect<AMPTracker, ComplexT>(*this, tentative_next_space));
 
                // copy the tentative vector into the current space vector;
                current_space = tentative_next_space;
                return AdjustAMPStepSuccess<ComplexT>();
            }
 

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

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

            template <typename ComplexT>
            SuccessCode AdjustAMPStepSuccess() const
            {
                // This mirrors Bertini 1's AMP2_update: pick the cost-minimizing (precision, stepsize)
                // after a successful step.  Faithful-to-B1 points, two of which were previously wrong:
                //  (1) The rule-B precision requirement CREDITS THE STEP.  digits_B is the right-hand side
                //      of  P + (p+1)*xi/N > digits_B  (xi = -log10(stepsize)); the floor is therefore
                //      digits_B - (p+1)*(-log10(min_stepsize))/N, crediting the smallest allowed step, NOT
                //      the raw digits_B.  Without this credit a tiny step (large size_proportion in the
                //      roundoff regime) spuriously forced precision up.  B1: digits_B2 = ceil(P0 - eta_minSS/N).
                //  (2) Precision decrease (and stepsize increase) is gated by BOTH B1's StepsForIncrease
                //      consecutive-successful-steps counter AND the digits-margin hysteresis
                //      (ExtraDigitsBeforePrecisionDecrease).  The single StepsForIncrease setting
                //      (SteppingConfig::consecutive_successful_steps_before_stepsize_increase) governs both
                //      gates -- the old, duplicate AMP-config precision-decrease setting was removed.
                real_mp min_stepsize = current_stepsize_ * real_mp(Get<Stepping>().step_size_fail_factor, current_precision_);
                real_mp max_stepsize = min( current_stepsize_ * real_mp(Get<Stepping>().step_size_success_factor, current_precision_),  real_mp(Get<Stepping>().max_step_size, current_precision_));
 
                const unsigned steps_for_increase = Get<Stepping>().consecutive_successful_steps_before_stepsize_increase; // B1 StepsForIncrease
 
                if (num_successful_steps_since_stepsize_increase_ < steps_for_increase)
                    max_stepsize = current_stepsize_; // not enough successes in a row: disallow stepsize increase
 
                const unsigned N = Get<NewtonConfig>().max_num_newton_iterations;
                const unsigned curr_digits = current_precision_;
                const unsigned extra = ExtraDigitsBeforePrecisionDecrease(curr_digits);
 
                // rule B, step-credited (credit the smallest allowed step = the most credit)
                int digits_B = static_cast<int>(ceil( double(B_RHS<ComplexT>())
                                   - (predictor_order_ + 1.0) * double(-log10(min_stepsize)) / N ));
                int digits_C = static_cast<int>(DigitsC<ComplexT>());
                int digits_stepsize = static_cast<int>(MinDigitsForStepsizeInterval(min_stepsize, max_stepsize, abs(current_time_ - endtime_)));
 
                // each requirement may only lower precision if it clears the hysteresis margin
                ApplyPrecisionDecreaseMargin(digits_B, curr_digits, extra);
                ApplyPrecisionDecreaseMargin(digits_C, curr_digits, extra);
                ApplyPrecisionDecreaseMargin(digits_stepsize, curr_digits, extra);
 
                int min_digits = std::max({ static_cast<int>(this->DigitsTrackingTolerance()),
                                            static_cast<int>(digits_final_),
                                            digits_B, digits_C, digits_stepsize,
                                            static_cast<int>(DoublePrecision()) });
                unsigned min_precision = static_cast<unsigned>(min_digits); // >= DoublePrecision(), so non-negative
 
                // StepsForIncrease gate + safety cap: precision may only be lowered after enough successful
                // steps in a row and while under the decrease budget; otherwise hold the current precision.
                const bool decrease_allowed =
                    (num_successful_steps_since_precision_decrease_ >= steps_for_increase)
                    &&
                    (num_precision_decreases_ < Get<PrecConf>().max_num_precision_decreases);
                if (!decrease_allowed)
                    min_precision = max(min_precision, current_precision_);
 
                unsigned max_precision = Get<PrecConf>().maximum_precision;
 
                try {
                    MinimizeTrackingCost(next_precision_, next_stepsize_,
                                min_precision, min_stepsize,
                                max_precision, max_stepsize,
                                DigitsB<ComplexT>(), // RAW digits_B (=P0) for the stepsize relation, per B1 minimize_cost
                                N,
                                predictor_order_);
                } catch (std::runtime_error const&) {
                    return SuccessCode::FailedToSelectPrecisionAndStepsize;
                }
 
                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_ = real_mp(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 ComplexT>
            SuccessCode AMPCriterionError() const
            {
 
                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
 
 
                real_mp min_stepsize = MinStepSizeForPrecision(current_precision_, abs(current_time_ - endtime_));
                real_mp max_stepsize = current_stepsize_ * real_mp(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+real_mp(Get<Stepping>().step_size_fail_factor, current_precision_))/2, min_stepsize);
                }
                else
                {
                    unsigned digits_B = DigitsB<ComplexT>(); // RAW digits_B (=P0); the stepsize relation in MinimizeTrackingCost needs it raw
 
                    // rule-B precision floor CREDITS THE STEP (B1 AMP2_update: digits_B2 = ceil(P0 - eta_minSS/N)):
                    // floor = digits_B - (p+1)*(-log10(min_stepsize))/N, crediting the smallest allowed step.
                    // This keeps a small step (large size_proportion in the roundoff regime) from forcing
                    // precision up by itself; genuine ill-conditioning still escalates via D / rules A,C.
                    const unsigned N = Get<NewtonConfig>().max_num_newton_iterations;
                    int digits_B_credited = static_cast<int>(ceil( double(digits_B)
                                                - (predictor_order_ + 1.0) * double(-log10(min_stepsize)) / N ));
                    if (digits_B_credited < 0)
                        digits_B_credited = 0;
 
                    unsigned min_precision = max(min_next_precision,
                                                 static_cast<unsigned>(digits_B_credited),
                                                 DigitsC<ComplexT>(),
                                                 MinDigitsForStepsizeInterval(min_stepsize, max_stepsize, abs(current_time_ - endtime_)),
                                                 digits_final_
                                                 );
 
                    try {
                        MinimizeTrackingCost(next_precision_, next_stepsize_,
                                min_precision, min_stepsize,
                                Get<PrecConf>().maximum_precision, max_stepsize,
                                digits_B,
                                N,
                                predictor_order_);
                    } catch (std::runtime_error const&) {
                        return SuccessCode::FailedToSelectPrecisionAndStepsize;
                    }
 
                }
 
                UpdatePrecisionAndStepsize();
                return SuccessCode::Success;
            }
 
 
 
 

            template<typename ComplexT>
            NumErrorT B_RHS() const
            {
                return max(amp::CriterionBRHS(this->last_step_.norm_J,
                                              this->last_step_.norm_J_inverse,
                                              Get<NewtonConfig>().max_num_newton_iterations,
                                              this->TrackingTolerance(),
                                              this->last_step_.norm_delta_z,
                                              Get<PrecConf>()), NumErrorT(0));
            }
 
 
 

            template<typename ComplexT>
            unsigned DigitsB() const
            {
                unsigned d = unsigned(B_RHS<ComplexT>());
                return d;
            }
 
 
 
 

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

            template<typename ComplexT>
            unsigned DigitsC() const
            {
                return unsigned(C_RHS<ComplexT>());
            }
 
 
 
 
 
            // (MinRequiredPrecision_BCTol was removed: it returned the UNCREDITED rule-B floor
            //  max(DigitsB, DigitsC, tol, double), which is what spuriously escalated precision on
            //  small steps.  AdjustAMPStepSuccess now computes the step-credited, margin-gated floor
            //  directly, faithful to Bertini 1's AMP2_update.)
 
 
 
 

            //
            //  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_stepsize_increase_ = 0;
                num_successful_steps_since_precision_decrease_ = 0;
                NotifyObservers(FailedStep<EmitterType>(*this));
            }
 
 
 
            void OnInfiniteTruncation() const override
            {
                NotifyObservers(InfinitePathTruncation<EmitterType>(*this));
            }
 

            //
            //       Predict and Correct functions
            //
 
 
 
 

            template<typename ComplexT, typename RealT>
            SuccessCode Predict(Vec<ComplexT> & predicted_space,
                                const Vec<ComplexT>& current_space,
                                ComplexT const& current_time, ComplexT const& delta_t) const
            {
 
 
 
 
                static_assert(std::is_same< typename Eigen::NumTraits<RealT>::Real,
                                        typename Eigen::NumTraits<ComplexT>::Real>::value,
                                        "underlying complex type and the type for comparisons must match");
 
 
                return predictor_.Predict(predicted_space,
                                this->last_step_,
                                tracked_system_,
                                current_space, current_time,
                                delta_t,
                                num_steps_since_last_condition_number_computation_,
                                Get<Stepping>().frequency_of_CN_estimation,
                                this->TrackingTolerance(),
                                &Get<PrecConf>());
            }
 
 

            template<typename ComplexT, typename RealT>
            SuccessCode Correct(Vec<ComplexT> & corrected_space,
                                Vec<ComplexT> const& current_space,
                                ComplexT const& current_time) const
            {
                static_assert(std::is_same< typename Eigen::NumTraits<RealT>::Real,
                                        typename Eigen::NumTraits<ComplexT>::Real>::value,
                                        "underlying complex type and the type for comparisons must match");
 
 
 
 
                return corrector_.Correct(corrected_space,
                                    this->last_step_,
                                    tracked_system_,
                                    current_space,
                                    current_time,
                                    this->TrackingTolerance(),
                                    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<complex_dbl>(std::get<Vec<complex_dbl> >(temporary_space_),std::get<Vec<complex_dbl> >(current_space_), complex_dbl(current_time_));
                    if (code == SuccessCode::Success)
                        std::get<Vec<complex_dbl> >(current_space_) = std::get<Vec<complex_dbl> >(temporary_space_);
                }
                else
                {
                    code = RefineImpl<complex_mp>(std::get<Vec<complex_mp> >(temporary_space_),std::get<Vec<complex_mp> >(current_space_), current_time_);
                    if (code == SuccessCode::Success)
                        std::get<Vec<complex_mp> >(current_space_) = std::get<Vec<complex_mp> >(temporary_space_);
                }
                return code;
            }
 
 

            template <typename ComplexT>
            SuccessCode RefineImpl(Vec<ComplexT> & new_space,
                                Vec<ComplexT> const& start_point, ComplexT const& current_time) const
            {
                using RealT = typename Eigen::NumTraits<ComplexT>::Real;
 
                static_assert(std::is_same< typename Eigen::NumTraits<RealT>::Real,
                                        typename Eigen::NumTraits<ComplexT>::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->last_step_,
                                           tracked_system_,
                                           start_point,
                                           current_time,
                                           this->TrackingTolerance(),
                                           Get<NewtonConfig>().min_num_newton_iterations,
                                           Get<NewtonConfig>().max_num_newton_iterations,
                                           &Get<PrecConf>());
            }
 
 
 
 
 

            template <typename ComplexT>
            SuccessCode RefineImpl(Vec<ComplexT> & new_space,
                                Vec<ComplexT> const& start_point, ComplexT 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->last_step_,
                                        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<complex_mp>() && "precision sanity check failed.  some internal variable is not in correct precision");
                    #endif
                }
                else
                {
                    MultipleToMultiple(new_precision);
                    #ifndef BERTINI_DISABLE_ASSERTS
                    assert(PrecisionSanityCheck<complex_mp>() && "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<complex_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();
                SetThreadPrecision(DoublePrecision());
 
 
                std::get<Vec<complex_dbl> >(current_space_) = source_point;
            }

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

            void MultipleToDouble(Vec<complex_mp> 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());
 
 
                // copy the current space in.
                if (std::get<Vec<complex_dbl> >(current_space_).size()!=source_point.size())
                    std::get<Vec<complex_dbl> >(current_space_).resize(source_point.size());
 
                for (unsigned ii=0; ii<source_point.size(); ii++)
                    std::get<Vec<complex_dbl> >(current_space_)(ii) = complex_dbl(source_point(ii));
 
                // The same internals every other conversion re-precisions: the step size, the
                // time, its increment, the end time, and the predictor's and corrector's work
                // space.  Re-precisioning only the end time here left the rest at whatever
                // multiple precision the track came down from, so a track's double-precision
                // steps depended on the precision history before them (#378).
                AdjustInternalsPrecision(DoublePrecision());
            }

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

            void DoubleToMultiple(unsigned new_precision, Vec<complex_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<complex_mp>();
                #endif
            }
 
 

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

            void MultipleToMultiple(unsigned new_precision, Vec<complex_mp> 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<complex_mp>();
                #endif
            }
 

            void MultipleToMultiple(unsigned new_precision) const
            {
                MultipleToMultiple( new_precision, std::get<Vec<complex_mp> >(current_space_));
            }
 
 
 
 
 
            void AdjustCurrentPrecision(unsigned new_precision) const
            {
                previous_precision_ = current_precision_;
                current_precision_ = new_precision;
                SetThreadPrecision(new_precision);
            }
 
 
            void CopyToCurrentSpace(Vec<complex_dbl> const& source_point) const
            {
                auto& space = std::get<Vec<complex_mp> >(current_space_);
                if (space.size()!=source_point.size())
                    space.resize(source_point.size());
                for (unsigned ii=0; ii<source_point.size(); ii++)
                    space(ii) = complex_mp(source_point(ii));
            }
 
            void CopyToCurrentSpace(Vec<complex_mp> const& source_point) const
            {
                auto& space = std::get<Vec<complex_mp> >(current_space_);
                if (space.size()!=source_point.size())
                    space.resize(source_point.size());
                space = source_point;
            }
 
            void AdjustInternalsPrecision(unsigned new_precision) const
            {
                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<complex_mp> >(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<complex_mp> >(tentative_space_).resize(static_cast<Eigen::Index>(num_vars));
                Precision(std::get<Vec<complex_mp> >(tentative_space_), new_precision);
 
                std::get<Vec<complex_mp> >(temporary_space_).resize(static_cast<Eigen::Index>(num_vars));
                Precision(std::get<Vec<complex_mp> >(temporary_space_), new_precision);
            }
 
 
 

            template<typename ComplexT>
            bool PrecisionSanityCheck() const
            {
 
                if constexpr (std::is_same<ComplexT, complex_dbl>::value){
                    return true;
                }
 
                if constexpr (std::is_same<ComplexT, complex_mp>::value){
                    assert(ThreadPrecision()==current_precision_ && "current precision differs from the thread-local default precision");
                    assert(GetSystem().precision() == current_precision_ && "tracked system is out of precision");
 
                    assert(std::get<Vec<complex_mp> >(current_space_)(0).precision() == current_precision_ && "current space out of precision");
                    assert(std::get<Vec<complex_mp> >(tentative_space_)(0).precision() == current_precision_ && "tentative space out of precision");
                    assert(std::get<Vec<complex_mp> >(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; 
            std::optional<unsigned> override_start_precision_; 
 
            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 complex_mp endtime_highest_precision_;  
 
        public:

            unsigned CurrentPrecision() const override
            {
                return current_precision_;
            }
        }; // re: class Tracker
 
    } // namespace tracking
} // namespace bertini
 
 
#endif