fixed_precision_tracker.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
//fixed_precision_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.
//
//fixed_precision_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 fixed_precision_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_FIXED_PRECISION_TRACKER_HPP
#define BERTINI_FIXED_PRECISION_TRACKER_HPP
 
#include "bertini2/trackers/base_tracker.hpp"
 
 
namespace bertini{
 
    namespace tracking{
 
        using std::max;
        using std::min;
        using std::pow;
 
        using bertini::max;
 
 
 
 
 
 

        template<class DerivedT>
        class FixedPrecisionTracker : public Tracker<FixedPrecisionTracker<DerivedT>>
        {
        public:
 
            using BaseComplexT = typename TrackerTraits<DerivedT>::BaseComplexT;  
            using BaseRealT = typename TrackerTraits<DerivedT>::BaseRealT;  
 
            using ComplexT = BaseComplexT;  
            using RealT = BaseRealT;  
 
            virtual ~FixedPrecisionTracker() = default;
 
            using EmitterType = FixedPrecisionTracker<DerivedT>;  
            using Base = Tracker<FixedPrecisionTracker<DerivedT>>;  
 
            using Config =  typename Base::Config;  
            FORWARD_GET_CONFIGURED
            using Stepping = typename Base::Stepping;  
            using Newton = typename Base::Newton;  

            FixedPrecisionTracker(System const& sys) : Base(sys){}

            void PrecisionSetup(FixedPrecisionConfig const&)
            { }
 

            Vec<ComplexT> CurrentPoint() const override
            {
                return std::get<Vec<ComplexT>>(this->current_space_);
            }
 
 
 
            void ResetCounters() const override
            {
                Base::ResetCountersBase();
 
                this->num_successful_steps_since_stepsize_increase_ = 0;
                // initialize to the frequency so guaranteed to compute it the first try
                this->num_steps_since_last_condition_number_computation_ = this->Get<Stepping>().frequency_of_CN_estimation;
            }
 

            SuccessCode PreIterationCheck() const override
            {
                // The budget counts EVERY step, not only the successful ones.  A path that
                // fails steps has spent the effort whether or not it advanced, and counting
                // only successes leaves a failing path un-budgeted: it can fail forever at
                // 0% of its allowance.  See issue #410.
                if (this->NumTotalStepsTaken() >= Get<Stepping>().max_num_steps)
                    return SuccessCode::MaxNumStepsTaken;
                if (this->current_stepsize_ < Get<Stepping>().min_step_size)
                    return SuccessCode::MinStepSizeReached;
 
                return SuccessCode::Success;
            }
 
 
 
 
 
 
 
 
 
 
            void PostTrackCleanup() const override
            {
                this->NotifyObservers(TrackingEnded<EmitterType>(*this));
            }

            void CopySolution(Vec<ComplexT> & solution_at_endtime) const override
            {
 
                // the current precision is the precision of the output solution point.
 
                unsigned num_vars = static_cast<unsigned>(this->GetSystem().NumVariables());
                solution_at_endtime.resize(num_vars);
                for (unsigned ii=0; ii<num_vars; ii++)
                {
                    solution_at_endtime(ii) = std::get<Vec<ComplexT> >(this->current_space_)(ii);
                }
 
            }
 
 
 

            SuccessCode TrackerIteration() const override
            {
                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");
 
                this->NotifyObservers(NewStep<EmitterType >(*this));
 
                Vec<ComplexT>& predicted_space = std::get<Vec<ComplexT> >(this->temporary_space_); // this will be populated in the Predict step
                Vec<ComplexT>& current_space = std::get<Vec<ComplexT> >(this->current_space_); // the thing we ultimately wish to update
                ComplexT current_time = ComplexT(this->current_time_);
                ComplexT delta_t = ComplexT(this->delta_t_);
 
                SuccessCode predictor_code = Predict(predicted_space, current_space, current_time, delta_t);
 
                if (predictor_code!=SuccessCode::Success)
                {
                    this->NotifyObservers(FirstStepPredictorMatrixSolveFailure<EmitterType >(*this));
 
                    this->next_stepsize_ = RealT(Get<Stepping>().step_size_fail_factor)*this->current_stepsize_;
 
                    UpdateStepsize();
 
                    return predictor_code;
                }
 
                this->NotifyObservers(SuccessfulPredict<EmitterType , ComplexT>(*this, predicted_space));
 
                Vec<ComplexT>& tentative_next_space = std::get<Vec<ComplexT> >(this->tentative_space_); // this will be populated in the Correct step
 
                ComplexT tentative_next_time = current_time + delta_t;
 
                SuccessCode corrector_code = Correct(tentative_next_space,
                                                     predicted_space,
                                                     tentative_next_time);
 
                if (corrector_code == SuccessCode::GoingToInfinity)
                {
                    // there is no corrective action possible...
                    return corrector_code;
                }
                else if (corrector_code!=SuccessCode::Success)
                {
                    this->NotifyObservers(CorrectorMatrixSolveFailure<EmitterType >(*this));
 
                    this->next_stepsize_ = RealT(Get<Stepping>().step_size_fail_factor)*this->current_stepsize_;
                    UpdateStepsize();
 
                    return corrector_code;
                }
 
 
                this->NotifyObservers(SuccessfulCorrect<EmitterType , ComplexT>(*this, tentative_next_space));
 
                // copy the tentative vector into the current space vector;
                current_space = tentative_next_space;
                return SuccessCode::Success;
            }
 

            SuccessCode CheckGoingToInfinity() const override
            {
                return Base::template CheckGoingToInfinity<ComplexT>();
            }
 
 
 
 

            SuccessCode UpdateStepsize() const
            {
                this->SetStepSize(this->next_stepsize_);
                return SuccessCode::Success;
            }
 
 
 
 

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

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

            void OnStepFail() const override
            {
                Base::IncrementBaseCountersFail();
                this->num_successful_steps_since_stepsize_increase_ = 0;
                this->NotifyObservers(FailedStep<EmitterType >(*this));
            }
 
 

            void OnInfiniteTruncation() const override
            {
                this->NotifyObservers(InfinitePathTruncation<EmitterType>(*this));
            }

            //
            //
 

            SuccessCode Predict(Vec<ComplexT> & predicted_space,
                                Vec<ComplexT> const& current_space,
                                ComplexT const& current_time, ComplexT const& delta_t) const
            {
 
                return this->predictor_.Predict(
                            predicted_space,
                            this->last_step_,
                            this->tracked_system_,
                            current_space, current_time,
                            delta_t,
                            this->num_steps_since_last_condition_number_computation_,
                            Get<Stepping>().frequency_of_CN_estimation,
                            this->TrackingTolerance());
            }
 
 

            SuccessCode Correct(Vec<ComplexT> & corrected_space,
                                Vec<ComplexT> const& current_space,
                                ComplexT const& current_time) const
            {
                return this->corrector_.Correct(corrected_space,
                                                this->last_step_,
                                                this->tracked_system_,
                                                current_space,
                                                current_time,
                                                this->TrackingTolerance(),
                                                Get<Newton>().min_num_newton_iterations,
                                                Get<Newton>().max_num_newton_iterations);
            }
 
 

            SuccessCode RefineImpl(Vec<ComplexT> & new_space,
                                Vec<ComplexT> const& start_point, ComplexT const& current_time) const
            {
                return this->corrector_.Correct(new_space,
                               this->last_step_,
                               this->tracked_system_,
                               start_point,
                               current_time,
                               this->TrackingTolerance(),
                               Get<Newton>().min_num_newton_iterations,
                               Get<Newton>().max_num_newton_iterations);
            }
 
 
 
 
 

            SuccessCode RefineImpl(Vec<ComplexT> & new_space,
                                Vec<ComplexT> const& start_point, ComplexT const& current_time,
                                NumErrorT const& tolerance, unsigned max_iterations) const
            {
                return this->corrector_.Correct(new_space,
                               this->last_step_,
                               this->tracked_system_,
                               start_point,
                               current_time,
                               tolerance,
                               1,
                               max_iterations);
            }
 

            //
 
            // no additional state variables needed for the FixedPrecision base tracker types
 
        };
 
 
 

        class DoublePrecisionTracker : public FixedPrecisionTracker<DoublePrecisionTracker>
        {
        public:
            using BaseComplexT = complex_dbl;  
            using BaseRealT = double;  
 
            using EmitterType = typename TrackerTraits<DoublePrecisionTracker>::EventEmitterType;  
 

            DoublePrecisionTracker(class System const& sys) : FixedPrecisionTracker<DoublePrecisionTracker>(sys)
            {
                // Report the precision honestly: double tracking is always DoublePrecision() digits.
                FixedPrecisionConfig c = this->template Get<FixedPrecisionConfig>();
                c.precision = DoublePrecision();
                this->template Set<FixedPrecisionConfig>(c);
            }
 
 
            DoublePrecisionTracker() = delete;
 
            virtual ~DoublePrecisionTracker() = default;
 

            unsigned CurrentPrecision() const override
            {
                return DoublePrecision();
            }

            void PrecisionSetup(FixedPrecisionConfig const& c)
            {
                if (c.precision != 0 && c.precision != DoublePrecision())
                    throw std::runtime_error("double-precision tracking is fixed at "
                        + std::to_string(DoublePrecision())
                        + " digits; cannot set precision to " + std::to_string(c.precision)
                        + " (use mptype 'multiple' or 'adaptive')");
            }
 

            SuccessCode TrackerLoopInitialization(BaseComplexT const& start_time,
                                           BaseComplexT const& end_time,
                                           Vec<BaseComplexT> const& start_point) const override
            {
                this->NotifyObservers(Initializing<EmitterType,BaseComplexT>(*this,start_time, end_time, start_point));
 
                // set up the master current time and the current step size
                this->current_time_ = start_time;
                this->endtime_ = end_time;
                std::get<Vec<BaseComplexT> >(this->current_space_) = start_point;
                if (this->reinitialize_stepsize_)
                    this->SetStepSize(min(BaseRealT(Get<Stepping>().initial_step_size),abs(start_time-end_time)/Get<Stepping>().min_num_steps));
 
                ResetCounters();
 
                return SuccessCode::Success;
            }
 
 
        private:
 
        }; // re: DoublePrecisionTracker
 

        class MultiplePrecisionTracker : public FixedPrecisionTracker<MultiplePrecisionTracker>
        {
        public:
            using BaseComplexT = complex_mp;  
            using BaseRealT = real_mp;  
 
            using EmitterType = FixedPrecisionTracker<MultiplePrecisionTracker>;  
 

            MultiplePrecisionTracker(class System const& sys) : FixedPrecisionTracker<MultiplePrecisionTracker>(sys), precision_(DefaultPrecision())
            {
                SyncPrecisionToConfig();
            }
 
 
            MultiplePrecisionTracker() = delete;
 
            virtual ~MultiplePrecisionTracker() = default;
 

            unsigned CurrentPrecision() const override
            {
                return precision_;
            }

            void SetPrecision(unsigned p)
            {
                precision_ = p;
                SyncPrecisionToConfig();
            }

            void PrecisionSetup(FixedPrecisionConfig const& c)
            {
                if (c.precision != 0)
                    SetPrecision(c.precision);
            }
 

            SuccessCode TrackerLoopInitialization(BaseComplexT const& start_time,
                                           BaseComplexT const& end_time,
                                           Vec<BaseComplexT> const& start_point) const override
            {
 
                // ThreadPrecision: thread-local read, correct both on the main thread
                // (where it equals the global) and on std::thread workers.
                if (start_point(0).precision()!=ThreadPrecision())
                {
                    std::stringstream err_msg;
                    err_msg << "start point for fixed multiple precision tracker has differing precision from default (" << start_point(0).precision() << "!=" << ThreadPrecision() << "), tracking cannot start";
                    throw std::runtime_error(err_msg.str());
                }
 
                if (start_point(0).precision()!=this->CurrentPrecision())
                {
                    std::stringstream err_msg;
                    err_msg << "start point for fixed multiple precision tracker has differing precision from tracker's precision (" << start_point(0).precision() << "!=" << this->CurrentPrecision() << "), tracking cannot start";
                    throw std::runtime_error(err_msg.str());
                }
 
                if (ThreadPrecision()!=this->CurrentPrecision())
                {
                    std::stringstream err_msg;
                    err_msg << "current default precision differs from tracker's precision (" << ThreadPrecision() << "!=" << this->CurrentPrecision() << "), tracking cannot start";
                    throw std::runtime_error(err_msg.str());
                }
 
 
                this->NotifyObservers(Initializing<EmitterType,BaseComplexT>(*this,start_time, end_time, start_point));
 
                // Reset precision of all persistent members before assignment so that
                // BMP preserve_related_precision doesn't propagate stale high precision
                // from a previous TrackPath call into the new one.
                this->current_time_.precision(precision_);
                this->endtime_.precision(precision_);
                this->delta_t_.precision(precision_);
                this->current_stepsize_.precision(precision_);
                this->next_stepsize_.precision(precision_);
                Precision(std::get<Vec<BaseComplexT>>(this->current_space_), precision_);
                Precision(std::get<Vec<BaseComplexT>>(this->temporary_space_), precision_);
                Precision(std::get<Vec<BaseComplexT>>(this->tentative_space_), precision_);
                this->predictor_.ChangePrecision(precision_);
                this->corrector_.ChangePrecision(precision_);
 
                // set up the master current time and the current step size
                this->current_time_ = start_time;
                this->endtime_ = end_time;
                std::get<Vec<BaseComplexT> >(this->current_space_) = start_point;
                if (this->reinitialize_stepsize_)
                    this->SetStepSize(min(real_mp(Get<Stepping>().initial_step_size),real_mp(abs(start_time-end_time)/Get<Stepping>().min_num_steps)));
 
                ResetCounters();
 
                return SuccessCode::Success;
            }

            bool PrecisionSanityCheck() const
            {
                // the System is no longer checked here: it carries no precision of its own, and
                // evaluation materializes it from the point it is handed (ADR-0057)
                return ThreadPrecision()==precision_ &&
                        std::get<Vec<complex_mp> >(current_space_)(0).precision() == precision_ &&
                        std::get<Vec<complex_mp> >(tentative_space_)(0).precision() == precision_ &&
                        std::get<Vec<complex_mp> >(temporary_space_)(0).precision() == precision_ &&
                        Precision(this->endtime_)==precision_
                                ;
            }
        private:
 
            // Keep the stored FixedPrecisionConfig.precision equal to the tracker's actual precision,
            // so reading the config tells you the precision in effect (not the 0 sentinel).
            void SyncPrecisionToConfig()
            {
                FixedPrecisionConfig c = this->template Get<FixedPrecisionConfig>();
                c.precision = precision_;
                this->template Set<FixedPrecisionConfig>(c);
            }
 
            unsigned precision_;  
        }; // re: MultiplePrecisionTracker
    } // namespace tracking
} // namespace bertini
 
 
#endif