base_endgame.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
//base_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.
//
//base_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 base_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.
 
#ifndef BERTINI_BASE_ENDGAME_HPP
#define BERTINI_BASE_ENDGAME_HPP
 
#pragma once
 
#include <iostream>
#include <limits>   // std::numeric_limits, for approximate_error_'s infinity sentinel
#include <typeinfo>
 
 
#include "bertini2/mpfr_complex.hpp"
 
#include "bertini2/system/system.hpp"
 
#include "bertini2/detail/enable_permuted_arguments.hpp"
 
#include "bertini2/trackers/config.hpp"
#include "bertini2/trackers/adaptive_precision_utilities.hpp"  // tracking::adaptive::SetPrecision, for the shared container-migration helpers
#include "bertini2/endgames/config.hpp"
#include "bertini2/endgames/interpolation.hpp"
#include "bertini2/endgames/events.hpp"
 
#include "bertini2/logging.hpp"
 
#include "bertini2/detail/configured.hpp"
 
 
namespace bertini{ namespace endgame {
 

template<class FlavorT, class PrecT>
class EndgameBase :
    public detail::Configured< typename AlgoTraits<FlavorT>::NeededConfigs >,
    public PrecT,
    public virtual Observable
{
public:
    using TrackerType = typename PrecT::TrackerType;  
 
    using BaseComplexT = typename tracking::TrackerTraits<TrackerType>::BaseComplexT;  
    using BaseRealT = typename tracking::TrackerTraits<TrackerType>::BaseRealT;  
 
    using EmitterType = FlavorT;  

    using Config = detail::Configured< typename AlgoTraits<FlavorT>::NeededConfigs >;
 
protected:
 
    using BCT = BaseComplexT;  
    using BRT = BaseRealT;  
 
 
    using Configured = detail::Configured< typename AlgoTraits<FlavorT>::NeededConfigs >;  
    using Configs = typename AlgoTraits<FlavorT>::NeededConfigs;  
    using ConfigsAsTuple = typename Configs::ToTuple;  
 
    // The list of arithmetic types this endgame may compute in.  Sourced from the tracker, so it
    // matches the tracker's own dual/single-slot tuple machinery exactly:
    //   AMP        -> TypeList<complex_dbl, complex_mp>  (hardware-double fast lane + mpfr authority)
    //   fixed dbl  -> TypeList<complex_dbl>
    //   fixed mp   -> TypeList<complex_mp>
    // For AMP this makes every container (TupOfVec/TupleOfTimes/TupleOfSamps) dual-slot; the existing
    // type-indexed std::get<...ComplexT...>(member) access then works for complex_dbl for free.
    using NeededTypes = typename tracking::TrackerTraits<TrackerType>::NeededTypes;  
    using TupOfVec = typename NeededTypes::ToTupleOfVec;  
    using TupOfReal = typename NeededTypes::ToTupleOfReal;  
    using TupleOfTimes = typename NeededTypes::template ToTupleOfCont<TimeCont>;  
    using TupleOfSamps = typename NeededTypes::template ToTupleOfCont<SampCont>;  
 
 
 
    // universal endgame state variables
    mutable Vec<BCT> final_approximation_;       
    mutable Vec<BCT> previous_approximation_;     
    mutable unsigned int cycle_number_ = 0;       
    mutable NumErrorT approximate_error_ = std::numeric_limits<NumErrorT>::infinity();
 
    // The adaptive-numeric-type state (current_endgame_precision_, adaptive_numeric_type_active_) lives in
    // the AMP precision policy, AMPEndgame -- the flavors reach it through this-> (it is a base via PrecT).
 
    BCT start_time_{};   
    BCT target_time_{};  

    mutable Vec<complex_dbl> c_over_k_probe_as_drawn_;

    mutable TupOfVec c_over_k_probe_;
 

    // Drawn per run, never lazily on first use: a lazy draw lands in whichever path first needs
    // the probe, so the direction a path used, and the draws it consumed, would depend on the
    // paths the same endgame ran before it, and on which thread ran them (#378, ADR-0071).
    void RefreshCOverKProbe(unsigned size) const
    {
        c_over_k_probe_as_drawn_.resize(size);
        for (unsigned ii = 0; ii < size; ++ii)
            c_over_k_probe_as_drawn_(ii) = RandomUnit<complex_dbl>();
        std::apply([](auto&... working){ (working.resize(0), ...); }, c_over_k_probe_);
    }
 

    template<typename ComplexT>
    Vec<ComplexT> const& GetCOverKProbe(unsigned size, unsigned prec) const
    {
        using bertini::Precision;
        if (static_cast<unsigned>(c_over_k_probe_as_drawn_.size()) != size)
            RefreshCOverKProbe(size);   // samples of another size than this run began with
 
        auto& probe = std::get<Vec<ComplexT> >(c_over_k_probe_);
        bool stale = static_cast<unsigned>(probe.size()) != size;
        if constexpr (!std::is_same<ComplexT, complex_dbl>::value)
            stale = stale || Precision(probe) != prec;
        if (stale)
        {
            probe.resize(size);
            for (unsigned ii = 0; ii < size; ++ii)
                probe(ii) = ComplexT(c_over_k_probe_as_drawn_(ii));
            if constexpr (!std::is_same<ComplexT, complex_dbl>::value)
                Precision(probe, prec);
        }
        return probe;
    }
 
 
public:

    template<typename ComplexT>
    auto ComputeCOverK(SampCont<ComplexT> const& samples) const -> typename Eigen::NumTraits<ComplexT>::Real
    {
        using RealT = typename Eigen::NumTraits<ComplexT>::Real;
        using std::abs;
        using std::log;
 
        assert(samples.size()>=3);
        const Vec<ComplexT> & sample0 = samples[0];
        const Vec<ComplexT> & sample1 = samples[1];
        const Vec<ComplexT> & sample2 = samples[2];
 
        const Vec<ComplexT> & rand_vector =
            GetCOverKProbe<ComplexT>(static_cast<unsigned>(sample0.size()), Precision(sample0));
 
        // //DO NOT USE Eigen .dot() it will do conjugate transpose which is not what we want.
        // //Also, the .transpose*rand_vector returns an expression template that we do .norm of since abs is not available for that expression type.
        RealT estimate = abs(log(abs((((sample2 - sample1).transpose()*rand_vector).template lpNorm<Eigen::Infinity>())/(((sample1 - sample0).transpose()*rand_vector).template lpNorm<Eigen::Infinity>()))));
        estimate = abs(log(RealT(this->EndgameSettings().sample_factor)))/estimate;
        if (estimate < 1)
            return RealT(1);
        else
            return estimate;
    }
 

    template<typename T>
    bool CheckForCOverKStabilization(std::deque<T> const& c_over_k_array) const
    {
        using RealT = typename Eigen::NumTraits<T>::Real;
        using std::abs;
 
        assert(c_over_k_array.size()>=this->EndgameSettings().num_needed_for_stabilization);
        for(unsigned ii = 1; ii < this->EndgameSettings().num_needed_for_stabilization ; ++ii)
        {
            RealT a = abs(c_over_k_array[ii-1]);
            RealT b = abs(c_over_k_array[ii]);
 
            RealT divide = a;
            if(a < b)
                divide = a/b;
            else
                divide = b/a;
 
            if(divide < this->EndgameSettings().minimum_for_c_over_k_stabilization)
                return false;
        }
        return true;
    }
 
protected:

    const FlavorT& AsFlavor() const
    {
        return dynamic_cast<const FlavorT&>(*this);
    }

    FlavorT& AsFlavor()
    {
        return dynamic_cast<FlavorT&>(*this);
    }
 
public:

    void SetBoundaryTime(BCT const& t) { start_time_ = t; }
    void SetTargetTime  (BCT const& t) { target_time_ = t; }
    BCT const& BoundaryTime() const { return start_time_; }
    BCT const& TargetTime()   const { return target_time_; }

    SuccessCode Run(Vec<BCT> const& start_point)
    {
        using bertini::Precision;
        // EVERY path enters here, so this is the per-run boundary an observer needs to tell
        // one path's samples from the next -- without it, a collector attached to a solver's
        // endgame concatenates every path into one undifferentiated heap.  Initializing was
        // declared, and handled by GoryDetailLogger, but emitted by nothing.
        NotifyObservers(Initializing<EmitterType>(this->AsFlavor()));
        // a run's c/k direction is drawn here, at its start, never inherited from an earlier run
        RefreshCOverKProbe(static_cast<unsigned>(start_point.size()));
        auto prec = Precision(start_point);
        BCT t  = start_time_;   Precision(t,  prec);
        BCT t0 = target_time_;  Precision(t0, prec);
        // Fixed precision runs today's straight-line RunImpl<BCT> verbatim (untouched).  Adaptive
        // precision takes the separate double-first driver, RunImplAMP, which lives in the flavor.
        // RunImplAMP is a member template (template<typename=void>) so the explicit fixed-precision
        // class instantiations never force-compile it -- hence the .template disambiguator here.
        if constexpr (tracking::TrackerTraits<TrackerType>::IsAdaptivePrec)
            return this->AsFlavor().template RunImplAMP<>(t, start_point, t0);
        else
            return this->AsFlavor().RunImpl(t, start_point, t0);
    }
 

    template<typename ComplexT>
    void EmitComputedSample(Vec<ComplexT> const& sample, ComplexT const& time)
    {
        if constexpr (std::is_same<ComplexT, BCT>::value)
            NotifyObservers(ComputedSamplePoint<EmitterType>(AsFlavor(), sample, time));
        else if (this->HasObservers())
        {
            Vec<BCT> ev_pt(sample.size());
            for (Eigen::Index i = 0; i < sample.size(); ++i) ev_pt(i) = BCT(sample(i));
            BCT ev_t(time);
            NotifyObservers(ComputedSamplePoint<EmitterType>(AsFlavor(), ev_pt, ev_t));
        }
    }

    template<typename ComplexT>
    SuccessCode RefineAllSamples(SampCont<ComplexT> & samples, TimeCont<ComplexT> & times)
    {
        for (size_t ii=0; ii<samples.size(); ++ii)
        {
            auto refine_success = this->RefineSample(samples[ii], samples[ii],  times[ii],
                                        this->FinalTolerance() * this->EndgameSettings().sample_point_refinement_factor,
                                        this->EndgameSettings().max_num_refinements);
            if (refine_success != SuccessCode::Success)
            {
                // BOOST_LOG_TRIVIAL(severity_level::trace) << "refining failed, code " << int(refine_success);
                return refine_success;
            }
            NotifyObservers(SampleRefined<EmitterType>(AsFlavor()));
        }
 
        if constexpr (tracking::TrackerTraits<TrackerType>::IsAdaptivePrec) // known at compile time
        {
            auto max_precision = this->EnsureAtUniformPrecision(times, samples);
        }
 
        return SuccessCode::Success;
    }
 

    template<typename ComplexT>
    SuccessCode RefineSample(Vec<ComplexT> & result, Vec<ComplexT> const& current_sample, ComplexT const& current_time, NumErrorT tol, unsigned max_iterations) const
    {
        return this->RefineSampleImpl(result, current_sample, current_time, tol, max_iterations);
    }
 

    template<typename ComplexT>
    SuccessCode EndgameTrackPath(Vec<ComplexT> & result, ComplexT const& start_time, ComplexT const& end_time, Vec<ComplexT> const& start_point) const
    {
        if constexpr (std::is_same<ComplexT, BCT>::value)
        {
            return this->GetTracker().TrackPath(result, start_time, end_time, start_point);
        }
        else
        {
            using bertini::Precision;
            Vec<BCT> mp_start(start_point.size());
            for (Eigen::Index i = 0; i < start_point.size(); ++i) mp_start(i) = BCT(start_point(i));
            Precision(mp_start, DoublePrecision());
 
            BCT mp_t0(start_time); Precision(mp_t0, DoublePrecision());
            BCT mp_t1(end_time);   Precision(mp_t1, DoublePrecision());
 
            Vec<BCT> mp_result(start_point.size());
            auto code = this->GetTracker().TrackPath(mp_result, mp_t0, mp_t1, mp_start);
            if (code != SuccessCode::Success)
                return code;
 
            result.resize(mp_result.size());
            for (Eigen::Index i = 0; i < mp_result.size(); ++i) result(i) = complex_dbl(mp_result(i));
            return code;
        }
    }
 
 
    // EscalateAndMigrate bridges the AMP precision policy (NextEscalatedPrecision and the Cross* helpers,
    // all in AMPEndgame) with the flavor's container list (its MigrateContainersToPrecision).  It lives
    // here, not in AMPEndgame, because it needs AsFlavor(), which the precision-policy base lacks.  Member
    // template so the explicit fixed-precision class instantiations never force-compile it.  Raises the
    // working precision and migrates every durable container up to it; bounded by a runaway guard.
    template<typename Dummy = void>
    SuccessCode EscalateAndMigrate(unsigned guard)
    {
        if (guard > 64)
            return SuccessCode::HigherPrecisionNecessary;
        unsigned newprec = this->NextEscalatedPrecision();
        this->AsFlavor().template MigrateContainersToPrecision<>(newprec);
        this->current_endgame_precision_ = newprec;
        SetThreadPrecision(newprec);
        return SuccessCode::Success;
    }
 
 
    // note: a ChangePrecision(unsigned) lived here until 2026-06-12; it called
    // flavor-level ChangePrecision methods that have never existed, so it could
    // not compile -- it just was never instantiated until explicit template
    // instantiation (ADR-0014) forced every member.  zero callers; deleted.
    // precision changes go through the PrecT policy (see prec_base.hpp).
 

    inline
    const auto & EndgameSettings() const
    {
        return Configured::template Get<EndgameConfig>();
    }

    inline
    const auto & SecuritySettings() const
    {
        return this->template Get<SecurityConfig>();
    }

    template<typename RealT>
    bool BeyondSecurityMaxNorm(RealT const& norm) const
    {
        using std::isnan;
        return isnan(norm) || norm > static_cast<RealT>(this->SecuritySettings().max_norm);
    }

    explicit EndgameBase(TrackerType const& tr, const ConfigsAsTuple& settings ) :
        EndgamePrecPolicyBase<TrackerType>(tr), Configured( settings ), PrecT(tr)
    {}
 

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

    inline unsigned CycleNumber() const { return cycle_number_;}
    inline void CycleNumber(unsigned c) { cycle_number_ = c;}
    inline void IncrementCycleNumber(unsigned inc) { cycle_number_ += inc;}
 

    inline
    const auto& FinalTolerance() const
    {
        return this->template Get<EndgameConfig>().final_tolerance;
    }
 
 

    inline
    void SetFinalTolerance(NumErrorT const& ft)
    {
        // pre-ETI this member was never instantiated, hiding two defects: Get<>
        // returns const& (assignment through it cannot compile), and the old
        // BRT parameter type didn't match final_tolerance's storage (NumErrorT)
        auto settings = this->template Get<EndgameConfig>();
        settings.final_tolerance = ft;
        this->Set(settings);
    }
 
 
 

    template<typename ComplexT>
    inline
    const Vec<ComplexT>& FinalApproximation() const
    {
        return final_approximation_;
    }

    template<typename ComplexT>
    inline
    const Vec<ComplexT>& PreviousApproximation() const
    {
        return previous_approximation_;
    }

    inline
    NumErrorT ApproximateError() const
    {
        return approximate_error_;
    }
 

    inline
    const BCT& LatestTime() const
    {
        return AsFlavor().LatestTimeImpl();
    }
 
    // /**
    // \brief Get the system being tracked on, which is referred to by the tracker.
    // */
    // inline
    // const System& GetSystem() const
    // {
    //  return GetTracker().GetSystem();
    // }
 

    template<typename ComplexT>
    SuccessCode ComputeInitialSamples(const ComplexT & start_time,const ComplexT & target_time, const Vec<ComplexT> & x_endgame_start, TimeCont<ComplexT> & times, SampCont<ComplexT> & samples) // passed by reference to allow times to be filled as well.
    {
        using RealT = typename Eigen::NumTraits<ComplexT>::Real;
        assert(this->template Get<EndgameConfig>().num_sample_points>0 && "number of sample points must be positive");
 
        if constexpr (tracking::TrackerTraits<TrackerType>::IsAdaptivePrec)
        {
            assert(Precision(start_time)==Precision(x_endgame_start) && "Computing initial samples requires input time and space with uniform precision");
        }
 
        samples.clear();
        times.clear();
 
        samples.push_back(x_endgame_start);
        times.push_back(start_time);
        // The endgame's own input point, AT the boundary time: the coarsest sample of
        // the approach, and the widest |t| any consumer will ever see.  Without these
        // the sequence starts at boundary * sample_factor^(num_sample_points-1) --
        // measured, 0.1 -> 1.25e-2, throwing away most of the span before anyone
        // looks at it.
        EmitComputedSample(samples[0], times[0]);
 
        auto num_vars = this->GetSystem().NumVariables();
        //start at 1, because the input point is the 0th element.
        for(unsigned ii=1; ii < this->template Get<EndgameConfig>().num_sample_points; ++ii)
        {
            times.emplace_back((times[ii-1] + target_time) * RealT(this->template Get<EndgameConfig>().sample_factor)); // next time is a point between the previous time and target time.
            samples.emplace_back(Vec<ComplexT>(num_vars));                                             // sample_factor gives us some point between the two, usually the midpoint.
 
            auto tracking_success = this->EndgameTrackPath(samples[ii],times[ii-1],times[ii],samples[ii-1]);
            this->EnsureAtPrecision(times[ii],Precision(samples[ii]));
 
            if (tracking_success!=SuccessCode::Success)
                return tracking_success;
 
            EmitComputedSample(samples[ii], times[ii]);
        }
 
        return SuccessCode::Success;
    }
 
    virtual ~EndgameBase() = default;
};
 
} }// end namespaces bertini
 
 
#endif