//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 . // // Copyright(C) 2015 - 2021 by Bertini2 Development Team // // See 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 /** \file amp_tracker.hpp \brief \brief Contains the AMPTracker type, and necessary functions. */ #ifndef BERTINI_AMP_TRACKER_HPP #define BERTINI_AMP_TRACKER_HPP #include #pragma once #include "bertini2/trackers/base_tracker.hpp" namespace bertini{ namespace tracking{ using std::max; using std::min; using std::pow; using bertini::max; /** The minimum time that can be represented effectively using a given precision \param precision The number of digits you want to use \param time_to_go The duration of remaining time to track \param safety_digits A buffer of extra digits to use \return The smallest permissible stepsize \todo This function has a hardcoded value which should be replaced */ 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; } /** \brief Just another name for MinTimeForCurrentPrecision \param precision The number of digits you want to use \param time_to_go The duration of remaining time to track \param safety_digits A buffer of extra digits to use \return The smallest permissible stepsize */ inline mpfr_float MinStepSizeForPrecision(unsigned precision, mpfr_float const& time_to_go, int safety_digits = 3) { return MinTimeForCurrentPrecision(precision, time_to_go, safety_digits); } /** \brief Compute the cost function for arithmetic versus precision. From \cite AMP2, \f$C(P)\f$. As currently implemented, this is \f$ 10.35 + 0.13 P \f$, where P is the precision. These numbers are stale, and need to be recomputed. Badly. Please do this. This function tells you the relative cost of arithmetic at a given precision. \todo Recompute this cost function for boost::multiprecision::mpfr_float \param precision An integral number of digits -- then this gives the cost for arithmetic \return A double indicating how expensive arithmetic at a given precision is. 1 is the base-line for double-precision. */ inline double ArithmeticCost(unsigned precision) { if (precision==DoublePrecision()) return 1; else return 10.35 + 0.13 * precision; } /** \brief Compute a stepsize satisfying AMP Criterion B with a given precision \param precision The current working precision \param digits_B The number of digits from Criterion B \param num_newton_iterations The number of remaining newton iterations allowed in the correction \param predictor_order The order of the predictor \return The stepsize \todo Remove the default value of the predictor order, as that seems weird to have */ 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)); } /** \brief The minimum number of digits based on a log10 of a stepsize. \param log_of_stepsize log10 of a stepsize \param time_to_go How much time you have left to track. \param safety_digits A buffer of extra digits to use, over the computed min. Yes, this is added in, in this function. \return An integral number of necessary digits to use */ 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(); } /** \brief For a given range of stepsizes under consideration, gives the minimum number of digits \param min_stepsize The smallest stepsize under consideration \param max_stepsize The largest stepsize under consideration \param time_to_go How much time is left to track. \return The min number of digits. This function assumes you are going to time=0, or that you have taken care of that difference. That is, time_to_go should be a duration, not the current time, unless you are tracking to t=0, in which case current time is the duration left to go. Dig it? */ 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)); } /** \brief Compute precision and stepsize minimizing the ArithmeticCost() of tracking. For a given range of precisions, an old stepsize, and a maximum stepsize, the Cost of tracking is computed, and a minimizer found. This function is used in the AMPTracker tracking loop, both in case of successful steps and in Criterion errors. \param[out] new_precision The minimizing precision. \param[out] new_stepsize The minimizing stepsize. \param[in] min_precision The minimum considered precision. \param[in] min_stepsize The minimum permitted stepsize. \param[in] max_precision The maximum considered precision. \param[in] max_stepsize The maximum permitted stepsize. \param[in] digits_B The number of digits required, according to CriterionB from \cite AMP1, \cite AMP2 \param[in] num_newton_iterations The number of allowed Newton corrector iterations. \param[in] predictor_order The order of the predictor being used. This is the order itself, not the order of the error estimate. \see ArithmeticCost */ template 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::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 \brief Functor-like class for tracking paths on a system ## Explanation The bertini::AMPTracker class enables tracking using Adaptive Multiple Precision on an arbitrary square homotopy. The intended usage is to: 1. Create a system, and instantiate some settings. 2. Create an AMPTracker, associating it to the system you are going to solve or track on. 3. Run AMPTracker::Setup and AMPTracker::PrecisionSetup, getting the settings in line for tracking. 4. Repeatedly, or as needed, call the AMPTracker::TrackPath function, feeding it a start point, and start and end times. The initial precision is that of the start point. Working precision and stepsize are adjusted automatically to get around nearby singularities to the path. If the endpoint is singular, this may very well fail, as prediction and correction get more and more difficult with proximity to singularities. The TrackPath method is intended to allow the user to track to nonsingular endpoints, or to an endgame boundary, from which an appropriate endgame will be called. ## Some notes The AMPTracker has internal state. That is, it stores the current state of the path being tracked as data members. After tracking has concluded, these statistics may be extracted. If you need additional accessors for these data, contact the software authors. The class additionally uses the Boost.Log library for logging. At time of this writing, the trivial logger is being used. As development continues we will move toward using a more sophisticated logging model. Suggestions are welcome. This class, like the other Tracker classes, uses mutable members to store the current state of the tracker. The tracking tolerances, settings, etc, remain constant throughout a track, but the internal state such as the current time or space values, will change. You can also expect the precision of the tracker to differ after a certain calls, too. ## Example Usage Below we demonstrate a basic usage of the AMPTracker class to track a single path. The pattern is as described above: create an instance of the class, feeding it the system to be tracked, and some configuration. Then, use the tracker to track paths of the system. \code{.cpp} DefaultPrecision(30); // set initial precision. This is not strictly necessary. using namespace bertini::tracking; // 1. Create the system Var x = Variable::Make("x"); Var y = Variable::Make("y"); Var t = Variable::Make("t"); System sys; VariableGroup v{x,y}; sys.AddFunction(pow(x,2) + (1-t)*x - 1); sys.AddFunction(pow(y,2) + (1-t)*x*y - 2); sys.AddPathVariable(t); sys.AddVariableGroup(v); auto AMP = bertini::tracking::AMPConfigFrom(sys); // 2. Create the Tracker object, associating the system to it. bertini::tracking::AMPTracker tracker(sys); SteppingConfig stepping_preferences; NewtonConfig newton_preferences; // 3. Get the settings into the tracker tracker.Setup(Predictor::Euler, mpfr_float("1e-5"), mpfr_float("1e5"), stepping_preferences, newton_preferences); tracker.PrecisionSetup(AMP); // 4. Create a start and end time. These are complex numbers. mpfr_complex t_start("1.0"); mpfr_complex t_end("0"); // 5. Create a start point, and container for the end point. Vec start_point(2); start_point << mpfr_complex("1"), mpfr_complex("1.414"); // set the value of the start point. This is Eigen syntax. Vec end_point; // 6. actually do the tracking SuccessCode tracking_success = tracker.TrackPath(end_point, t_start, t_end, start_point); // 7. and then onto whatever processing you are doing to the computed point. \endcode If this documentation is insufficient, please contact the authors with suggestions, or get involved! Pull requests welcomed. ## Testing * Test suite driving this class: AMP_tracker_basics. * File: test/tracking_basics/tracker_test.cpp * Functionality tested: Can use an AMPTracker to track in various situations, including tracking on nonsingular paths. Also track to a singularity on the square root function. Furthermore test that tracking fails to start from a singular start point, and tracking fails if a singularity is directly on the path being tracked. */ class AMPTracker : public Tracker { friend class Tracker; public: typedef Tracker Base; typedef typename TrackerTraits::EventEmitterType EmitterType; enum UpsampleRefinementOption { upsample_refine_off = 0, upsample_refine_on = 1 }; /** \brief Construct an Adaptive Precision tracker, associating to it a System. */ AMPTracker(class System const& sys) : Tracker(sys), current_precision_(DefaultPrecision()) { Set(AMPConfigFrom(sys)); } /** \brief Special additional setup call for the AMPTracker, selecting the config for adaptive precision. */ void PrecisionSetup(AdaptiveMultiplePrecisionConfig const& AMP_config) { Set(AMP_config); } const unsigned GetCurrentPrecision() const { return current_precision_; } /** \brief Switch preservation of precision after tracking on / off By default, precision is preserved after tracking, so the precision of the ambient workspace is returned to its previous state once Adaptive Precision tracking is done. */ void PrecisionPreservation(bool should_preseve_precision) { preserve_precision_ = should_preseve_precision; } virtual ~AMPTracker() = default; Vec CurrentPoint() const override { if (this->CurrentPrecision()==DoublePrecision()) { const auto& curr_vector = std::get>(this->current_space_); Vec returnme(NumVariables()); for (unsigned ii = 0; ii < NumVariables(); ++ii) { returnme(ii) = mpfr_complex(curr_vector(ii)); } return returnme; } else return std::get>(this->current_space_); } private: /** \brief Set up the internals of the tracker for a fresh start. Copies the start time, current stepsize, and start point. Adjusts the current precision to match the precision of the start point. Zeros counters. \param start_time The time at which to start tracking. \param end_time The time to which to track. \param start_point The space values from which to start tracking. */ SuccessCode TrackerLoopInitialization(mpfr_complex const& start_time, mpfr_complex const& end_time, Vec 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(*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().min_num_steps; SetStepSize(min(NumTraits::FromRational(Get().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(initial_precision_); ResetCounters(); ChangePrecision(start_point(0).precision()); auto initial_refinement_code = InitialRefinement(); #ifndef BERTINI_DISABLE_ASSERTS if (initial_precision_==DoublePrecision()){ PrecisionSanityCheck(); } else{ PrecisionSanityCheck(); } #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().frequency_of_CN_estimation; } /** \brief Run an initial refinement of the start point, to ensure in high enough precision to start. \return Whether initial refinement was successful. */ SuccessCode InitialRefinement() const { SuccessCode initial_refinement_code = RefineStoredPoint(); if (initial_refinement_code!=SuccessCode::Success) { do { if (current_precision_ > Get().maximum_precision) { NotifyObservers(SingularStartPoint(*this)); return SuccessCode::SingularStartPoint; } if (current_precision_==DoublePrecision()) initial_refinement_code = ChangePrecision(LowestMultiplePrecision()); else initial_refinement_code = ChangePrecision(current_precision_+PrecisionIncrement()); } while (initial_refinement_code!=SuccessCode::Success); } return SuccessCode::Success; } /** \brief Ensure that number of steps, stepsize, and precision still ok. \return Success if ok to keep going, and a different code otherwise. */ SuccessCode PreIterationCheck() const override { if (num_successful_steps_taken_ >= Get().max_num_steps) return SuccessCode::MaxNumStepsTaken; if (current_stepsize_ < Get().min_step_size) return SuccessCode::MinStepSizeReached; if (current_precision_ > Get().maximum_precision) return SuccessCode::MaxPrecisionReached; return SuccessCode::Success; } void PostTrackCleanup() const override { if (preserve_precision_) ChangePrecision(initial_precision_); NotifyObservers(TrackingEnded(*this)); } /** \brief Copy from the internally stored current solution into a final solution. If preservation of precision is on, this function first returns to the initial precision. \param[out] solution_at_endtime The solution at the end time */ void CopyFinalSolution(Vec & 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 >(current_space_)(ii)); } else { unsigned num_vars = GetSystem().NumVariables(); solution_at_endtime.resize(num_vars); for (unsigned ii=0; ii >(current_space_)(ii); solution_at_endtime(ii).precision(current_precision_); } } } /** \brief Run an iteration of the tracker loop. Predict and correct, adjusting precision and stepsize as necessary. \return Success if the step was successful, and a non-success code if something went wrong, such as a linear algebra failure or AMP Criterion violation. */ SuccessCode TrackerIteration() const override { if (current_precision_==DoublePrecision()) return TrackerIteration(); else return TrackerIteration(); } /** \brief Run an iteration of AMP tracking. \return SuccessCode indicating whether the iteration was successful. \tparam ComplexType The complex number type. \tparam RealType The real number type. */ template SuccessCode TrackerIteration() const // not an override, because it is templated { using RealType = typename Eigen::NumTraits::Real; #ifndef BERTINI_DISABLE_ASSERTS assert(PrecisionSanityCheck() && "precision sanity check failed. some internal variable is not in correct precision"); #endif NotifyObservers(NewStep(*this)); Vec& predicted_space = std::get >(temporary_space_); // this will be populated in the Predict step Vec& current_space = std::get >(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(); #endif SuccessCode predictor_code = Predict(predicted_space, current_space, current_time, delta_t); if (predictor_code==SuccessCode::MatrixSolveFailureFirstPartOfPrediction) { NotifyObservers(FirstStepPredictorMatrixSolveFailure(*this)); InitialMatrixSolveError(); return predictor_code; } else if (predictor_code==SuccessCode::MatrixSolveFailure) { NotifyObservers(PredictorMatrixSolveFailure(*this)); NewtonConvergenceError();// decrease stepsize, and adjust precision as necessary return predictor_code; } else if (predictor_code==SuccessCode::HigherPrecisionNecessary) { NotifyObservers(PredictorHigherPrecisionNecessary(*this)); AMPCriterionError(); return predictor_code; } NotifyObservers(SuccessfulPredict(*this, predicted_space)); Vec& tentative_next_space = std::get >(tentative_space_); // this will be populated in the Correct step ComplexType tentative_next_time = current_time + delta_t; SuccessCode corrector_code = Correct(tentative_next_space, predicted_space, tentative_next_time); #ifndef BERTINI_DISABLE_ASSERTS assert(PrecisionSanityCheck() && "precision sanity check failed. some internal variable is not in correct precision"); #endif if (corrector_code==SuccessCode::MatrixSolveFailure || corrector_code==SuccessCode::FailedToConverge) { NotifyObservers(CorrectorMatrixSolveFailure(*this)); NewtonConvergenceError(); return corrector_code; } else if (corrector_code == SuccessCode::HigherPrecisionNecessary) { NotifyObservers(CorrectorHigherPrecisionNecessary(*this)); AMPCriterionError(); return corrector_code; } else if (corrector_code == SuccessCode::GoingToInfinity) { // there is no corrective action possible... return corrector_code; } NotifyObservers(SuccessfulCorrect(*this, tentative_next_space)); // copy the tentative vector into the current space vector; current_space = tentative_next_space; return AdjustAMPStepSuccess(); } /** Check whether the path is going to infinity. */ SuccessCode CheckGoingToInfinity() const override { if (current_precision_ == DoublePrecision()) return Base::CheckGoingToInfinity(); else return Base::CheckGoingToInfinity(); } /** \brief Commit the next precision and stepsize, and adjust internals. \tparam refine_if_necessary Flag indicating whether to refine if precision if increasing. \see UpsampleRefinementOption */ template SuccessCode UpdatePrecisionAndStepsize() const { SetStepSize(next_stepsize_); return ChangePrecision(next_precision_); } /** \brief Increase stepsize or decrease precision, because of consecutive successful steps. \tparam ComplexType The complex number type. \tparam RealType The real number type. If the most recent step was successful, maybe adjust down precision and up stepsize. The number of consecutive successful steps is recorded as state in this class, and if this number exceeds a user-determine threshold, the precision or stepsize are allowed to favorably change. If not, then precision can only go up or remain the same. Stepsize can only decrease. These changes depend on the AMP criteria and current tracking tolerance. */ template 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::FromRational(Get().step_size_fail_factor, current_precision_); mpfr_float max_stepsize = min( current_stepsize_ * NumTraits::FromRational(Get().step_size_success_factor, current_precision_), NumTraits::FromRational(Get().max_step_size, current_precision_)); unsigned min_precision = MinRequiredPrecision_BCTol(); unsigned max_precision = max(min_precision,current_precision_); if (num_successful_steps_since_stepsize_increase_ < Get().consecutive_successful_steps_before_stepsize_increase) max_stepsize = current_stepsize_; // disallow stepsize changing if ( (num_successful_steps_since_precision_decrease_ < Get().consecutive_successful_steps_before_precision_decrease) || (num_precision_decreases_ >= Get().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(), Get().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(); } /** \brief The convergence_error function from \cite AMP2. Adjust precision and stepsize due to Newton's method failure to converge in the maximum number of steps. Decrease stepsize, then increase precision until the new stepsize is greater than a precision-dependent minimum stepsize. This function is used in the AMPTracker loop, when Newton's method fails to converge. \see MinStepSizeForPrecision */ void NewtonConvergenceError() const { next_precision_ = current_precision_; next_stepsize_ = mpfr_float(Get().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(); } /** \brief Adjust step size and precision due to AMP Criterion violation This function adjusts internals to the tracker object. Precision is REQUIRED to increase at least one increment. Stepsize is REQUIRED to decrease. \tparam ComplexType The complex number type. */ template void AMPCriterionError() const { using RealT = typename Eigen::NumTraits::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::FromRational(Get().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::FromRational(Get().step_size_fail_factor,current_precision_))/2, min_stepsize); } else { unsigned digits_B = DigitsB(); unsigned min_precision = max(min_next_precision, digits_B, DigitsC(), MinDigitsForStepsizeInterval(min_stepsize, max_stepsize, abs(current_time_ - endtime_)), digits_final_ ); unsigned a = ceil(digits_B - (predictor_order_+1)* -log10(max_stepsize)/Get().max_num_newton_iterations).convert_to(); unsigned max_precision = max(min_precision, a); MinimizeTrackingCost(next_precision_, next_stepsize_, min_precision, min_stepsize, Get().maximum_precision, max_stepsize, digits_B, Get().max_num_newton_iterations, predictor_order_); } UpdatePrecisionAndStepsize(); } /** \brief Get the raw right-hand side of Criterion B based on current state. */ template NumErrorT B_RHS() const { return max(amp::CriterionBRHS(this->norm_J_, this->norm_J_inverse_, Get().max_num_newton_iterations, tracking_tolerance_, this->size_proportion_, Get()), NumErrorT(0)); } /** \brief Get the right hand side of Criterion B based on current state. Returns the larger of the right hand side of B, or 1. Uses: * the most recent estimate on the norm of \f$J\f$, * the most recent estimate of the norm of \f$J^{-1}\f$, * the number of allowed newton iterations, * the tracking tolerance, * the size_proportion of the latest prediction, * the AMP configuration. */ template unsigned DigitsB() const { return unsigned(B_RHS()); } /** \brief Get the raw right-hand side of Criterion B based on current state. */ template NumErrorT C_RHS() const { return max(amp::CriterionCRHS(this->norm_J_inverse_, NumErrorT(std::get > (current_space_).norm()), tracking_tolerance_, Get()), NumErrorT(0)); } /** \brief Get the right hand side of Criterion C based on current state. Returns the larger of the right hand side of C, or 1. Uses: * the most recent estimate of the norm of \f$J^{-1}\f$, * the most recent current space point, * the current tracking tolerance, * the AMP configuration. */ template unsigned DigitsC() const { return unsigned(C_RHS()); } /** \brief Get the minimum required precision based on current state. The current state determining the minimum required precision is: * the norm of the Jacobian, * the norm of the inverse of the Jacobian, * the size_proportion, related to stepsize and predictor order, * the norm of the current solution. * the tracking tolerance. The min digits needed is the maximum of * Criterion B, * Criterion C, and * the digits required by the tracking tolerance. \tparam ComplexType The complex number type. \tparam RealType The real number type. */ template unsigned MinRequiredPrecision_BCTol() const { return max(DigitsB(), DigitsC(), digits_tracking_tolerance_, DoublePrecision()); } /////////// // // overrides for counter adjustment after a TrackerIteration() // //////////////// /** \brief Increment and reset counters after a successful TrackerIteration() */ void OnStepSuccess() const override { Tracker::IncrementBaseCountersSuccess(); NotifyObservers(SuccessfulStep(*this)); } /** \brief Increment and reset counters after a failed TrackerIteration() */ void OnStepFail() const override { Tracker::IncrementBaseCountersFail(); num_successful_steps_since_precision_decrease_ = 0; num_successful_steps_since_stepsize_increase_ = 0; NotifyObservers(FailedStep(*this)); } void OnInfiniteTruncation() const override { NotifyObservers(InfinitePathTruncation(*this)); } //////////////// // // Predict and Correct functions // ///////////////// /** \brief Wrapper function for calling the correct predictor. This function computes the next predicted space value, and sets some internals based on the prediction, such as the norm of the Jacobian. The real type and complex type must be commensurate. \param[out] predicted_space The result of the prediction \param current_space The current space point. \param current_time The current time value. \param delta_t The time differential for this step. Allowed to be complex. \tparam ComplexType The complex number type. \tparam RealType The real number type. */ template SuccessCode Predict(Vec & predicted_space, const Vec& current_space, ComplexType const& current_time, ComplexType const& delta_t) const { static_assert(std::is_same< typename Eigen::NumTraits::Real, typename Eigen::NumTraits::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().frequency_of_CN_estimation, tracking_tolerance_, Get()); 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().frequency_of_CN_estimation, tracking_tolerance_, Get()); } /** \brief Run Newton's method. Wrapper function for calling Correct and getting the error estimates etc directly into the tracker object. \tparam ComplexType The complex number type. \tparam RealType The real number type. \param corrected_space[out] The spatial result of the correction loop. \param current_space The start point in space for running the corrector loop. \param current_time The current time value. \return A SuccessCode indicating whether the loop was successful in converging in the max number of allowable newton steps, to the current path tolerance. */ template SuccessCode Correct(Vec & corrected_space, Vec const& current_space, ComplexType const& current_time) const { static_assert(std::is_same< typename Eigen::NumTraits::Real, typename Eigen::NumTraits::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().min_num_newton_iterations, Get().max_num_newton_iterations, Get()); } /** \brief Run Newton's method from the currently stored time and space value, in current precision. This overwrites the value of the current space, if successful, with the refined value. \return Whether the refinement was successful. */ SuccessCode RefineStoredPoint() const { SuccessCode code; if (current_precision_==DoublePrecision()) { code = RefineImpl(std::get >(temporary_space_),std::get >(current_space_), dbl(current_time_)); if (code == SuccessCode::Success) std::get >(current_space_) = std::get >(temporary_space_); } else { code = RefineImpl(std::get >(temporary_space_),std::get >(current_space_), current_time_); if (code == SuccessCode::Success) std::get >(current_space_) = std::get >(temporary_space_); } return code; } /** \brief Run Newton's method from a start point with a current time. Returns new space point by reference, as new_space. Operates at current precision. The tolerance is the tracking tolerance specified during Setup(...). \tparam ComplexType The complex number type. \tparam RealType The real number type. \param[out] new_space The result of running the refinement. \param start_point The base point for running Newton's method. \param current_time The current time value. \return Code indicating whether was successful or not. Regardless, the value of new_space is overwritten with the correction result. */ template SuccessCode RefineImpl(Vec & new_space, Vec const& start_point, ComplexType const& current_time) const { using RealType = typename Eigen::NumTraits::Real; static_assert(std::is_same< typename Eigen::NumTraits::Real, typename Eigen::NumTraits::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().min_num_newton_iterations, Get().max_num_newton_iterations, Get()); } /** \brief Run Newton's method from a start point with a current time. Returns new space point by reference, as new_space. Operates at current precision. \tparam ComplexType The complex number type. \tparam RealType The real number type. \param[out] new_space The result of running the refinement. \param start_point The base point for running Newton's method. \param current_time The current time value. \param tolerance The tolerance for convergence. This is a tolerance on \f$\Delta x\f$, not on function residuals. \param max_iterations The maximum allowable number of iterations to perform. \return Code indicating whether was successful or not. Regardless, the value of new_space is overwritten with the correction result. */ template SuccessCode RefineImpl(Vec & new_space, Vec 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()); } ///////////////// // // Functions for converting between precision types // /////////////////////// public: /** Change precision of tracker to next_precision. Converts the internal temporaries, and adjusts precision of system. Then refines if necessary. If the new precision is higher than current precision, a refine step will be called, which runs Newton's method. This may fail, leaving the tracker in a state with higher precision internals, but garbage digits after the previously known digits. \param new_precision The precision to change to. \return SuccessCode indicating whether the change was successful. If the precision increases, and the refinement loop fails, this could be not Success. Changing down is guaranteed to succeed. */ template SuccessCode ChangePrecision(unsigned new_precision) const { if (new_precision==current_precision_) // no op return SuccessCode::Success; NotifyObservers(PrecisionChanged(*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().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() && "precision sanity check failed. some internal variable is not in correct precision"); #endif } else { MultipleToMultiple(new_precision); #ifndef BERTINI_DISABLE_ASSERTS assert(PrecisionSanityCheck() && "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: /** \brief Converts from double to double Copies a multiple-precision into the double storage vector, and changes precision of the time and delta_t. \param source_point The point into which to copy to the internally stored current space point. */ void DoubleToDouble(Vec 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 >(current_space_) = source_point; } /** \brief Converts from multiple to double Changes the precision of the internal temporaries to double precision */ void DoubleToDouble() const { DoubleToDouble(std::get >(current_space_)); } /** \brief Converts from multiple to double Copies a multiple-precision into the double storage vector, and changes precision of the time and delta_t. \param source_point The point into which to copy to the internally stored current space point. */ void MultipleToDouble(Vec 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 >(current_space_).size()!=source_point.size()) std::get >(current_space_).resize(source_point.size()); for (unsigned ii=0; ii >(current_space_)(ii) = dbl(source_point(ii)); endtime_.precision(DoublePrecision()); // i question this one 2021-04-12 } /** \brief Converts from multiple to double Changes the precision of the internal temporaries to double precision */ void MultipleToDouble() const { MultipleToDouble(std::get >(current_space_)); } /** \brief Converts from double to multiple Copies a double-precision into the multiple-precision storage vector. You should call Refine after this to populate the new digits with non-garbage data. \param new_precision The new precision. \param source_point The point into which to copy to the internally stored current space point. */ void DoubleToMultiple(unsigned new_precision, Vec 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(); #endif } /** \brief Converts from double to multiple Copies the double-precision temporaries into the multiple-precision temporaries. You should call Newton after this to populate the new digits with non-garbage data. \param new_precision The new precision. */ void DoubleToMultiple(unsigned new_precision) const { DoubleToMultiple( new_precision, std::get >(current_space_)); } /** \brief Converts from multiple to different precision multiple precision Copies a multiple-precision into the multiple-precision storage vector, and changes precision of the time and delta_t. Also resets counter so have to re-compute the condition number on next step attempt. \param new_precision The new precision. \param source_point The point into which to copy to the internally stored current space point. */ void MultipleToMultiple(unsigned new_precision, Vec 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(); #endif } /** \brief Converts from multiple to different precision multiple precision Changes the precision of the internal temporaries to desired precision \param new_precision The new precision. */ void MultipleToMultiple(unsigned new_precision) const { MultipleToMultiple( new_precision, std::get >(current_space_)); } void AdjustCurrentPrecision(unsigned new_precision) const { previous_precision_ = current_precision_; current_precision_ = new_precision; DefaultPrecision(new_precision); } void CopyToCurrentSpace(Vec const& source_point) const { auto& space = std::get >(current_space_); if (space.size()!=source_point.size()) space.resize(source_point.size()); for (unsigned ii=0; ii const& source_point) const { auto& space = std::get >(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 >(current_space_),new_precision); } /** \brief Change precision of all temporary internal state variables. This excludes those which cannot be re-written without copying -- the current space point most notably. \brief new_precision The new precision to adjust to. */ 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 >(tentative_space_).resize(num_vars); Precision(std::get >(tentative_space_), new_precision); std::get >(temporary_space_).resize(num_vars); Precision(std::get >(temporary_space_), new_precision); } /** \brief Ensure that all internal state is in uniform precision. \return True if all internal state variables are in the same precision as current_precision_, will fail on assertion otherwise. */ template bool PrecisionSanityCheck() const { if constexpr (std::is_same::value){ return true; } if constexpr (std::is_same::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 >(current_space_)(0).precision() == current_precision_ && "current space out of precision"); assert(std::get >(tentative_space_)(0).precision() == current_precision_ && "tentative space out of precision"); assert(std::get >(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; } } ///////////////////////////////////////////// ////////////////////////////////////// ///////////////////////////// //////////////////// data members stored in this class //////////// ////// // //////////// // state variables ///////////// bool preserve_precision_ = false; ///< Whether the tracker should change back to the initial precision after tracking paths. mutable unsigned previous_precision_; ///< The previous precision of the tracker. mutable unsigned current_precision_; ///< The current precision of the tracker, the system, and all temporaries. mutable unsigned next_precision_; ///< The next precision mutable unsigned num_precision_decreases_; ///< The number of times precision has decreased this track. mutable unsigned initial_precision_; ///< The precision at the start of tracking. mutable unsigned num_successful_steps_since_precision_decrease_; ///< The number of successful steps since decreased precision. mutable mpfr_complex endtime_highest_precision_; public: unsigned CurrentPrecision() const override { return current_precision_; } }; // re: class Tracker } // namespace tracking } // namespace bertini #endif