base_tracker.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
//base_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.
//
//base_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 base_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_BASE_TRACKER_HPP
#define BERTINI_BASE_TRACKER_HPP
 
#include <algorithm>
#include <chrono>
#include <optional>
//#include "bertini2/tracking/step.hpp"
#include "bertini2/trackers/ode_predictors.hpp"
#include "bertini2/trackers/newton_corrector.hpp"
#include "bertini2/logging.hpp"
 
#include "bertini2/detail/observable.hpp"
 
// Must be at the end of the include list
#include "bertini2/trackers/events.hpp"
 
#include "bertini2/detail/is_template_parameter.hpp"
#include "bertini2/detail/configured.hpp"
#include "bertini2/common/stop_request.hpp"
 
namespace bertini{
 
    namespace tracking{

        inline
        bool StepFailureIsTerminal(SuccessCode c)
        {
            return c == SuccessCode::FailedToSelectPrecisionAndStepsize;
        }
 

        template<class D>
        class Tracker :
            public Observable,
            public detail::Configured<
                typename TrackerTraits< D >::NeededConfigs
                    >
        {
            using NeededTypes = typename TrackerTraits< D >::NeededTypes;
            using BaseComplexT = typename TrackerTraits<D>::BaseComplexT;
            using BaseRealT = typename TrackerTraits<D>::BaseRealT;
 
            using ComplexT = BaseComplexT;
            using RealT = BaseRealT;
 
 
        public:
            using Config = detail::Configured< typename TrackerTraits< D >::NeededConfigs >;  
            using Stepping = SteppingConfig;  
            using Newton = NewtonConfig;  
            using PrecConf = typename TrackerTraits< D >::PrecisionConfig;  

            Tracker(System const& sys) : tracked_system_(std::ref(sys)),
                predictor_(predict::DefaultPredictor(), sys),
                corrector_(sys)
            {
                SetPredictor(predict::DefaultPredictor());
            }
 
 
 

            void Setup(Predictor new_predictor_choice,
                       NumErrorT const& tracking_tolerance,
                        SteppingConfig const& stepping,
                        NewtonConfig const& newton)
            {
                SetPredictor(new_predictor_choice);
                corrector_.Settings(newton);
 
                SetTrackingTolerance(tracking_tolerance);
 
                this->template Set<SteppingConfig>(stepping);
                this->template Set<NewtonConfig>(newton);
 
                current_stepsize_ = BaseRealT(stepping.initial_step_size);
            }
 

            void Setup(Predictor new_predictor_choice,
                       NumErrorT const& tracking_tolerance,
                        NumErrorT const& path_truncation_threshold,
                        SteppingConfig const& stepping,
                        NewtonConfig const& newton)
            {
                Setup(new_predictor_choice, tracking_tolerance, stepping, newton);
                SetInfiniteTruncationTolerance(path_truncation_threshold);
            }
 

            using Config::Get;
 
 

            void SetTrackingTolerance(NumErrorT const& tracking_tolerance)
            {
                if (tracking_tolerance <= 0)
                    throw std::runtime_error("tracking tolerance must be strictly positive");
 
                auto settings = this->template Get<TrackerConfig>();
                settings.tracking_tolerance = tracking_tolerance;
                this->Set(settings);
            }
 

            void SetInfiniteTruncationTolerance(NumErrorT const& tol)
            {
                if (tol <= 0)
                    throw std::runtime_error("truncation threshold must be strictly positive");
 
                auto settings = this->template Get<TrackerConfig>();
                settings.path_truncation_threshold = tol;
                this->Set(settings);
            }
 

            void RefreshConditionDirection()
            {
                corrector_.RefreshRandomDirection();
                // Share the SAME per-path probe with the predictor, so the predictor's and
                // corrector's ||J^{-1}|| estimates use one consistent direction (ADR-0024).
                predictor_.SetConditionProbe(corrector_.ConditionProbe());
            }
 

            SuccessCode TrackPath(Vec<ComplexT> & solution_at_endtime,
                                    ComplexT const& start_time, ComplexT const& endtime,
                                    Vec<ComplexT> const& start_point
                                    ) const
            {
                if (start_point.size()!=static_cast<Eigen::Index>(GetSystem().NumVariables()))
                    throw std::runtime_error("start point size must match the number of variables in the system to be tracked");
 
                // whatever route the settings were set by, they take effect here
                SyncSettings();
 
                SuccessCode initialization_code = TrackerLoopInitialization(start_time, endtime, start_point);
                if (initialization_code!=SuccessCode::Success)
                {
                    PostTrackCleanup();
                    return initialization_code;
                }
 
                NotifyObservers(TrackingStarted<typename TrackerTraits<D>::EventEmitterType>(static_cast<typename TrackerTraits<D>::EventEmitterType const&>(*this)));
 
                // as precondition to this while loop, the correct container, either complex_dbl or mpfr, must have the correct data.
                while (!IsSymmRelDiffSmall(current_time_,endtime_, Eigen::NumTraits<ComplexT>::epsilon()))
                {
                    // Has somebody outside asked us to stop?  Checked FIRST, before the
                    // tracker's own budgets, because it is the only one whose answer is not
                    // about this path: a keyboard interrupt is a statement about the process.
                    // Between steps is the only place we can honour it without threads, so a
                    // stop overruns by at most the step in flight -- which at high precision
                    // is not free, and is the price of not preempting mid-step.
                    if (StopRequested())
                    {
                        PostTrackCleanup();
                        return SuccessCode::ExternallyTerminated;
                    }
 
                    // Has our wall-clock deadline passed?  A deadline rather than a duration,
                    // because it must survive across TrackPath calls: an endgame issues hundreds
                    // of them for one path, and a per-call budget would restart with each.
                    // Same overrun as the stop above: at most the step in flight.
                    if (wall_clock_deadline_ && std::chrono::steady_clock::now() >= *wall_clock_deadline_)
                    {
                        PostTrackCleanup();
                        return SuccessCode::WallClockLimitReached;
                    }
 
                    SuccessCode pre_iteration_code = PreIterationCheck();
                    if (pre_iteration_code!=SuccessCode::Success)
                    {
                        PostTrackCleanup();
                        return pre_iteration_code;
                    }
 
                    using std::abs;
                    // compute the next delta_t
                    if (abs(endtime_-current_time_) < abs(current_stepsize_))
                        delta_t_ = endtime_-current_time_;
                    else
                        delta_t_ = current_stepsize_ * (endtime_ - current_time_)/abs(endtime_ - current_time_);
 
                    step_success_code_ = TrackerIteration();
 
                    if (infinite_path_truncation_ && (CheckGoingToInfinity()==SuccessCode::GoingToInfinity))
                    {
                        OnInfiniteTruncation();
                        PostTrackCleanup();
                        return SuccessCode::GoingToInfinity;
                    }
                    else if (step_success_code_==SuccessCode::Success)
                        OnStepSuccess();
                    else
                    {
                        OnStepFail();
                        // A terminal failure adjusted nothing, so the next iteration would be
                        // bit-identical to the one that just failed -- an infinite loop that
                        // burns a core forever instead of failing the path.
                        if (StepFailureIsTerminal(step_success_code_))
                        {
                            PostTrackCleanup();
                            return step_success_code_;
                        }
                    }
 
                }// re: while
 
 
                CopySolution(solution_at_endtime);
                PostTrackCleanup();
                return SuccessCode::Success;
            }
 
 
 

            template<typename C>
            SuccessCode Refine(Vec<C> & new_space,
                                Vec<C> const& start_point, C const& current_time) const
            {
 
                static_assert(detail::IsTemplateParameter<C,NeededTypes>::value,"complex type for refinement must be a used type for the tracker");
                return this->AsDerived().RefineImpl(new_space, start_point, current_time);
            }
 
 
 

            template<typename C>
            SuccessCode Refine(Vec<C> & new_space,
                                Vec<C> const& start_point, C const& current_time, NumErrorT const& tolerance, unsigned max_iterations) const
            {
                static_assert(detail::IsTemplateParameter<C,NeededTypes>::value,"complex type for refinement must be a used type for the tracker");
 
                return this->AsDerived().RefineImpl(new_space, start_point, current_time, tolerance, max_iterations);
            }
 

            void SetPredictor(Predictor new_predictor_choice)
            {
                auto settings = this->template Get<TrackerConfig>();
                settings.predictor = new_predictor_choice;
                this->Set(settings);
                SyncSettings();
            }
 

            Predictor GetPredictor() const
            {
                return this->template Get<TrackerConfig>().predictor;
            }
 

            void SyncSettings() const
            {
                auto const wanted = this->template Get<TrackerConfig>().predictor;
                if (predictor_.PredictorMethod() == wanted)
                    return;
                predictor_.PredictorMethod(wanted);
                predictor_order_ = predictor_.Order();
            }
 

            void SetSystem(const System & new_sys)
            {
                tracked_system_ = std::ref(new_sys);
                predictor_.ChangeSystem(tracked_system_);
                corrector_.ChangeSystem(tracked_system_);
            }

            const System& GetSystem() const
            {
                return tracked_system_.get();
            }

            unsigned NumTotalStepsTaken () const
            {
                return num_failed_steps_taken_ + num_successful_steps_taken_;
            }

            void ResetCumulativeStepCounts() const
            {
                num_successful_steps_cumulative_ = 0;
                num_failed_steps_cumulative_ = 0;
            }

            unsigned CumulativeSuccessfulSteps() const
            {
                return num_successful_steps_cumulative_;
            }

            unsigned CumulativeFailedSteps() const
            {
                return num_failed_steps_cumulative_;
            }

            void SetMaxWallClockTime(std::chrono::steady_clock::time_point deadline) const
            {
                wall_clock_deadline_ = deadline;
            }

            template <typename Rep, typename Period>
            void SetMaxWallClockDuration(std::chrono::duration<Rep, Period> duration) const
            {
                using clock = std::chrono::steady_clock;
                wall_clock_deadline_ = clock::now() + std::chrono::duration_cast<clock::duration>(duration);
            }

            void ClearMaxWallClockTime() const
            {
                wall_clock_deadline_.reset();
            }

            std::optional<std::chrono::steady_clock::time_point> MaxWallClockTime() const
            {
                return wall_clock_deadline_;
            }

            void SetStepSize(RealT const& new_stepsize) const
            {
                current_stepsize_ = new_stepsize;
            }
 

            void ReinitializeInitialStepSize(bool should_reinitialize_stepsize)
            {
                reinitialize_stepsize_ = should_reinitialize_stepsize;
            }
 
            virtual ~Tracker() = default;

            auto TrackingTolerance() const
            {
                return this->template Get<TrackerConfig>().tracking_tolerance;
            }

            unsigned DigitsTrackingTolerance() const
            {
                return NumTraits<double>::TolToDigits(TrackingTolerance());
            }

            auto InfiniteTruncationTolerance() const
            {
                return this->template Get<TrackerConfig>().path_truncation_threshold;
            }
 
        private:
 
            // convert the base tracker into the derived type.
            const D& AsDerived() const
            {
                return static_cast<const D&>(*this);
            }

            virtual
            SuccessCode TrackerLoopInitialization(ComplexT const& start_time, ComplexT const& end_time, Vec<ComplexT> const& start_point) const = 0;
 

            virtual
            SuccessCode PreIterationCheck() const = 0;

            virtual
            SuccessCode TrackerIteration() const = 0;

            virtual
            void CopySolution(Vec<ComplexT> & solution_at_endtime) const = 0;
 
            // virtual
            // void CopySolution(Vec<complex_dbl> & solution_at_endtime) const = 0;
 
 
        protected:
 
 
 
 

            template <typename ComplexT>
            SuccessCode CheckGoingToInfinity() const
            {
                if (!GetSystem().IsFinite(std::get<Vec<ComplexT> >(current_space_), InfiniteTruncationTolerance()))
                    return SuccessCode::GoingToInfinity;
                else
                    return SuccessCode::Success;
            }
 
 

            virtual
            void PostTrackCleanup() const
            {}

            virtual
            void ResetCounters() const = 0;
 
 
 

            void ResetCountersBase() const
            {
                // reset a bunch of counters to 0.
                num_consecutive_successful_steps_ = 0;
                num_successful_steps_taken_ = 0;
                num_failed_steps_taken_ = 0;
                num_consecutive_failed_steps_ = 0;
            }
 

            void IncrementBaseCountersSuccess() const
            {
                num_successful_steps_taken_++;
                num_successful_steps_cumulative_++;
                num_consecutive_successful_steps_++;
                current_time_ += delta_t_;
                num_consecutive_failed_steps_ = 0;
            }

            virtual
            void OnStepSuccess() const = 0;
 

            void IncrementBaseCountersFail() const
            {
                num_consecutive_successful_steps_=0;
                num_failed_steps_taken_++;
                num_failed_steps_cumulative_++;
                num_consecutive_failed_steps_++;
            }
 

            virtual
            void OnStepFail() const = 0;

            virtual
            SuccessCode CheckGoingToInfinity() const = 0;

            virtual
            void OnInfiniteTruncation() const = 0;
 
 
            std::reference_wrapper<const System> tracked_system_; 
 
            bool infinite_path_truncation_ = true; 
            bool reinitialize_stepsize_ = true; 
 
            // tracking the numbers of things
            // zeroed here as well as in ResetCountersBase: that reset runs at the top of
            // TrackPath, so until a path has been tracked these would otherwise be read
            // uninitialized -- NumTotalStepsTaken() on a fresh tracker returned whatever
            // was in the memory
            mutable unsigned num_successful_steps_taken_ = 0;  
            mutable unsigned num_consecutive_successful_steps_ = 0; 
            mutable unsigned num_consecutive_failed_steps_ = 0; 
            mutable unsigned num_failed_steps_taken_ = 0; 
            mutable unsigned num_successful_steps_cumulative_ = 0; 
            mutable unsigned num_failed_steps_cumulative_ = 0; 
            mutable std::optional<std::chrono::steady_clock::time_point> wall_clock_deadline_; 
 
 
            // configuration for tracking
            //
            // predictor and corrector are held BY VALUE so that copying a tracker
            // deep-copies them.  (They were previously shared_ptr, which made every
            // tracker copy share one predictor/corrector with its source — unusable
            // from multiple threads.)  Both hold only work buffers and settings; they
            // store no reference to the System, so plain memberwise copy is correct.
            // They are mutable for the same reason as the state members above: they
            // hold scratch space mutated during the logically-const TrackPath.
            mutable predict::ExplicitRKPredictor predictor_; 
            mutable unsigned predictor_order_; 
 
            mutable correct::NewtonCorrector corrector_;  
 
 
 
            unsigned digits_final_ = 0; 
 
            mutable ComplexT endtime_; 
            mutable ComplexT current_time_; 
            mutable ComplexT delta_t_; 
            mutable RealT current_stepsize_; 
 
 
            // permanent temporaries
            mutable RealT next_stepsize_; 
            mutable SuccessCode step_success_code_; 
 
 
 
            mutable unsigned num_steps_since_last_condition_number_computation_; 
            mutable unsigned num_successful_steps_since_stepsize_increase_; 
 
            using TupOfVec = typename NeededTypes::ToTupleOfVec;  
            using TupOfReal = typename NeededTypes::ToTupleOfReal;  
 
            mutable TupOfVec current_space_; 
            mutable TupOfVec tentative_space_; 
            mutable TupOfVec temporary_space_; 
 

            mutable StepMetadata last_step_;
 
 
 
 
            public:
 

            NumErrorT LatestConditionNumber() const
            {
                return this->last_step_.condition_number_estimate;
            }
 

            NumErrorT LatestErrorEstimate() const
            {
                return this->last_step_.error_estimate;
            }
 

            NumErrorT LatestNormOfStep() const
            {
                return this->last_step_.norm_delta_z;
            }

            void SetInfiniteTruncation(bool b)
            {
                infinite_path_truncation_ = b;
            }

            auto InfiniteTruncation()
            {
                return infinite_path_truncation_;
            }

            unsigned NumVariables() const
            {
                return static_cast<unsigned>(GetSystem().NumVariables());
            }

            auto CurrentTime() const
            {
                return current_time_;
            }

            auto DeltaT() const
            {
                return delta_t_;
            }

            auto CurrentStepsize() const
            {
                return current_stepsize_;
            }
 

            virtual Vec<ComplexT> CurrentPoint() const = 0;
 

            virtual unsigned CurrentPrecision() const = 0;
        };
 
 
 
    } // re: namespace tracking
} // re: namespace bertini
 
 
#endif