powerseries.hpp

Go to the documentation of this file.

//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) 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.
 
#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 BaseComplexT = typename BaseEGT::BaseComplexT;  
    using BaseRealT = typename BaseEGT::BaseRealT;  
 
    using EmitterType = PowerSeriesEndgame<PrecT>;  
 
protected:
 
    using EndgameBase<PowerSeriesEndgame<PrecT>, PrecT>::NotifyObservers;
 
    using TupleOfTimes = typename BaseEGT::TupleOfTimes;  
    using TupleOfSamps = typename BaseEGT::TupleOfSamps;  
    using TupOfVec = typename BaseEGT::TupOfVec;  
 
    using BCT = BaseComplexT;  
    using BRT = BaseRealT;  
 
    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 TupleOfTimes c_over_k_;

    mutable bool announced_operating_zone_ = false;

    mutable TupOfVec rand_vector_;

    template<typename ComplexT>
    void AssertSizesTimeSpace() const
    {
#ifndef NDEBUG
        const auto num_sample_points = this->EndgameSettings().num_sample_points;
        assert(std::get<SampCont<ComplexT> >(samples_).size()==std::get<TimeCont<ComplexT> >(times_).size() && "must have same number of samples in times and spaces");
        assert(std::get<SampCont<ComplexT> >(samples_).size()>=num_sample_points && "must have sufficient number of samples");
#endif
    }

    template<typename ComplexT>
    void AssertSizesTimeSpaceDeriv() const
    {
#ifndef NDEBUG
        const auto num_sample_points = this->EndgameSettings().num_sample_points;
        assert(std::get<SampCont<ComplexT> >(samples_).size()==std::get<TimeCont<ComplexT> >(times_).size() && "must have same number of samples in times and spaces");
        assert(std::get<SampCont<ComplexT> >(samples_).size()==std::get<SampCont<ComplexT> >(derivatives_).size() && "must have same number of samples in derivatives and spaces");
        assert(std::get<SampCont<ComplexT> >(samples_).size()>=num_sample_points && "must have sufficient number of samples");
#endif
    }
 
public:

    auto UpperBoundOnCycleNumber() const { return upper_bound_on_cycle_number_;}
 

    template<typename ComplexT>
    void ClearTimesAndSamples()
    {
        std::get<TimeCont<ComplexT> >(times_).clear();
        std::get<SampCont<ComplexT> >(samples_).clear();
        // the operating-zone question is asked afresh for each path
        std::get<TimeCont<ComplexT> >(c_over_k_).clear();
        announced_operating_zone_ = false;
    }

    template<typename ComplexT>
    void SetTimes(TimeCont<ComplexT> const& times_to_set) { std::get<TimeCont<ComplexT> >(times_) = times_to_set;}

    template<typename ComplexT>
    const auto& GetTimes() const {return std::get<TimeCont<ComplexT> >(times_);}
 
 
    mutable BCT latest_time_cache_;  

    const BCT& LatestTimeImpl() const
    {
        // In the adaptive-numeric-type endgame the latest time may live in the complex_dbl slot, with
        // the BCT slot empty.  Dispatch on which slot holds data so this is correct in either lane and
        // both during and after the run.  Fixed precision compiles to the original BCT read.
        if constexpr (tracking::TrackerTraits<TrackerType>::IsAdaptivePrec)
        {
            if (GetTimes<BCT>().empty())
            {
                auto const& dbl_times = GetTimes<complex_dbl>();
                latest_time_cache_ = dbl_times.empty() ? BCT(0) : BCT(dbl_times.back());
                return latest_time_cache_;
            }
        }
        return GetTimes<BCT>().back();
    }

    template<typename ComplexT>
    void SetSamples(SampCont<ComplexT> const& samples_to_set) { std::get<SampCont<ComplexT> >(samples_) = samples_to_set;}

    template<typename ComplexT>
    const auto& GetSamples() const {return std::get<SampCont<ComplexT> >(samples_);}

    template<typename ComplexT>
    void SetRandVec(int size)
    {
        auto& rv = std::get<Vec<ComplexT> >(rand_vector_);
        rv.resize(size);
        for (int ii = 0; ii < size; ++ii)
            rv(ii) = RandomUnit<ComplexT>();
    }
 
 
 
 

    explicit PowerSeriesEndgame(TrackerType const& tr,
                                const ConfigsAsTuple& settings )
      : EndgamePrecPolicyBase<TrackerType>(tr), BaseEGT(tr, settings)
    {}

    template< typename... Ts >
    explicit
    PowerSeriesEndgame(TrackerType const& tr, const Ts&... ts ) : PowerSeriesEndgame(tr, Configs::Unpermute( ts... ) )
        {}
 
 

    template<typename ComplexT>
    unsigned ComputeBoundOnCycleNumber()
    {
        using RealT = typename Eigen::NumTraits<ComplexT>::Real;
        using std::log; using std::abs;
 
        const auto& samples = std::get<SampCont<ComplexT> >(samples_);
 
        AssertSizesTimeSpace<ComplexT>();
 
        auto num_samples = samples.size();
        const Vec<ComplexT> & sample0 = samples[num_samples-3];
        const Vec<ComplexT> & sample1 = samples[num_samples-2];
        const Vec<ComplexT> & sample2 = samples[num_samples-1]; // most recent sample.  oldest samples at front of the container
 
 
// should this only be if the system is homogenized?
        const auto& rand_vector = std::get<Vec<ComplexT> >(rand_vector_);
        ComplexT rand_sum1 = ((sample1 - sample0).transpose()*rand_vector).sum();
        ComplexT 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_;
        }
 
        RealT estimate = log(static_cast<RealT>(this->EndgameSettings().sample_factor))/log(abs(rand_sum2/rand_sum1));
 
        const auto& ps_config = this->template Get<PowerSeriesConfig>();
        if (estimate < 1)
            upper_bound_on_cycle_number_ = 1;
        else
        {
            // max_cycle_number is a CEILING on the candidate search (each candidate costs a full
            // Hermite solve).  Clamp before any conversion to unsigned: near-unity sample ratios
            // (slow convergence -- high multiplicity, or a slow diverger) drive the estimate toward
            // +inf, and unsigned(inf) is undefined behavior.  The !(a < b) form also routes NaN to
            // the ceiling.  This was `max` -- the ceiling was a floor, and the search was unbounded.
            using std::round;
            RealT amplified = round(estimate) * static_cast<RealT>(ps_config.cycle_number_amplification);
            if (!(amplified < static_cast<RealT>(ps_config.max_cycle_number)))
                upper_bound_on_cycle_number_ = ps_config.max_cycle_number;
            else
                upper_bound_on_cycle_number_ = std::max(1u, static_cast<unsigned>(amplified));
        }
 
        return upper_bound_on_cycle_number_;
    }//end ComputeBoundOnCycleNumber
 
 
 

 
    template<typename ComplexT>
    unsigned ComputeCycleNumber(ComplexT const& t0)
    {
        using RealT = typename Eigen::NumTraits<ComplexT>::Real;
 
        const auto& samples = std::get<SampCont<ComplexT> >(samples_);
        const auto& times   = std::get<TimeCont<ComplexT> >(times_);
 
        AssertSizesTimeSpaceDeriv<ComplexT>();
 
        const Vec<ComplexT> &most_recent_sample = samples.back();
        const ComplexT& most_recent_time = times.back();
 
        //Compute upper bound for cycle number.
        ComputeBoundOnCycleNumber<ComplexT>();
 
 
        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<RealT>::highest();
 
        // Default the selection before the candidate loop: if every candidate's difference is
        // NaN (poisoned samples), the comparisons below are all false and nothing is assigned --
        // a fresh endgame would otherwise carry cycle number 0 into TransformToSPlane and throw.
        this->cycle_number_ = 1;
 
        TimeCont<ComplexT> s_times(num_pts);
        SampCont<ComplexT> s_derivatives(num_pts);
 
        for(unsigned int candidate = 1; candidate <= upper_bound_on_cycle_number_; ++candidate)
        {
            using std::pow;
 
            std::tie(s_times, s_derivatives) = TransformToSPlane(static_cast<int>(candidate), t0, num_pts, ContStart::Front);
            RealT cand_power{1/static_cast<RealT>(candidate)};
            RealT curr_diff = (HermiteInterpolateAndSolve<ComplexT>(
                                  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 ComplexT>
    void ComputeAllDerivatives()
    {
        auto& samples = std::get<SampCont<ComplexT> >(samples_);
        auto& times   = std::get<TimeCont<ComplexT> >(times_);
        auto& derivatives = std::get<SampCont<ComplexT> >(derivatives_);
 
        assert((samples.size() == times.size()) && "must have same number of times and samples");
 
        if constexpr (tracking::TrackerTraits<TrackerType>::IsAdaptivePrec) // known at compile time
        {
            auto max_precision = this->EnsureAtUniformPrecision(times, samples);
        }
 
        //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 ComplexT>
    std::tuple<TimeCont<ComplexT>, SampCont<ComplexT>> TransformToSPlane(int cycle_num, ComplexT 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<ComplexT>();
 
 
        using RealT = typename Eigen::NumTraits<ComplexT>::Real;
 
        const auto& times   = std::get<TimeCont<ComplexT> >(times_);
        const auto& derivatives  = std::get<SampCont<ComplexT> >(derivatives_);
 
        RealT c = static_cast<RealT>(cycle_num);
        RealT one_over_c = 1/c;
 
        size_t 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<ComplexT> s_times(num_pts);
        SampCont<ComplexT> s_derivatives(num_pts);
 
        ComplexT 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 ComplexT>
    SuccessCode ComputeApproximationOfXAtT0(Vec<ComplexT>& result, const ComplexT & t0)
    {
        const auto c = ComputeCycleNumber<ComplexT>(t0);
 
        auto num_pts = this->EndgameSettings().num_sample_points;
 
        TimeCont<ComplexT> s_times;
        SampCont<ComplexT> s_derivatives;
 
        std::tie(s_times, s_derivatives) = TransformToSPlane(static_cast<int>(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(ComplexT(0), num_pts, s_times, std::get<SampCont<ComplexT> >(samples_), s_derivatives, ContStart::Back);
        // A NaN extrapolation must be a FAILURE code.  The run loop converges on
        // `approx_error > FinalTolerance()` becoming false, and every IEEE comparison against
        // NaN is false -- so a NaN approximation would exit the loop down the SUCCESS path,
        // reporting Converged with a poisoned answer.
        if (bertini::ContainsNaN(result))
            return SuccessCode::FailedToConverge;
        return SuccessCode::Success;
    }//end ComputeApproximationOfXAtT0
 
 

    template<typename ComplexT>
    SuccessCode AdvanceTime(const ComplexT & target_time)
    {
        using RealT = typename Eigen::NumTraits<ComplexT>::Real;
 
        auto& samples = std::get<SampCont<ComplexT> >(samples_);
        auto& times   = std::get<TimeCont<ComplexT> >(times_);
 
        AssertSizesTimeSpaceDeriv<ComplexT>();
 
        Vec<ComplexT> next_sample;
        ComplexT next_time = (times.back() + target_time) * static_cast<RealT>(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->EndgameTrackPath(next_sample,times.back(),next_time,samples.back());
            if (tracking_success != SuccessCode::Success)
                return tracking_success;
 
        // Pure-(i) escalation: return BEFORE pushing the new sample so the window stays an untouched
        // checkpoint; the driver migrates to mpfr and retries this advance.  Elided for fixed precision.
        if constexpr (tracking::TrackerTraits<TrackerType>::IsAdaptivePrec)
            if (this->adaptive_numeric_type_active_ &&
                this->GetTracker().GetCurrentPrecision() > this->current_endgame_precision_)
                return SuccessCode::HigherPrecisionNecessary;
 
        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_refinements);
        if (refine_success != SuccessCode::Success)
        {
            // Adaptive-numeric-type: a refine that double cannot satisfy is a request to cross to mpfr.
            // Roll back the just-pushed sample so the window is a clean checkpoint, then signal the driver.
            if constexpr (tracking::TrackerTraits<TrackerType>::IsAdaptivePrec)
                if (this->adaptive_numeric_type_active_ &&
                    (refine_success == SuccessCode::HigherPrecisionNecessary || refine_success == SuccessCode::FailedToConverge))
                {
                    times.pop_back();
                    samples.pop_back();
                    return SuccessCode::HigherPrecisionNecessary;
                }
            NotifyObservers(RefiningFailed<EmitterType>(*this));
            return refine_success;
        }
 
        this->EnsureAtPrecision(times.back(),Precision(samples.back()));
 
        NotifyObservers(SampleRefined<EmitterType>(*this));
        // the sample is complete once the new sample has been refined -- emit the
        // REFINED value, not the freshly-tracked one pushed above
        // ComputedSamplePoint carries the sample/time at BaseComplexT, like CircleAdvanced.
        // The fixed/mpfr lane emits directly; the complex_dbl fast lane would have to convert
        // the point to mpfr for the event, so we only pay that when something is observing
        // (temporaries live through the synchronous NotifyObservers).
        if constexpr (std::is_same<ComplexT, BCT>::value)
            NotifyObservers(ComputedSamplePoint<EmitterType>(*this, samples.back(), times.back()));
        else if (this->HasObservers())
        {
            Vec<BCT> ev_pt(samples.back().size());
            for (Eigen::Index i = 0; i < samples.back().size(); ++i) ev_pt(i) = BCT(samples.back()(i));
            BCT ev_t(times.back());
            NotifyObservers(ComputedSamplePoint<EmitterType>(*this, ev_pt, ev_t));
        }
 
        // 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();
        }
 
        // Has the path entered the endgame operating zone?  The same test the Cauchy endgame
        // uses, over the same kind of window, so the event asserts the same thing in both
        // (b2#402).  It used to fire on every advance here, which made it mean "advanced".
        NoteWhetherInOperatingZone<ComplexT>(samples);
 
        return SuccessCode::Success;
    }
 

    template<typename ComplexT>
    void NoteWhetherInOperatingZone(SampCont<ComplexT> const& samples) const
    {
        if (announced_operating_zone_ || samples.size() < 3)
            return;
 
        auto const needed = this->EndgameSettings().num_needed_for_stabilization;
        auto& c_over_k = std::get<TimeCont<ComplexT> >(c_over_k_);
        c_over_k.push_back(this->template ComputeCOverK<ComplexT>(samples));
        while (c_over_k.size() > needed + 1)
            c_over_k.pop_front();
 
        if (c_over_k.size() > needed && this->CheckForCOverKStabilization(c_over_k))
        {
            announced_operating_zone_ = true;
            NotifyObservers(InEGOperatingZone<EmitterType>(*this));
        }
    }
 

    template<typename ComplexT>
    SuccessCode RunImpl(const ComplexT & start_time, const Vec<ComplexT> & start_point, ComplexT const& target_time)
    {
        if (start_point.size()!=static_cast<Eigen::Index>(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());
        }
 
        // thread-local only: the endgame may run on a std::thread worker.
        SetThreadPrecision(Precision(start_point));
 
        using RealT = typename Eigen::NumTraits<ComplexT>::Real;
        //Set up for the endgame.
        ClearTimesAndSamples<ComplexT>();
 
 
        // unpack some references for easy use
        auto& samples = std::get<SampCont<ComplexT> >(samples_);
        auto& times   = std::get<TimeCont<ComplexT> >(times_);
        Vec<ComplexT>& latest_approx = this->final_approximation_;
        Vec<ComplexT>& prev_approx = this->previous_approximation_;
 
        // this is for estimating a ... norm?
        SetRandVec<ComplexT>(static_cast<int>(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<ComplexT>(samples, times);
        ComputeAllDerivatives<ComplexT>();
 
 
 
        auto extrapolation_code = ComputeApproximationOfXAtT0(prev_approx, target_time);
        latest_approx = prev_approx;
 
        if (extrapolation_code != SuccessCode::Success)
            return extrapolation_code;
 
        RealT norm_of_dehom_of_latest_approx(0); // initialized to 0 so the security check never reads indeterminate values
        RealT norm_of_dehom_of_prev_approx(0);
        if (this->SecuritySettings().level <= 0)
            norm_of_dehom_of_prev_approx = this->GetSystem().InfinityNormOfDehomogenized(prev_approx);
 
 
        NumErrorT& approx_error = this->approximate_error_;
        approx_error = 1;
 
        while (approx_error > this->FinalTolerance())
        {
            auto advance_code = AdvanceTime<ComplexT>(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<ComplexT>(samples, times);
            ComputeAllDerivatives<ComplexT>();
 
            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().InfinityNormOfDehomogenized(latest_approx);
                if(this->BeyondSecurityMaxNorm(norm_of_dehom_of_latest_approx) && this->BeyondSecurityMaxNorm(norm_of_dehom_of_prev_approx))
                {
                    NotifyObservers(SecurityMaxNormReached<EmitterType>(*this));
                    return SuccessCode::SecurityMaxNormReached;
                }
                norm_of_dehom_of_prev_approx = norm_of_dehom_of_latest_approx;
            }
 
 
 
 
            // CONVERGED.  Break BEFORE the assignment below so previous_approximation_ keeps
            // holding the genuine PREDECESSOR of final_approximation_.  The Cauchy flavor
            // already behaves this way -- it returns from its acceptance gate before the
            // corresponding assignment -- and the two must agree, because the public
            // accessors FinalApproximation()/PreviousApproximation()/ApproximateError() are
            // only coherent as a triple: overwriting here made previous == final on every
            // successful run, so the reported error described a pair the caller could not
            // see.  The loop's own condition is left in place so a FinalTolerance() >= 1
            // still refuses to enter at all, exactly as before.
            if (approx_error <= this->FinalTolerance())
                break;
 
            Precision(prev_approx, Precision(latest_approx));
            prev_approx = latest_approx;
        } //end while
 
        NotifyObservers(Converged<EmitterType>(*this));
        return SuccessCode::Success;
 
    } //end PSEG
 
 
    // ================================================================================================
    //   Adaptive-numeric-type (double-first) PowerSeries endgame.  Mirrors the Cauchy flavor: compute in
    //   the hardware complex_dbl fast lane while the tracker's authoritative precision stays double, and
    //   cross to complex_mp only when the tracker escalates (pure-(i)).  Member templates so the explicit
    //   fixed-precision class instantiations never force-compile them.  Fixed precision uses RunImpl<BCT>.
    // ================================================================================================

    template<typename Dummy = void>
    SuccessCode RunImplAMP(BCT const& start_time, Vec<BCT> const& start_point, BCT const& target_time)
    {
        using bertini::Precision;
        using RealT = typename Eigen::NumTraits<BCT>::Real;
 
        if (start_point.size()!=static_cast<Eigen::Index>(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());
        }
 
        this->adaptive_numeric_type_active_ = true;
        struct Disarmer { bool& flag; ~Disarmer(){ flag = false; } } disarm{this->adaptive_numeric_type_active_};
 
        // Both lanes start empty.  A segment clears only the lane it runs in, so a path whose
        // endgame stays in double would otherwise leave an earlier path's multiprecision times in
        // place -- and LatestTime reads that lane first whenever it is non-empty, reporting the
        // earlier path's last time as this one's (#378).
        ClearTimesAndSamples<complex_dbl>();
        ClearTimesAndSamples<complex_mp>();
 
        this->current_endgame_precision_ = std::max(DoublePrecision(), Precision(start_point));
 
        // ---- SETUP (restart-in-mpfr on escalation; bounded) ----
        while (true)
        {
            SuccessCode code = (this->current_endgame_precision_==DoublePrecision())
                ? SetupSegmentT<complex_dbl>(complex_dbl(start_time), this->DowncastToDouble(start_point), complex_dbl(target_time))
                : SetupSegmentT<complex_mp>(this->AtActivePrecisionScalar(start_time), this->AtActivePrecisionVec(start_point), this->AtActivePrecisionScalar(target_time));
            if (code == SuccessCode::HigherPrecisionNecessary)
            {
                auto const previous_precision = this->current_endgame_precision_;
                this->current_endgame_precision_ = this->NextEscalatedPrecision();
                SetThreadPrecision(this->current_endgame_precision_);
                // the samples announced at the lower precision are superseded (the window is
                // tracked again at the new one); say so, so a collector drops them as we do
                NotifyObservers(PrecisionChanged<EmitterType>(*this, previous_precision, this->current_endgame_precision_));
                NotifyObservers(SamplesRecomputedAtHigherPrecision<EmitterType>(*this));
                continue;
            }
            if (code != SuccessCode::Success)
                return code;
            break;
        }
 
        // previous_approximation_ (BCT) holds the first extrapolation, set by SetupSegmentT.
        RealT norm_prev(0), norm_latest(0);
        if (this->SecuritySettings().level <= 0)
            norm_prev = this->GetSystem().InfinityNormOfDehomogenized(this->previous_approximation_);
 
        this->approximate_error_ = 1;
 
        // ---- MAIN LOOP: each tracker-touching phase self-heals (migrate-and-retry in mpfr) ----
        while (this->approximate_error_ > this->FinalTolerance())
        {
            auto adv = AdvanceTimeAMP(target_time);
            if (adv != SuccessCode::Success) { NotifyObservers(EndgameFailure<EmitterType>(*this)); return adv; }
 
            auto ref = RefineAllSamplesAMP();
            if (ref != SuccessCode::Success) return ref;
 
            ComputeAllDerivativesAMP();
 
            auto ext = ComputeApproxAMP(target_time);   // -> final_approximation_ (BCT)
            if (ext != SuccessCode::Success) { NotifyObservers(EndgameFailure<EmitterType>(*this)); return ext; }
 
            Precision(this->previous_approximation_, Precision(this->final_approximation_));
            this->approximate_error_ = static_cast<NumErrorT>((this->final_approximation_ - this->previous_approximation_).template lpNorm<Eigen::Infinity>());
            NotifyObservers(ApproximatedRoot<EmitterType>(*this));
 
            if (this->SecuritySettings().level <= 0)
            {
                norm_latest = this->GetSystem().InfinityNormOfDehomogenized(this->final_approximation_);
                if (this->BeyondSecurityMaxNorm(norm_latest) && this->BeyondSecurityMaxNorm(norm_prev))
                {
                    NotifyObservers(SecurityMaxNormReached<EmitterType>(*this));
                    return SuccessCode::SecurityMaxNormReached;
                }
                norm_prev = norm_latest;
            }
 
            // CONVERGED.  Break BEFORE the assignment below so previous_approximation_ keeps
            // holding the genuine PREDECESSOR of final_approximation_.  The Cauchy flavor
            // already behaves this way -- it returns from its acceptance gate before the
            // corresponding assignment -- and the two must agree, because the public
            // accessors FinalApproximation()/PreviousApproximation()/ApproximateError() are
            // only coherent as a triple: overwriting here made previous == final on every
            // successful run, so the reported error described a pair the caller could not
            // see.  The loop's own condition is left in place so a FinalTolerance() >= 1
            // still refuses to enter at all, exactly as before.
            if (this->approximate_error_ <= this->FinalTolerance())
                break;
 
            this->previous_approximation_ = this->final_approximation_;
        }
 
        NotifyObservers(Converged<EmitterType>(*this));
        return SuccessCode::Success;
    }
 
 
    // Pre-loop work + first extrapolation, at one numeric type.  Reports HigherPrecisionNecessary up to
    // the driver (which restarts setup in mpfr) on escalation.
    template<typename ComplexT>
    SuccessCode SetupSegmentT(ComplexT const& start_time, Vec<ComplexT> const& start_point, ComplexT const& target_time)
    {
        SetThreadPrecision(Precision(start_point));
        ClearTimesAndSamples<ComplexT>();
 
        auto& samples = std::get<SampCont<ComplexT> >(samples_);
        auto& times   = std::get<TimeCont<ComplexT> >(times_);
 
        SetRandVec<ComplexT>(static_cast<int>(start_point.size()));
 
        auto init = this->ComputeInitialSamples(start_time, target_time, start_point, times, samples);
        if (init != SuccessCode::Success) { NotifyObservers(EndgameFailure<EmitterType>(*this)); return init; }
        if (this->GetTracker().GetCurrentPrecision() > this->current_endgame_precision_)
            return SuccessCode::HigherPrecisionNecessary;
 
        auto ref = this->template RefineAllSamples<ComplexT>(samples, times);
        if (ref == SuccessCode::HigherPrecisionNecessary || ref == SuccessCode::FailedToConverge)
            return SuccessCode::HigherPrecisionNecessary;
        if (ref != SuccessCode::Success) return ref;
 
        ComputeAllDerivatives<ComplexT>();
 
        Vec<ComplexT> first_approx;
        auto ext = ComputeApproximationOfXAtT0<ComplexT>(first_approx, target_time);
        if (ext != SuccessCode::Success) return ext;
        this->previous_approximation_ = this->ToBCT(first_approx, this->current_endgame_precision_);
        return SuccessCode::Success;
    }
 

    template<typename Dummy = void>
    SuccessCode AdvanceTimeAMP(BCT const& target_time)
    {
        unsigned guard = 0;
        while (true)
        {
            SuccessCode code = (this->current_endgame_precision_==DoublePrecision())
                ? AdvanceTime<complex_dbl>(complex_dbl(target_time))
                : AdvanceTime<complex_mp>(this->AtActivePrecisionScalar(target_time));
            if (code == SuccessCode::HigherPrecisionNecessary)
            {
                if (this->template EscalateAndMigrate<>(++guard) != SuccessCode::Success) return SuccessCode::HigherPrecisionNecessary;
                continue;
            }
            return code;
        }
    }

    template<typename Dummy = void>
    SuccessCode RefineAllSamplesAMP()
    {
        unsigned guard = 0;
        while (true)
        {
            SuccessCode code = (this->current_endgame_precision_==DoublePrecision())
                ? this->template RefineAllSamples<complex_dbl>(std::get<SampCont<complex_dbl> >(samples_), std::get<TimeCont<complex_dbl> >(times_))
                : this->template RefineAllSamples<complex_mp>(std::get<SampCont<complex_mp> >(samples_), std::get<TimeCont<complex_mp> >(times_));
            if (code == SuccessCode::HigherPrecisionNecessary || code == SuccessCode::FailedToConverge)
            {
                if (this->template EscalateAndMigrate<>(++guard) != SuccessCode::Success) return SuccessCode::HigherPrecisionNecessary;
                continue;
            }
            return code;
        }
    }

    template<typename Dummy = void>
    void ComputeAllDerivativesAMP()
    {
        if (this->current_endgame_precision_==DoublePrecision())
            ComputeAllDerivatives<complex_dbl>();
        else
            ComputeAllDerivatives<complex_mp>();
    }

    template<typename Dummy = void>
    SuccessCode ComputeApproxAMP(BCT const& target_time)
    {
        if (this->current_endgame_precision_==DoublePrecision())
        {
            Vec<complex_dbl> r;
            auto code = ComputeApproximationOfXAtT0<complex_dbl>(r, complex_dbl(target_time));
            if (code == SuccessCode::Success)
                this->final_approximation_ = this->ToBCT(r, this->current_endgame_precision_);
            return code;
        }
        else
        {
            return ComputeApproximationOfXAtT0<complex_mp>(this->final_approximation_, this->AtActivePrecisionScalar(target_time));
        }
    }
 
 
    // Flavor-specific: cross every PowerSeries container from the complex_dbl slot to complex_mp (or, if
    // already mpfr, raise its precision uniformly), via the shared AMP-policy Cross* / SetPrecision
    // helpers.  Called by the base EscalateAndMigrate.
    template<typename Dummy = void>
    void MigrateContainersToPrecision(unsigned newprec)
    {
        using bertini::Precision;
        if (this->current_endgame_precision_ == DoublePrecision())
        {
            this->CrossTimesUp(times_,       newprec);
            this->CrossSampsUp(samples_,     newprec);
            this->CrossSampsUp(derivatives_, newprec);
            this->CrossVecUp  (rand_vector_, newprec);
        }
        else
        {
            tracking::adaptive::SetPrecision(std::get<TimeCont<complex_mp> >(times_),       newprec);
            tracking::adaptive::SetPrecision(std::get<SampCont<complex_mp> >(samples_),     newprec);
            tracking::adaptive::SetPrecision(std::get<SampCont<complex_mp> >(derivatives_), newprec);
            auto& rv = std::get<Vec<complex_mp> >(rand_vector_);
            if (rv.size() > 0) Precision(rv, newprec);
        }
        if (this->final_approximation_.size()    > 0) Precision(this->final_approximation_,    newprec);
        if (this->previous_approximation_.size() > 0) Precision(this->previous_approximation_, newprec);
    }
 
 
    virtual ~PowerSeriesEndgame() = default;
}; // end powerseries class
 
 
 
 
}} // re: namespaces