powerseries.hpp

//This file is part of Bertini 2.
//
//powerseries_endgame.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.
//
//powerseries_endgame.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 powerseries_endgame.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
// Tim Hodges, Colorado State University
 
 
#pragma once
 
#include "bertini2/endgames/base_endgame.hpp"
 
 
namespace bertini{ namespace endgame{
 

 
template<typename PrecT> 
class PowerSeriesEndgame : 
    public virtual EndgameBase<PowerSeriesEndgame<PrecT>, PrecT>
{
public:
    using BaseEGT = EndgameBase<PowerSeriesEndgame<PrecT>, PrecT>;
    using FinalEGT = PowerSeriesEndgame<PrecT>;
    using TrackerType = typename BaseEGT::TrackerType;
 
    using BaseComplexType = typename BaseEGT::BaseComplexType;
    using BaseRealType = typename BaseEGT::BaseRealType;
 
    using EmitterType = PowerSeriesEndgame<PrecT>;
 
protected:
    
    using EndgameBase<PowerSeriesEndgame<PrecT>, PrecT>::NotifyObservers;
 
    using TupleOfTimes = typename BaseEGT::TupleOfTimes;
    using TupleOfSamps = typename BaseEGT::TupleOfSamps;
 
    using BCT = BaseComplexType;
    using BRT = BaseRealType;
 
    using Configs = typename AlgoTraits<FinalEGT>::NeededConfigs;
    using ConfigsAsTuple = typename Configs::ToTuple;

    mutable unsigned upper_bound_on_cycle_number_;
  
    mutable TupleOfTimes times_;
          
    mutable TupleOfSamps samples_;
          
    mutable TupleOfSamps derivatives_;

    mutable Vec<BCT> rand_vector_;
 
    template<typename CT>
    void AssertSizesTimeSpace() const
    {
        const auto num_sample_points = this->EndgameSettings().num_sample_points;
        assert(std::get<SampCont<CT> >(samples_).size()==std::get<TimeCont<CT> >(times_).size() && "must have same number of samples in times and spaces");
        assert(std::get<SampCont<CT> >(samples_).size()>=num_sample_points && "must have sufficient number of samples");
    }
 
    template<typename CT>
    void AssertSizesTimeSpaceDeriv() const
    {
        const auto num_sample_points = this->EndgameSettings().num_sample_points;
        assert(std::get<SampCont<CT> >(samples_).size()==std::get<TimeCont<CT> >(times_).size() && "must have same number of samples in times and spaces");
        assert(std::get<SampCont<CT> >(samples_).size()==std::get<SampCont<CT> >(derivatives_).size() && "must have same number of samples in derivatives and spaces");
        assert(std::get<SampCont<CT> >(samples_).size()>=num_sample_points && "must have sufficient number of samples");
    }
 
public:
 
    auto UpperBoundOnCycleNumber() const { return upper_bound_on_cycle_number_;}
 
  
    template<typename CT>
    void ClearTimesAndSamples()
    {
        std::get<TimeCont<CT> >(times_).clear(); 
        std::get<SampCont<CT> >(samples_).clear();
    }
  
    template<typename CT>
    void SetTimes(TimeCont<CT> const& times_to_set) { std::get<TimeCont<CT> >(times_) = times_to_set;}
  
    template<typename CT>
    const auto& GetTimes() const {return std::get<TimeCont<CT> >(times_);}
 
 
    const BCT& LatestTimeImpl() const
    {
        return GetTimes<BCT>().back();
    }
  
    template<typename CT>
    void SetSamples(SampCont<CT> const& samples_to_set) { std::get<SampCont<CT> >(samples_) = samples_to_set;}
  
    template<typename CT>
    const auto& GetSamples() const {return std::get<SampCont<CT> >(samples_);}
  
    template<typename CT>
    void SetRandVec(int size) {rand_vector_ = Vec<CT>::Random(size);}
 
 
 
 
 
    explicit PowerSeriesEndgame(TrackerType const& tr, 
                                const ConfigsAsTuple& settings )
      : BaseEGT(tr, settings), EndgamePrecPolicyBase<TrackerType>(tr)
    {}
 
    template< typename... Ts >
    explicit
    PowerSeriesEndgame(TrackerType const& tr, const Ts&... ts ) : PowerSeriesEndgame(tr, Configs::Unpermute( ts... ) ) 
        {}
 
 

    template<typename CT>
    unsigned ComputeBoundOnCycleNumber()
    { 
        using RT = typename Eigen::NumTraits<CT>::Real;
        using std::log; using std::abs;
 
        const auto& samples = std::get<SampCont<CT> >(samples_);
        
        AssertSizesTimeSpace<CT>();
 
        auto num_samples = samples.size();
        const Vec<CT> & sample0 = samples[num_samples-3];
        const Vec<CT> & sample1 = samples[num_samples-2];
        const Vec<CT> & sample2 = samples[num_samples-1]; // most recent sample.  oldest samples at front of the container
 
 
// should this only be if the system is homogenized?
        CT rand_sum1 = ((sample1 - sample0).transpose()*rand_vector_).sum();
        CT rand_sum2 = ((sample2 - sample1).transpose()*rand_vector_).sum();
 
        if ( abs(rand_sum1)==0 || abs(rand_sum2)==0) // avoid division by 0
        {
            upper_bound_on_cycle_number_ = 1;
            return upper_bound_on_cycle_number_;
        }
 
        RT estimate = log(static_cast<RT>(this->EndgameSettings().sample_factor))/log(abs(rand_sum2/rand_sum1));
 
        if (estimate < 1) // would be nan if sample points are same as each other
            upper_bound_on_cycle_number_ = 1;
        else
        {
            using std::max;
            auto upper_bound = unsigned(round(estimate)*this->template Get<PowerSeriesConfig>().cycle_number_amplification);
            upper_bound_on_cycle_number_ = max(upper_bound,this->template Get<PowerSeriesConfig>().max_cycle_number);
        }
 
        return upper_bound_on_cycle_number_;
    }//end ComputeBoundOnCycleNumber
 
 
 

 
    template<typename CT>
    unsigned ComputeCycleNumber(CT const& t0)
    {
        using RT = typename Eigen::NumTraits<CT>::Real;
 
        const auto& samples = std::get<SampCont<CT> >(samples_);
        const auto& times   = std::get<TimeCont<CT> >(times_);
        const auto& derivatives = std::get<SampCont<CT> >(derivatives_);
 
        AssertSizesTimeSpaceDeriv<CT>();
        
        const Vec<CT> &most_recent_sample = samples.back();  
        const CT& most_recent_time = times.back();
 
        //Compute upper bound for cycle number.
        ComputeBoundOnCycleNumber<CT>();
 
 
        unsigned num_pts;
        if (samples.size() > this->EndgameSettings().num_sample_points)
            num_pts = this->EndgameSettings().num_sample_points;
        else 
            num_pts = this->EndgameSettings().num_sample_points-1;
 
 
        auto min_found_difference = Eigen::NumTraits<RT>::highest();
 
        TimeCont<CT> s_times(num_pts);
        SampCont<CT> s_derivatives(num_pts);
 
        auto offset = samples.size() - num_pts - 1; // -1 here to shift away from the back of the container
        for(unsigned int candidate = 1; candidate <= upper_bound_on_cycle_number_; ++candidate)
        {           
            using std::pow;
 
            std::tie(s_times, s_derivatives) = TransformToSPlane(candidate, t0, num_pts, ContStart::Front);
            RT cand_power{1/static_cast<RT>(candidate)};
            RT curr_diff = (HermiteInterpolateAndSolve<CT>(
                                  pow((most_recent_time-t0)/(times[0]-t0),cand_power), // the target time
                                  num_pts,s_times,samples,s_derivatives, ContStart::Front) // the input data
                             - 
                             most_recent_sample).template lpNorm<Eigen::Infinity>();
 
            if (curr_diff < min_found_difference)
            {
                min_found_difference = curr_diff;
                this->cycle_number_ = candidate;
            }
 
        }// end cc loop over cycle number possibilities
 
        return this->cycle_number_;
    }//end ComputeCycleNumber
 
 

    template<typename CT>
    void ComputeAllDerivatives()
    {
        auto& samples = std::get<SampCont<CT> >(samples_);
        auto& times   = std::get<TimeCont<CT> >(times_);
        auto& derivatives = std::get<SampCont<CT> >(derivatives_);
 
        assert((samples.size() == times.size()) && "must have same number of times and samples");
 
        if (tracking::TrackerTraits<TrackerType>::IsAdaptivePrec) // known at compile time
        {
            auto max_precision = this->EnsureAtUniformPrecision(times, samples);
            this->GetSystem().precision(max_precision);
        }
 
        //Compute dx_dt for each sample.
        derivatives.clear(); derivatives.resize(samples.size());
        for(unsigned ii = 0; ii < samples.size(); ++ii)
        {   
            derivatives[ii] = -this->GetSystem().Jacobian(samples[ii],times[ii]).lu().solve(this->GetSystem().TimeDerivative(samples[ii],times[ii]));
        }
    }
 

    template <typename CT>
    std::tuple<TimeCont<CT>, SampCont<CT>> TransformToSPlane(int cycle_num, CT const& t0, unsigned num_pts, ContStart shift_from)
    {
        if (cycle_num==0)
            throw std::runtime_error("cannot transform to s plane with cycle number 0");
        AssertSizesTimeSpaceDeriv<CT>();
 
 
        using RT = typename Eigen::NumTraits<CT>::Real;
 
        const auto& times   = std::get<TimeCont<CT> >(times_);
        const auto& derivatives  = std::get<SampCont<CT> >(derivatives_);
 
        RT c = static_cast<RT>(cycle_num);
        RT one_over_c = 1/c;
 
        unsigned offset_t, offset_d;
        if (shift_from == ContStart::Back)
        {
            offset_t = times.size()-num_pts;
            offset_d = derivatives.size()-num_pts;
        }
        else
            offset_t = offset_d = 0;
 
 
        TimeCont<CT> s_times(num_pts);
        SampCont<CT> s_derivatives(num_pts);
 
        CT time_shift = times[offset_t] - t0;
 
        for(unsigned ii = 0; ii < num_pts; ++ii){
            s_times[ii] = pow((times[ii+offset_t]-t0)/time_shift, one_over_c); 
            s_derivatives[ii] = derivatives[ii+offset_d]*( c*pow(s_times[ii],cycle_num-1))*time_shift;
        }
 
        return std::make_tuple(s_times, s_derivatives);
    }
 
 

    template<typename CT>
    SuccessCode ComputeApproximationOfXAtT0(Vec<CT>& result, const CT & t0)
    {   
        const auto c = ComputeCycleNumber<CT>(t0);
 
        auto num_pts = this->EndgameSettings().num_sample_points;
 
        TimeCont<CT> s_times;
        SampCont<CT> s_derivatives;
 
        std::tie(s_times, s_derivatives) = TransformToSPlane(c, t0, num_pts, ContStart::Back);
        // the data was transformed to be on the interval [0 1] so we can hard-code the time-to-solve as 0 here.
 
        Precision(result, Precision(s_derivatives.back()));
        result = HermiteInterpolateAndSolve(CT(0), num_pts, s_times, std::get<SampCont<CT> >(samples_), s_derivatives, ContStart::Back);
        return SuccessCode::Success;
    }//end ComputeApproximationOfXAtT0
 
 

    template<typename CT>
    SuccessCode AdvanceTime(const CT & target_time)
    {
        using RT = typename Eigen::NumTraits<CT>::Real;
        
        auto& samples = std::get<SampCont<CT> >(samples_);
        auto& times   = std::get<TimeCont<CT> >(times_);
        auto& derivatives  = std::get<SampCont<CT> >(derivatives_);
 
        AssertSizesTimeSpaceDeriv<CT>();
 
        Vec<CT> next_sample;
        CT next_time = (times.back() + target_time) * static_cast<RT>(this->EndgameSettings().sample_factor); //setting up next time value using the midpoint formula, sample_factor will give us some 
 
        if (abs(next_time - target_time) < this->EndgameSettings().min_track_time) // generalized for target_time not equal to 0.
        {
            NotifyObservers(MinTrackTimeReached<EmitterType>(*this));
            return SuccessCode::MinTrackTimeReached;
        }
 
        SuccessCode tracking_success = this->GetTracker().TrackPath(next_sample,times.back(),next_time,samples.back());
            if (tracking_success != SuccessCode::Success)
                return tracking_success;
 
        NotifyObservers(InEGOperatingZone<EmitterType>(*this));
 
        this->EnsureAtPrecision(next_time,Precision(next_sample));
    
 
        times.push_back(next_time);
        samples.push_back(next_sample);
 
    
 
 
        auto refine_success = this->RefineSample(samples.back(), next_sample,  times.back(), 
                                        this->FinalTolerance() * this->EndgameSettings().sample_point_refinement_factor,
                                        this->EndgameSettings().max_num_newton_iterations);
        if (refine_success != SuccessCode::Success)
        {
            NotifyObservers(RefiningFailed<EmitterType>(*this));
            return refine_success;
        }
        
        this->EnsureAtPrecision(times.back(),Precision(samples.back()));
 
        NotifyObservers(SampleRefined<EmitterType>(*this));
 
        // we keep one more samplepoint than needed around, for estimating the cycle number
        if (times.size() > this->EndgameSettings().num_sample_points+1)
        {
            times.pop_front();
            samples.pop_front();
        }
 
        return SuccessCode::Success;
    }
 
      
    template<typename CT>
    SuccessCode RunImpl(const CT & start_time, const Vec<CT> & start_point, CT const& target_time)
    {
        if (start_point.size()!=this->GetSystem().NumVariables())
        {
            std::stringstream err_msg;
            err_msg << "number of variables in start point for PSEG, " << start_point.size() << ", must match the number of variables in the system, " << this->GetSystem().NumVariables();
            throw std::runtime_error(err_msg.str());
        }
 
        DefaultPrecision(Precision(start_point));
 
        using RT = typename Eigen::NumTraits<CT>::Real;
        //Set up for the endgame.
        ClearTimesAndSamples<CT>();
 
 
        // unpack some references for easy use
        auto& samples = std::get<SampCont<CT> >(samples_);
        auto& times   = std::get<TimeCont<CT> >(times_);
        auto& derivatives  = std::get<SampCont<CT> >(derivatives_);
        Vec<CT>& latest_approx = this->final_approximation_;
        Vec<CT>& prev_approx = this->previous_approximation_;
 
        // this is for estimating a ... norm?
        SetRandVec<CT>(start_point.size());       
        
        
        auto initial_sample_success = this->ComputeInitialSamples(start_time, target_time, start_point, times, samples);
 
        if (initial_sample_success!=SuccessCode::Success)
        {
            NotifyObservers(EndgameFailure<EmitterType>(*this));
            return initial_sample_success;
        }
 
        this->template RefineAllSamples<CT>(samples, times);
        ComputeAllDerivatives<CT>();
 
 
 
        auto extrapolation_code = ComputeApproximationOfXAtT0(prev_approx, target_time);
        latest_approx = prev_approx;
 
        if (extrapolation_code != SuccessCode::Success)
            return extrapolation_code;
 
        RT norm_of_dehom_of_latest_approx;
        RT norm_of_dehom_of_prev_approx;
        if (this->SecuritySettings().level <= 0)
            norm_of_dehom_of_prev_approx = this->GetSystem().DehomogenizePoint(prev_approx).template lpNorm<Eigen::Infinity>();
 
 
        NumErrorT& approx_error = this->approximate_error_;
        approx_error = 1;
 
        while (approx_error > this->FinalTolerance())
        {
            auto advance_code = AdvanceTime<CT>(target_time);
            if (advance_code!=SuccessCode::Success)
            {
                NotifyObservers(EndgameFailure<EmitterType>(*this));
                return advance_code;
            }
 
            // this code is what bertini1 does... it refines all samples, like, all the time.
            this->template RefineAllSamples<CT>(samples, times);
            ComputeAllDerivatives<CT>();
 
            extrapolation_code = ComputeApproximationOfXAtT0(latest_approx, target_time);
            if (extrapolation_code!=SuccessCode::Success)
            {
                NotifyObservers(EndgameFailure<EmitterType>(*this));
                return extrapolation_code;
            }
 
            approx_error = static_cast<NumErrorT>((latest_approx - prev_approx).template lpNorm<Eigen::Infinity>());
            NotifyObservers(ApproximatedRoot<EmitterType>(*this));
 
 
            if(this->SecuritySettings().level <= 0)
            {
                norm_of_dehom_of_latest_approx = this->GetSystem().DehomogenizePoint(latest_approx).template lpNorm<Eigen::Infinity>();
                if(norm_of_dehom_of_latest_approx > this->SecuritySettings().max_norm && norm_of_dehom_of_prev_approx > this->SecuritySettings().max_norm)
                {
                    NotifyObservers(SecurityMaxNormReached<EmitterType>(*this));
                    return SuccessCode::SecurityMaxNormReached;
                }
                norm_of_dehom_of_prev_approx = norm_of_dehom_of_latest_approx;
            }
 
            
 
 
            Precision(prev_approx, Precision(latest_approx));
            prev_approx = latest_approx;
        } //end while   
 
        NotifyObservers(Converged<EmitterType>(*this));
        return SuccessCode::Success;
 
    } //end PSEG
 
 
    virtual ~PowerSeriesEndgame() = default;
}; // end powerseries class
 
 
 
 
}} // re: namespaces