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) 2015 - 2021 by Bertini2 Development Team
//
// See <http://www.gnu.org/licenses/> for a copy of the license, 
// as well as COPYING.  Bertini2 is provided with permitted 
// additional terms in the b2/licenses/ directory.
 
// individual authors of this file include:
// silviana amethyst, university of wisconsin eau claire

 
#ifndef BERTINI_BASE_TRACKER_HPP
#define BERTINI_BASE_TRACKER_HPP
 
#include <algorithm>
//#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"
 
namespace bertini{
 
    namespace tracking{

        template<class D>
        class Tracker : 
            public Observable,
            public detail::Configured<
                typename TrackerTraits< D >::NeededConfigs
                    >
        {
            using NeededTypes = typename TrackerTraits< D >::NeededTypes;
            using BaseComplexType = typename TrackerTraits<D>::BaseComplexType;
            using BaseRealType = typename TrackerTraits<D>::BaseRealType;
 
            using CT = BaseComplexType;
            using RT = BaseRealType;
            
            
        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_ = std::make_shared< predict::ExplicitRKPredictor >(predict::DefaultPredictor(), tracked_system_);
                corrector_ = std::make_shared< correct::NewtonCorrector >(tracked_system_);
                SetPredictor(predict::DefaultPredictor());
            }
 
 
 

            void Setup(Predictor new_predictor_choice,
                       double const& tracking_tolerance,
                        double const& path_truncation_threshold,
                        SteppingConfig const& stepping,
                        NewtonConfig const& newton)
            {
                SetPredictor(new_predictor_choice);
                corrector_->Settings(newton);
                
                SetTrackingTolerance(tracking_tolerance);
 
                path_truncation_threshold_ = path_truncation_threshold;
 
                this->template Set(stepping);
                this->template Set(newton);
 
                current_stepsize_ = BaseRealType(stepping.initial_step_size);
            }
 
 
            using Config::Get;
 
 

            void SetTrackingTolerance(double const& tracking_tolerance)
            {
                if (tracking_tolerance <= 0)
                    throw std::runtime_error("tracking tolerance must be strictly positive");
 
                tracking_tolerance_ = tracking_tolerance;
                digits_tracking_tolerance_ = NumTraits<double>::TolToDigits(tracking_tolerance);
            }
 
 
            void SetInfiniteTruncationTolerance(double const& tol)
            {
                if (tol <= 0)
                    throw std::runtime_error("truncation threshold must be strictly positive");
 
                path_truncation_threshold_ = tol;
            }
 

            SuccessCode TrackPath(Vec<CT> & solution_at_endtime,
                                    CT const& start_time, CT const& endtime,
                                    Vec<CT> const& start_point
                                    ) const
            {   
                if (start_point.size()!=GetSystem().NumVariables())
                    throw std::runtime_error("start point size must match the number of variables in the system to be tracked");
 
                
                
                SuccessCode initialization_code = TrackerLoopInitialization(start_time, endtime, start_point);
                if (initialization_code!=SuccessCode::Success)
                {
                    PostTrackCleanup();
                    return initialization_code;
                }
 
                // as precondition to this while loop, the correct container, either dbl or mpfr, must have the correct data.
                while (!IsSymmRelDiffSmall(current_time_,endtime_, Eigen::NumTraits<CT>::epsilon()))
                {   
                    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();
 
                }// re: while
 
 
                CopyFinalSolution(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, double 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)
            {
                predictor_->PredictorMethod(new_predictor_choice);
                predictor_order_ = predictor_->Order();
            }
 

            Predictor GetPredictor() const
            {
                return predictor_->PredictorMethod();
            }
 

            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 SetStepSize(RT 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 tracking_tolerance_;
            }
 
            auto InfiniteTruncationTolerance() const
            {
                return 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(CT const& start_time, CT const& end_time, Vec<CT> const& start_point) const = 0;
            

            virtual 
            SuccessCode PreIterationCheck() const = 0;

            virtual 
            SuccessCode TrackerIteration() const = 0;

            virtual
            void CopyFinalSolution(Vec<CT> & solution_at_endtime) const = 0;
 
            // virtual
            // void CopyFinalSolution(Vec<dbl> & solution_at_endtime) const = 0;
 
 
        protected:
 
 
            
 
 
            template <typename ComplexType>
            SuccessCode CheckGoingToInfinity() const
            {
                if (GetSystem().DehomogenizePoint(std::get<Vec<ComplexType> >(current_space_)).norm() > path_truncation_threshold_)
                    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;
                num_total_steps_taken_ = 0;
            }
 

            void IncrementBaseCountersSuccess() const
            {
                num_successful_steps_taken_++; 
                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_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
            mutable unsigned num_total_steps_taken_; 
            mutable unsigned num_successful_steps_taken_;  
            mutable unsigned num_consecutive_successful_steps_; 
            mutable unsigned num_consecutive_failed_steps_; 
            mutable unsigned num_failed_steps_taken_; 
 
            
            // configuration for tracking
            std::shared_ptr<predict::ExplicitRKPredictor > predictor_; // The predictor to use while tracking
            unsigned predictor_order_; 
 
            std::shared_ptr<correct::NewtonCorrector> corrector_;
 
 
 
            unsigned digits_final_ = 0; 
            unsigned digits_tracking_tolerance_ = 5; 
            NumErrorT tracking_tolerance_ = 1e-5; 
            NumErrorT path_truncation_threshold_ = 1e5; 
 
            mutable CT endtime_; 
            mutable CT current_time_; 
            mutable CT delta_t_; 
            mutable RT current_stepsize_; 
 
 
            // permanent temporaries
            mutable RT next_stepsize_; 
            mutable SuccessCode step_success_code_; 
 
 
 
            mutable unsigned num_steps_since_last_condition_number_computation_; 
            mutable unsigned num_successful_steps_since_stepsize_increase_; 
            mutable unsigned num_successful_steps_since_precision_decrease_; 
 
            using TupOfVec = typename NeededTypes::ToTupleOfVec;
            using TupOfReal = typename NeededTypes::ToTupleOfReal;
            
            mutable TupOfVec current_space_; 
            mutable TupOfVec tentative_space_; 
            mutable TupOfVec temporary_space_; 
 
 
            mutable NumErrorT condition_number_estimate_; 
            mutable NumErrorT error_estimate_; 
            mutable NumErrorT norm_J_; 
            mutable NumErrorT norm_J_inverse_;
            mutable NumErrorT norm_delta_z_; 
            mutable NumErrorT size_proportion_; 
 
 
 
 
            public: 
 
            
            NumErrorT LatestConditionNumber() const
            {
                return this->condition_number_estimate_;
            }
 
            
            NumErrorT LatestErrorEstimate() const
            {
                return this->error_estimate_;
            }
 
            
            NumErrorT LatestNormOfStep() const
            {
                return this->norm_delta_z_;
            }
 
            void SetInfiniteTruncation(bool b)
            {
                infinite_path_truncation_ = b;
            }
    
            auto InfiniteTruncation()
            {
                return infinite_path_truncation_;
            }
 
            unsigned NumVariables() const
            {
                return GetSystem().NumVariables();
            }
 
            auto CurrentTime() const
            {
                return current_time_;
            }
 
            auto DeltaT() const
            {
                return delta_t_;
            }
 
            auto CurrentStepsize() const
            {
                return current_stepsize_;
            }
 
 
            virtual Vec<CT> CurrentPoint() const = 0;
 
 
            virtual unsigned CurrentPrecision() const = 0;
        };
 
 
 
    } // re: namespace tracking
} // re: namespace bertini
 
 
#endif