zero_dim_solve.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
//bertini2/nag_algorithms/zero_dim_solve.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.
//
//bertini2/nag_algorithms/zero_dim_solve.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 bertini2/nag_algorithms/zero_dim_solve.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/num_traits.hpp"
#include "bertini2/eigen_extensions.hpp"   // bertini::IsDistinct (issues #302, #304)
 
#include "bertini2/detail/visitable.hpp"
#include "bertini2/tracking.hpp"
#include "bertini2/nag_algorithms/midpath_check.hpp"
#include "bertini2/io/generators.hpp"
 
#include "bertini2/detail/configured.hpp"
#include "bertini2/detail/observable.hpp"
 
#include "bertini2/nag_algorithms/common/algorithm_base.hpp"
#include "bertini2/nag_algorithms/common/config.hpp"
#include "bertini2/nag_algorithms/events.hpp"
#include "bertini2/system/start_base.hpp"   // start_system::StartSystem + StartSystemFactory / MakeStartFactory
#include "bertini2/parallel.hpp"
#include "bertini2/records/output_directory.hpp"
#include "bertini2/records/producer.hpp"
#include "bertini2/records/solver_recording.hpp"
#include "bertini2/records/config_encoding.hpp"
#include "bertini2/io/classic_writer.hpp"
#include "bertini2/io/json_writer.hpp"
#include <algorithm>
#include <chrono>
 
#include "bertini2/common/stop_request.hpp"
#include <cstdlib>
#include <limits>
#include <mutex>
#include <iostream>
#include <map>
 
 
namespace bertini {
 
    // forward-declare the interim default start system so ZeroDimSolver's default factory argument
    // (MakeStartFactory<TotalDegreeBinomial>) can name it; the concrete type rides in via start_systems.hpp
    // at every call site that actually constructs a ZeroDimSolver.
    namespace start_system { class TotalDegreeBinomial; }
 
    namespace algorithm {
 

template<typename TrackerType, typename EndgameType, typename SystemType>
struct HomotopySolver;

template<typename TrackerType, typename EndgameType, typename SystemType>
struct ZeroDimSolver;
 
 

template<typename TrackerType, typename EndgameType, typename SystemType>
struct AlgoTraits <HomotopySolver<TrackerType, EndgameType, SystemType>>
{
    using BaseRealT = typename tracking::TrackerTraits<TrackerType>::BaseRealT;        
    using BaseComplexT = typename tracking::TrackerTraits<TrackerType>::BaseComplexT;  

    using NeededConfigs = detail::TypeList<
                                TolerancesConfig,
                                PostProcessingConfig,
                                ZeroDimConfig,
                                AutoRetrackConfig,
                                RecordsConfig,
                                MidPathConfig
                                >;
};
 

struct AnyZeroDim : public virtual AnyAlgorithm
{
    virtual void WriteMainData(std::ostream& out) const = 0;
    virtual void WriteRawData(std::ostream& out)  const = 0;
    // Bertini 1.7-compatible machine-readable solution files (count-led coordinate blocks).
    virtual void WriteFiniteSolutions(std::ostream& out)      const = 0;
    virtual void WriteRealFiniteSolutions(std::ostream& out)  const = 0;
    virtual void WriteNonsingularSolutions(std::ostream& out) const = 0;
    virtual void WriteSingularSolutions(std::ostream& out)    const = 0;
    virtual void WriteRawSolutions(std::ostream& out)         const = 0;
    virtual void ApplyParsedConfigs(std::string const& config_str) = 0;
    virtual void RecordToPath(std::string const& /*path*/) {}
    virtual std::string RecordsRunIdentity() const { return {}; }
    virtual void AppendRecordJson(std::string const& /*record_json*/) {}
    virtual std::string PutRecordsDefinition(std::string const& /*content*/,
                                             std::string const& /*kind*/,
                                             std::string const& /*label*/ = {}) { return {}; }
    virtual ~AnyZeroDim() = default;
};
 
 
 

using SolnIndT = typename SolnCont<complex_dbl>::size_type;
 


struct AlgorithmMetaData
{
 
    SolnIndT number_path_failures = 0;   
    SolnIndT number_path_successes = 0;  
    SolnIndT number_paths_tracked = 0;   
 
    std::chrono::system_clock::time_point start_time;  
    std::chrono::microseconds elapsed_time;            
};
 

template<typename ComplexT>
struct SolutionMetaData
{
    using SolnIndT = typename SolnCont<ComplexT>::size_type;  
 
    // only vaguely metadata.  artifacts of randomness or ordering
    SolnIndT path_index;            
    SolnIndT solution_index;        

    bool precision_changed = false;          
    ComplexT time_of_first_prec_increase;    
    decltype(DefaultPrecision()) max_precision_used = 0;  
    double path_time_seconds = 0.0;          
    unsigned num_successful_steps = 0;       
    unsigned num_failed_steps = 0;           
    double wall_clock_limit_seconds = 0;     
    Vec<ComplexT> latest_path_point;

    SuccessCode pre_endgame_success_code = SuccessCode::NeverStarted;     
 

    NumErrorT condition_number;                 
    Vec<NumErrorT> singular_values;
    NumErrorT newton_residual;              
    ComplexT final_time_used;               
    // Two accuracy estimates, and each says in its NAME which coordinates it is in (ADR-0069).
    // Both are the distance between the endgame's last two approximations of the root.  The
    // internal one is what the endgame's convergence test compares against final_tolerance;
    // it reads as a number of digits.  The user one is the absolute error in the units of
    // the caller's own variables.  There is deliberately no field named plain
    // `accuracy_estimate`: a reader takes that to be in their own coordinates, and it was not.
    NumErrorT accuracy_estimate_internal_coords;    
    NumErrorT accuracy_estimate_user_coords;    
    unsigned cycle_num;                         
    // Honest precision/accuracy of the computed solution, as DIGIT COUNTS.  precision_digits is the
    // working precision the endgame actually finished in (DoublePrecision() ~16 for a path that stayed in
    // the adaptive-numeric-type endgame's hardware-double fast lane; the mpfr precision for one that
    // escalated).  accuracy_digits is how many of those digits are trustworthy, from the convergence
    // agreement: floor(-log10(accuracy_estimate_internal_coords)), clamped to [0, precision_digits].  Read together:
    // "computed in precision_digits digits, good to accuracy_digits of them."
    unsigned precision_digits = 0;  
    unsigned accuracy_digits = 0;   
    SuccessCode endgame_success_code = SuccessCode::NeverStarted;      
 

    NumErrorT function_residual;    
 
    int multiplicity = 1;       
    bool multiplicity_representative = true; 
    bool is_real = false;               
    bool is_finite = false;             
    bool is_singular = false;               
    bool crossing_unresolved = false;
    // nonsolution flag: a finite, successful endpoint that is NOT a solution of the actual target
    // system -- a nonsolution.  ZeroDimSolver sets this when it squares up an over-determined system:
    // the randomized square system has extraneous roots that satisfy the random combinations but not
    // the original equations.  Orthogonal to is_finite (a nonsolution is finite); the finite / real /
    // singular accessors exclude nonsolutions, and they are exposed on their own (Nonsolutions()).
    // Load-bearing for the regeneration cascade, which must identify and discard nonsolutions.
    bool is_nonsolution = false;  

    bool operator==(const SolutionMetaData<ComplexT> & other){
        bool result =
            this->path_index == other.path_index
             && this->solution_index == other.solution_index
             && this->precision_changed == other.precision_changed
             && this->time_of_first_prec_increase == other.time_of_first_prec_increase
             && this->max_precision_used == other.max_precision_used
             && this->pre_endgame_success_code == other.pre_endgame_success_code
             && this->condition_number == other.condition_number
             && this->newton_residual == other.newton_residual
             && this->final_time_used == other.final_time_used
             && this->accuracy_estimate_internal_coords == other.accuracy_estimate_internal_coords
             && this->accuracy_estimate_user_coords == other.accuracy_estimate_user_coords
             && this->cycle_num == other.cycle_num
             && this->precision_digits == other.precision_digits
             && this->accuracy_digits == other.accuracy_digits
             && this->endgame_success_code == other.endgame_success_code
             && this->function_residual == other.function_residual
             && this->multiplicity == other.multiplicity
             && this->multiplicity_representative == other.multiplicity_representative
             && this->is_real == other.is_real
             && this->is_finite == other.is_finite
             && this->is_singular == other.is_singular
             && this->is_nonsolution == other.is_nonsolution
             && this->crossing_unresolved == other.crossing_unresolved
             && this->num_successful_steps == other.num_successful_steps
             && this->num_failed_steps == other.num_failed_steps
             && this->wall_clock_limit_seconds == other.wall_clock_limit_seconds
             && this->singular_values.size() == other.singular_values.size()
             && (this->singular_values.size() == 0
                 || (this->singular_values.array() == other.singular_values.array()).all())
             && this->latest_path_point.size() == other.latest_path_point.size()
             && (this->latest_path_point.size() == 0
                 || (this->latest_path_point.array() == other.latest_path_point.array()).all())
        ;
 
        return result; }
};

// this is for interoperability with vectors of these in the Python bindings, for better or for worse.
template<typename NumT>
std::ostream& operator<<(std::ostream & out, const SolutionMetaData<NumT> & meta){
    out << "path_index = " << meta.path_index << std::endl;
    out << "solution_index = " << meta.solution_index << std::endl;
 
    out << "precision_changed = " << meta.precision_changed << std::endl;
    out << "time_of_first_prec_increase = " << meta.time_of_first_prec_increase << std::endl;
    out << "max_precision_used = " << meta.max_precision_used << std::endl;
    out << "path_time_seconds = " << meta.path_time_seconds << std::endl;
    out << "num_successful_steps = " << meta.num_successful_steps << std::endl;
    out << "num_failed_steps = " << meta.num_failed_steps << std::endl;
    out << "wall_clock_limit_seconds = " << meta.wall_clock_limit_seconds << std::endl;
    out << "latest_path_point = " << meta.latest_path_point.transpose() << std::endl;
 
    out << "pre_endgame_success_code = " << meta.pre_endgame_success_code << std::endl;
 
    out << "condition_number = " << meta.condition_number << std::endl;
    out << "singular_values = " << meta.singular_values.transpose() << std::endl;
    out << "newton_residual = " << meta.newton_residual << std::endl;
    out << "final_time_used = " << meta.final_time_used << std::endl;
    out << "accuracy_estimate_internal_coords = " << meta.accuracy_estimate_internal_coords << std::endl;
    out << "accuracy_estimate_user_coords = " << meta.accuracy_estimate_user_coords << std::endl;
    out << "precision_digits = " << meta.precision_digits << std::endl;
    out << "accuracy_digits = " << meta.accuracy_digits << std::endl;
    out << "cycle_num = " << meta.cycle_num << std::endl;
    out << "endgame_success_code = " << meta.endgame_success_code << std::endl;
 
    out << "function_residual = " << meta.function_residual << std::endl;
 
    out << "multiplicity = " << meta.multiplicity << std::endl;
    out << "multiplicity_representative = " << meta.multiplicity_representative << std::endl;
    out << "is_real = " << meta.is_real << std::endl;
    out << "is_finite = " << meta.is_finite << std::endl;
    out << "is_singular = " << meta.is_singular << std::endl;
    out << "is_nonsolution = " << meta.is_nonsolution << std::endl;
    out << "crossing_unresolved = " << meta.crossing_unresolved << std::endl;
 
    return out;
}
 

template<typename ComplexT>
struct EGBoundaryMetaData
{
    using RealT = typename NumTraits<ComplexT>::Real;  
 
    Vec<ComplexT> path_point;                            
    SuccessCode success_code = SuccessCode::NeverStarted; 
    RealT last_used_stepsize;                            
    // The precision the tracker was actually using when it reached the endgame boundary.
    // Carried explicitly (rather than inferred from Precision(path_point)) because the
    // path_point is always stored as the tracker's BaseComplexT (multiprecision for AMP);
    // a path that tracked in double gets widened on output, so its mantissa precision no
    // longer reflects the precision that was in use.  The endgame should resume at this
    // precision.  See zero_dim_solve TrackSinglePathDuringEG.
    unsigned precision = DoublePrecision();  

    EGBoundaryMetaData() = default;
    EGBoundaryMetaData(EGBoundaryMetaData const&) = default;
    EGBoundaryMetaData& operator=(EGBoundaryMetaData const&) = default;
    EGBoundaryMetaData(Vec<ComplexT> const& pt, SuccessCode const& code, RealT const& ss, unsigned prec) :
        path_point(pt), success_code(code), last_used_stepsize(ss), precision(prec)
    {}

    bool operator==(const EGBoundaryMetaData<ComplexT> & other){
        bool result =
            this->path_point == other.path_point
            && this->success_code == other.success_code
            && this->last_used_stepsize == other.last_used_stepsize
            && this->precision == other.precision
        ;
 
        return result;
    }
};

// this is for interoperability with vectors of these in the Python bindings, for better or for worse.
template<typename NumT>
std::ostream& operator<<(std::ostream & out, const EGBoundaryMetaData<NumT> & meta){
    out << "path_point = " << meta.path_point << std::endl;
    out << "success_code = " << meta.success_code << std::endl;
    out << "last_used_stepsize = " << meta.last_used_stepsize << std::endl;
    out << "precision = " << meta.precision << std::endl;
    return out;
}

struct MidpathCheckReport
{
    bool passed = true;                                    
    unsigned num_crossings_detected = 0;                   
    unsigned num_resolve_attempts = 0;                     
    std::vector<unsigned long long> crossed_path_indices;  
};

inline
std::ostream& operator<<(std::ostream & out, const MidpathCheckReport & r){
    out << "passed = " << std::boolalpha << r.passed << std::endl;
    out << "num_crossings_detected = " << r.num_crossings_detected << std::endl;
    out << "num_resolve_attempts = " << r.num_resolve_attempts << std::endl;
    out << "crossed_path_indices = [";
    for (size_t i = 0; i < r.crossed_path_indices.size(); ++i)
        out << (i ? ", " : "") << r.crossed_path_indices[i];
    out << "]" << std::endl;
    return out;
}

struct SolveReport
{
    unsigned long long num_paths_tracked = 0;    
    unsigned long long num_finite_solutions = 0; 
    unsigned long long num_finite_endpoints = 0; 
    unsigned long long num_diverged = 0;         
    unsigned long long num_failed = 0;           
    unsigned long long num_singular = 0;         
    unsigned long long num_real = 0;             
    unsigned long long num_nonsolutions = 0;     
    std::map<SuccessCode, unsigned long long> failures_by_reason; 
    double max_condition_number = 0;             
    unsigned max_precision_used = 0;             
    MidpathCheckReport midpath;                  
    bool all_paths_resolved = false;             
};

template<typename ComplexT>
SolveReport SummarizeSolve(std::vector<SolutionMetaData<ComplexT>> const& metadata,
                           MidpathCheckReport const& midpath)
{
    SolveReport r;
    r.num_paths_tracked = metadata.size();
    r.midpath = midpath;
    double finite_distinct = 0;
    for (auto const& m : metadata)
    {
        if (m.max_precision_used > r.max_precision_used)
            r.max_precision_used = m.max_precision_used;
 
        // A path the endgame's security check truncated near infinity (SecurityMaxNormReached) is
        // diverging, not failing -- it is the security/max-norm heuristic catching a path on its way
        // to infinity, exactly like a clean GoingToInfinity.  (The endgame tests classify the two
        // together too.)  With infinite-path truncation on by default, these are expected, so they
        // must not inflate num_failed / clear all_paths_resolved.
        bool diverged = (m.endgame_success_code == SuccessCode::GoingToInfinity
                      || m.endgame_success_code == SuccessCode::SecurityMaxNormReached);
        if (m.endgame_success_code == SuccessCode::Success)
        {
            if (m.is_nonsolution)
                ++r.num_nonsolutions;       // a nonsolution: finite, but not a solution of the target
            else if (m.is_finite)
            {
                ++r.num_finite_endpoints;
                finite_distinct += 1.0 / m.multiplicity;
                if (m.is_singular) ++r.num_singular;
                if (m.is_real)     ++r.num_real;
                if (static_cast<double>(m.condition_number) > r.max_condition_number)
                    r.max_condition_number = static_cast<double>(m.condition_number);
            }
            else
                diverged = true;
        }
 
        if (diverged)
            ++r.num_diverged;
        else if (m.endgame_success_code != SuccessCode::Success)   // neither a solution nor a clean divergence
        {
            ++r.num_failed;
            ++r.failures_by_reason[m.endgame_success_code];
        }
    }
    r.num_finite_solutions = static_cast<unsigned long long>(finite_distinct + 0.5);
    r.all_paths_resolved = (r.num_failed == 0) && midpath.passed;
    return r;
}

inline
std::ostream& operator<<(std::ostream & out, const SolveReport & r)
{
    out << "zero-dim solve -- " << r.num_paths_tracked << " paths tracked\n";
    out << "  finite solutions    " << r.num_finite_solutions << "   (distinct)\n";
    out << "  diverged            " << r.num_diverged << "\n";
    out << "  FAILED              " << r.num_failed;
    if (r.num_failed)
    {
        out << "    ";
        bool first = true;
        for (auto const& kv : r.failures_by_reason)
        {
            out << (first ? "" : ", ") << kv.second << " x " << kv.first;
            first = false;
        }
    }
    out << "\n  ----\n";
    out << "  singular solutions  " << r.num_singular << "\n";
    out << "  real solutions      " << r.num_real << "\n";
    if (r.num_nonsolutions)
        out << "  nonsolutions        " << r.num_nonsolutions << "   (finite, but not solutions of the target)\n";
    out << "  path crossings      " << r.midpath.num_crossings_detected
        << (r.midpath.passed ? " (resolved)" : " (UNRESOLVED)") << "\n";
    out << "  max condition num   " << r.max_condition_number << "\n";
    out << "  max precision used  " << r.max_precision_used << " digits\n";
    out << "  all paths resolved? " << (r.all_paths_resolved ? "yes" : "NO") << "\n";
    return out;
}

        template<typename TrackerType, typename EndgameType, typename SystemType>
        struct HomotopySolver :
                            public virtual AnyZeroDim,
                            public detail::Configured<
                                typename AlgoTraits< HomotopySolver<TrackerType, EndgameType, SystemType>>::NeededConfigs>
        {
            // these usings are for getters in python
            using TrackerT          = TrackerType;   
            using EndgameT          = EndgameType;   
            using SystemT           = SystemType;    
            // the engine sees the start system only through the polymorphic base
            using StartSystemT       = bertini::start_system::StartSystem;  
            using StartSystemBaseT   = bertini::start_system::StartSystem;  
 
            // This engine emits its lifecycle events on the AnyZeroDim base, so a
            // single observer type can watch any templated solver.  Accept observers
            // declared for AnyZeroDim (in addition to the exact concrete type).
            bool ObservableIsA(std::type_index t) const override
            {
                return Observable::ObservableIsA(t) || t == std::type_index(typeid(AnyZeroDim));
            }
 
 
 

            using BaseComplexT  = typename tracking::TrackerTraits<TrackerType>::BaseComplexT;  
            using BaseRealT     = typename tracking::TrackerTraits<TrackerType>::BaseRealT;     
 
            using PrecisionConfig   = typename tracking::TrackerTraits<TrackerType>::PrecisionConfig;  
 
            using SolnIndT          = typename SolnCont<BaseComplexT>::size_type;  
 

            using Config = detail::Configured<
                                typename AlgoTraits<HomotopySolver<TrackerType, EndgameType, SystemType>>::NeededConfigs>;
            using Config::Get;
 
 
            using Tolerances = TolerancesConfig;          
            using PostProcessing = PostProcessingConfig;  
            using ZeroDimConf = ZeroDimConfig;            
            using AutoRetrack = AutoRetrackConfig;        
 
 
            using EGBoundaryMetaDataT = EGBoundaryMetaData<BaseComplexT>;  
            using SolutionMetaDataT = SolutionMetaData<BaseComplexT>;      
 
 
// a few more using statements
 
            using MidpathType = MidpathChecker<BaseRealT, BaseComplexT, EGBoundaryMetaData<BaseComplexT>>;  
 
            // --- system storage ---------------------------------------------------------------
            // The engine holds the homotopy, start system, and target by REFERENCE; the caller
            // (a user, or ZeroDimSolver) owns them and guarantees they outlive the engine.  This
            // references, not owned: the engine never forms its systems (ZeroDimSolver / the user does).
        protected:
            std::reference_wrapper<const SystemType>        target_system_;  
            std::reference_wrapper<const StartSystemBaseT>  start_system_;   
            std::reference_wrapper<const SystemType>        homotopy_;       
 
            // Re-seat the system references.  ZeroDimSolver calls this after an MPI broadcast
            // re-materializes rank 0's authoritative systems into its owned storage.
            void ResetSystems(SystemType const& hom, StartSystemBaseT const& start, SystemType const& target)
            {
                homotopy_      = std::cref(hom);
                start_system_  = std::cref(start);
                target_system_ = std::cref(target);
            }
 
#ifdef BERTINI2_HAVE_MPI
            // Hook: install rank 0's authoritative systems on every rank before a distributed solve.
            // A user-supplied homotopy (plain HomotopySolver) is the user's responsibility across
            // ranks, so this is a no-op here; ZeroDimSolver, which owns its systems, overrides it.
            virtual void DistributeSystems(MPI_Comm /*comm*/) {}
#endif
 
        public:
            const SystemType&       TargetSystem() const { return target_system_.get(); }
            const StartSystemBaseT& StartSystem()  const { return start_system_.get();  }
            const SystemType&       Homotopy()     const { return homotopy_.get();      }
 
 


            HomotopySolver(SystemType const& target, StartSystemBaseT const& start, SystemType const& homotopy)
             : target_system_(std::cref(target)), start_system_(std::cref(start)), homotopy_(std::cref(homotopy)),
               tracker_(homotopy), endgame_(tracker_)
            {
                DefaultSetup();
            }
 

            struct ScopedSessionState
            {
                unsigned precision_;               
                std::mt19937 engine_;              
                records::DrawStream stream_;       

                ScopedSessionState()
                    : precision_(DefaultPrecision()),
                      engine_(ThreadEngine()),
                      stream_(records::ThreadDrawStream())
                {}

                ~ScopedSessionState()
                {
                    DefaultPrecision(precision_);
                    ThreadEngine() = engine_;
                    records::ThreadDrawStream() = stream_;
                }
                ScopedSessionState(ScopedSessionState const&) = delete;             
                ScopedSessionState& operator=(ScopedSessionState const&) = delete;  
            };

            void Run() override
            {
                // TODO(Phase 2 Python): AnyZeroDim lacks GetSolutions() -- solutions are only
                // accessible via the concrete type's member arrays or WriteMainData().
                // Resolve when designing Python bindings; the right approach depends on whether
                // Python holds AnyZeroDim or the concrete ZeroDim<...> type directly.
#ifdef BERTINI2_HAVE_MPI
                if (parallel::Size() > 1)
                {
                    RunParallel(parallel::WorldComm());
                    return;
                }
#endif
                Solve();
            }
 
#ifdef BERTINI2_HAVE_MPI
            void RunParallel(MPI_Comm comm)
            {
                using Result = parallel::FullPathResult<BaseComplexT>;
                using Task   = parallel::StartPointTask<BaseComplexT>;
 
                ScopedSessionState const session_state_guard{};   // a solve leaves the session as it found it
 
                mpfr_free_cache(); // reproducibility: clear this rank's mpfr constant cache (see Solve())
 
                solutions_user_coords_fresh_ = false;
 
                // Each rank built its homotopy independently in the constructor -- with its OWN
                // random patch, start-system coefficients, and gamma -- so a given path index would
                // mean a different path on each worker, and the manager would gather a scrambled,
                // mostly-wrong solution set.  Make rank 0 the single authoritative source: broadcast
                // its actual target system, start system, and homotopy to every rank.  These are all
                // exact (rational coefficients + integer/rational patch), so serialization carries
                // them bit-for-bit -- there is no per-rank re-derivation that could drift, and the
                // distributed homotopy is identical to rank 0's serial one.  We also broadcast rank
                // 0's RNG seed so any randomness drawn later (e.g. in the endgame) is consistent; the
                // per-path tracking RNG is already deterministic in the path index (ReseedThisThread).
                {
                    unsigned long seed = parallel::IsManager() ? GetGlobalSeed() : 0ul;
                    MPI_Bcast(&seed, 1, MPI_UNSIGNED_LONG, 0, comm);
                    SetGlobalSeed(seed);
 
                    // Install rank 0's authoritative systems on every rank.  The plain engine over a
                    // user-supplied homotopy does nothing here (the user owns cross-rank consistency);
                    // ZeroDimSolver, which owns its systems, overrides DistributeSystems to broadcast
                    // rank 0's target / start / homotopy and re-seat the engine's references.
                    DistributeSystems(comm);
 
                    num_start_points_ = StartSystem().NumStartPoints();
                    GetTracker().SetSystem(Homotopy());
                }
 
                PreSolveChecks();
                PreSolveSetup();
 
                // threads per rank: the same rule as the standalone solve (BERTINI_NUM_THREADS,
                // then num_threads, then every CPU this rank may run on)
                const int n_threads = static_cast<int>(
                    parallel::EffectiveThreadCount(this->template Get<ZeroDimConf>().num_threads));
 
                // One round = dispatch a set of path indices, each executed as a WHOLE path
                // (pre-endgame + endgame) by a worker.  Round 0 is every path; later rounds re-run
                // only the crossed paths with escalated settings -- the re-track is path-independent,
                // so it rides the same demand-driven worker pool (no idle-manager bottleneck).
                if (parallel::IsManager())
                {
                    // Rank 0 computes every start point once, authoritatively, and ships each to the
                    // worker that tracks it.  Cached so re-track rounds reuse the identical points.
                    std::vector<Vec<BaseComplexT>> start_points(num_start_points_);
                    for (decltype(num_start_points_) ii{0}; ii < num_start_points_; ++ii)
                        start_points[ii] = ComputeStartPoint(static_cast<SolnIndT>(ii));
 
                    auto run_round = [&](std::queue<Task>& q)
                    {
                        parallel::RunManagerLoop<Task, Result>(comm, q,
                            [this](Result const& r){ StoreFullPathResult(r); });
                    };
 
                    // the records seam (ADR-0046): the manager is the sole writer; recall what
                    // is already recorded and dispatch only the rest
                    MaybeAttachAmbientRecords();
                    std::vector<SolnIndT> indices_to_run;
                    for (decltype(num_start_points_) ii{0}; ii < num_start_points_; ++ii)
                        indices_to_run.push_back(static_cast<SolnIndT>(ii));
                    if (records_)
                    {
                        EnsureRunRecorded();
                        indices_to_run = RecallRecordedPaths(indices_to_run);
                    }
 
                    std::queue<Task> queue;
                    for (auto idx : indices_to_run)
                        queue.push(Task{ idx, start_points[idx] });
                    run_round(queue);
 
                    // Midpath/crossing check on the collected boundary points, then bounded parallel
                    // resolve rounds.  A continue/stop flag is broadcast to the workers each round so
                    // they know whether to run again (and escalate their settings first).
                    auto passed = midpath_.Check(solutions_at_endgame_boundary_, StartSystem());
                    midpath_report_ = MidpathCheckReport{};
                    midpath_report_.passed = passed;
                    for (auto const& v : midpath_.GetCrossedPaths())
                        midpath_report_.crossed_path_indices.push_back(v.index());
                    midpath_report_.num_crossings_detected =
                        static_cast<unsigned>(midpath_report_.crossed_path_indices.size());
 
                    const auto max_attempts = this->template Get<ZeroDimConf>().max_num_crossed_path_resolve_attempts;
                    unsigned num_resolve_attempts = 0;
                    while (true)
                    {
                        int keep_going = (!passed && num_resolve_attempts < max_attempts) ? 1 : 0;
                        MPI_Bcast(&keep_going, 1, MPI_INT, 0, comm);
                        if (!keep_going)
                            break;
 
                        std::queue<Task> redo;
                        for (auto const& v : midpath_.GetCrossedPaths())
                            if (v.rerun())
                            {
                                auto idx = static_cast<SolnIndT>(v.index());
                                redo.push(Task{ idx, start_points[idx] });
                            }
                        run_round(redo);
 
                        passed = midpath_.Check(solutions_at_endgame_boundary_, StartSystem());
                        ++num_resolve_attempts;
                    }
                    midpath_report_.num_resolve_attempts = num_resolve_attempts;
                    midpath_report_.passed = passed;
                    FlagUnresolvedCrossings(passed);
                    RecordMidpathOutcome();
                    if (!passed)
                        std::cerr << "warning: " << midpath_report_.num_crossings_detected
                                  << " path crossing(s) detected at the endgame boundary remained unresolved "
                                  << "after " << num_resolve_attempts << " re-track attempt(s); "
                                  << "the affected solutions may be wrong.  Consider a higher-order predictor "
                                  << "or a tighter tracking tolerance.  See EndgameBoundaryMetadata()." << std::endl;
 
                    PostEGAction();
 
                    DeclareFiniteSolutionsResult();
                }
                else // worker
                {
                    auto run_worker_round = [&]()
                    {
                        if (n_threads <= 1)
                        {
                            // Single-threaded worker: execute whole paths on the rank's member
                            // tracker/endgame directly, from rank 0's authoritative start point.
                            parallel::RunWorkerLoop<Task, Result>(comm,
                                [this](Task const& task){ ExecuteOnePath(MemberDuringEGContext(), static_cast<SolnIndT>(task.path_index), task.start_point); },
                                [this](Task const& task) -> Result { return PackFullPathResult(static_cast<SolnIndT>(task.path_index)); });
                        }
                        else
                        {
                            // Threaded worker: each thread owns a self-contained PathThreadState clone
                            // so concurrent whole-path execution shares no mutable state.  Cloned per
                            // round so it inherits any escalation (predictor bump) applied to the rank's
                            // member tracker between rounds.
                            auto state_factory = [this]() -> std::unique_ptr<PathThreadState>
                            {
                                // Clone(), not copy: System copies are SHALLOW (shared expression-tree
                                // nodes whose value caches mutate on Eval); Clone() deep-copies the tree.
                                // Trackers/Endgames have no default ctor, so aggregate-init via new.
                                std::unique_ptr<PathThreadState> s(
                                    new PathThreadState{ Clone(GetTracker().GetSystem()), Clone(TargetSystem()),
                                                         GetTracker(), GetEndgame(), {} });
                                s->tracker.SetSystem(s->sys);
                                s->endgame.SetTracker(s->tracker);
                                SetThreadPrecision(this->template Get<ZeroDimConf>().initial_ambient_precision);
                                return s;
                            };
 
                            auto track_fn = [this](std::unique_ptr<PathThreadState>& state, Task const& task) -> Result
                            {
                                auto idx = static_cast<SolnIndT>(task.path_index);
                                ExecuteOnePath(state->Context(), idx, task.start_point);
                                return PackFullPathResult(idx);
                            };
 
                            parallel::RunWorkerLoopThreaded<Task, Result>(
                                comm, state_factory, track_fn, n_threads);
                        }
                    }; // run_worker_round
 
                    while (true)
                    {
                        run_worker_round();
 
                        int keep_going = 0;
                        MPI_Bcast(&keep_going, 1, MPI_INT, 0, comm);
                        if (!keep_going)
                            break;
 
                        // Match the escalation the serial RunMidpathResolution applies between
                        // re-track passes: tighten midpath_retrack_tolerance_ (read by ExecuteOnePath)
                        // and bump a low-order predictor on this rank's member tracker, so the next
                        // round's thread clones inherit it.
                        EscalateRetrackSettings();
                    }
                }
            }
#endif // BERTINI2_HAVE_MPI
 
            // Definitions are out-of-line at the bottom of this file, after
            // output.hpp is included (avoiding a circular-include chicken-and-egg).
            void WriteMainData(std::ostream& out) const override;
            void WriteRawData(std::ostream& out)  const override;
            void WriteFiniteSolutions(std::ostream& out)      const override;
            void WriteRealFiniteSolutions(std::ostream& out)  const override;
            void WriteNonsingularSolutions(std::ostream& out) const override;
            void WriteSingularSolutions(std::ostream& out)    const override;
            void WriteRawSolutions(std::ostream& out)         const override;
            void ApplyParsedConfigs(std::string const& config_str) override;
 
            virtual ~HomotopySolver() = default;
            // NOTE: feasibility checking (ConsistencyCheck) and start-system construction live in
            // ZeroDimSolver, which owns the systems; the engine is handed a ready homotopy.
 
 
            void DefaultSetup()
            {
                DefaultSettingsSetup();
                DefaultSystemSetup();
                DefaultTrackerSetup();
                DefaultMidpathSetup();
            }
 
 
 
 
 

            void DefaultSystemSetup()
            {
                // The systems are already built and referenced (the engine does not form them);
                // just read off the number of start points the start system will produce.
                num_start_points_ = StartSystem().NumStartPoints(); // populate the internal variable
            }
 
 

            void DefaultSettingsSetup()
            {
                // this code can be made generic using Boost.Hana.
                // see https://stackoverflow.com/questions/28764085/how-to-create-an-element-for-each-type-in-a-typelist-and-add-it-to-a-vector-in-c,
                // for example
                //
                // all that would need to be done is to extract a hana::tuple_t from the
                // typelist contained in this::Config, and then hana::for_each() over it.
 
                this->template Set<Tolerances>(Tolerances());
                this->template Set<PostProcessing>(PostProcessing());
                // ZeroDimConfig is not templated on the complex type, so it cannot pick its own
                // per-type ambient-precision default.  Set it here, where the tracking type is known.
                ZeroDimConf zdc;
                zdc.initial_ambient_precision = DefaultInitialAmbientPrecision<BaseComplexT>();
                this->template Set<ZeroDimConf>(zdc);
                this->template Set<AutoRetrack>(AutoRetrack());
            }

            void SetMidpathRetrackTol(NumErrorT const& rt)
            {
                midpath_retrack_tolerance_ = rt;
            }

            const auto& MidpathRetrackTol() const
            {
                return midpath_retrack_tolerance_;
            }

            void SetStartPointPrecision(unsigned p)
            {
                start_point_precision_ = p;
                start_point_precision_set_by_user_ = true;
            }

            unsigned StartPointPrecision() const
            {
                return start_point_precision_;
            }
 

            void DefaultMidpathSetup()
            {
                midpath_ = MidpathType(MidPathConfig());
            }
 

            void SetMidpath(MidPathConfig const& mp)
            {
                this->template Set<MidPathConfig>(mp);
                midpath_.Set(mp);
            }

            void DefaultTrackerSetup()
            {
                tracker_.SetSystem(Homotopy());
                // the solver sets what it owns -- the phase tracking tolerance -- and leaves the
                // tracker's own settings alone, so a predictor or truncation threshold set on the
                // tracker survives DefaultSetup
                tracker_.Setup(tracker_.GetPredictor(),
                                this->template Get<Tolerances>().newton_before_endgame,
                                tracking::SteppingConfig(), tracking::NewtonConfig());
 
                tracker_.PrecisionSetup(PrecisionConfig(Homotopy()));
            }
 

            void SetTracker(TrackerType const& new_tracker)
            {
                tracker_ = new_tracker;
                endgame_.SetTracker(tracker_);
            }

            const TrackerType & GetTracker() const
            {
                return tracker_;
            }

            TrackerType & GetTracker()
            {
                return tracker_;
            }
 


            void SetEndgame(EndgameType const& new_endgame)
            {
                endgame_ = new_endgame;
            }

            const EndgameType & GetEndgame() const
            {
                return endgame_;
            }

            EndgameType & GetEndgame()
            {
                return endgame_;
            }
 

            template<typename T>
            const T & GetFromEndgame() const
            {
                return endgame_.template Get<T>();
            }

            template<typename T>
            void SetToEndgame(T const& t)
            {
                endgame_.Set(t);
            }

 
 

            void Solve()
            {
                // Reproducibility: clear this thread's mpfr constant cache (pi, etc.) at the solve
                // boundary.  mpfr caches transcendental constants at the precision last requested; the
                // Cauchy endgame needs pi (roots of unity), so a high-precision solve leaves pi cached
                // at high precision.  A subsequent lower-precision solve would otherwise reuse that
                // cached value rounded down -- differing by ~1 ULP from a freshly-computed
                // low-precision pi -- making a solve's behaviour depend on what ran before it in the
                // process (non-reproducible).  Clearing here makes each solve independent of history.
                // Cheap (one constant recompute per solve); thread-local (each worker thread clears its
                // own when it starts a solve).
                mpfr_free_cache();
 
                ScopedSessionState const session_state_guard{};   // a solve leaves the session as it found it
 
                solutions_user_coords_fresh_ = false;
 
                PreSolveChecks();
 
                PreSolveSetup();
 
                // the whole-solve wall-clock budget, counted from here; none unless asked for
                solve_deadline_.reset();
                if (auto const limit = this->template Get<ZeroDimConf>().max_solve_wall_clock_duration; limit > 0)
                    solve_deadline_ = std::chrono::steady_clock::now()
                        + std::chrono::duration_cast<std::chrono::steady_clock::duration>(std::chrono::duration<double>(limit));
 
                this->NotifyObservers(AlgorithmStarted<AnyZeroDim>(*this));
 
                // Speculative-full-path model: carry every path all the way through (pre-endgame +
                // endgame) as one unit, then detect crossings from the collected boundary points and
                // re-run any crossed path in full.  This is the same shape the distributed solve uses
                // (one worker per whole path), so serial and distributed share the per-path primitive
                // -- including computing the start point via the same ComputeStartPoint.
 
                // Compute every start point once, up front (deterministic per index, as the MPI
                // manager does).  The threaded dispatch ships each worker its start point, and the
                // serial path / crossed-path re-tracks reuse the identical points from here.
                std::vector<Vec<BaseComplexT>> start_points(num_start_points_);
                std::vector<SolnIndT>          all_indices(num_start_points_);
                for (decltype(num_start_points_) ii{0}; ii < num_start_points_; ++ii)
                {
                    auto idx = static_cast<SolnIndT>(ii);
                    start_points[ii] = ComputeStartPoint(idx);
                    all_indices[ii]  = idx;
                }
 
                // the records seam (ADR-0046): ensure the run is on record, recall any paths
                // already recorded (solve() is ensure-answered), and compute only the rest
                MaybeAttachAmbientRecords();
                std::vector<SolnIndT> indices_to_run = all_indices;
                if (records_)
                {
                    EnsureRunRecorded();
                    // which recorded outcomes are reused is RecallPolicy's business, applied inside
                    indices_to_run = RecallRecordedPaths(all_indices);
                }
 
                // num_threads: 0 = auto (every CPU this process may run on), 1 = serial, N = N
                // threads; BERTINI_NUM_THREADS overrides.  n_threads <= 1 takes the pool-free
                // serial path.
                const unsigned n_threads =
                    parallel::EffectiveThreadCount(this->template Get<ZeroDimConf>().num_threads);
 
                if (n_threads <= 1)
                {
                    for (auto idx : indices_to_run)
                    {
                        if (StopRequested() || SolveDeadlinePassed())
                            break;          // ExecuteOnePath would no-op; this just stops asking
                        ExecuteOnePath(MemberDuringEGContext(), idx, start_points[idx]);
                        // the serial path installs results directly (no StoreFullPathResult),
                        // so emit its record here
                        RecordCompletedPath(PackFullPathResult(idx));
                    }
                }
                else
                {
                    RunPathsThreaded(indices_to_run, start_points, n_threads);
                }
 
                TallyStoppedSolve(indices_to_run);
 
                // Crossed-path re-tracks run on the main thread against the (escalated) member
                // tracker.  They are typically rare (often none); running them in parallel too is a
                // deferred optimization.
                //
                // Skipped entirely when the solve was stopped: the midpath check compares paths
                // against each other at the endgame boundary, and a set where an arbitrary
                // number of paths never ran has nothing meaningful to compare.  Re-tracking
                // under a standing stop request would also just abandon each path again.
                if (!stopped_early_)
                    RunMidpathResolution([this, &start_points](SolnIndT idx){
                        ExecuteOnePath(MemberDuringEGContext(), idx, start_points[idx]);
                        RecordCompletedPath(PackFullPathResult(idx));   // a re-track appends a fresh record (last wins)
                    });
 
                PostEGAction();
 
                DeclareFiniteSolutionsResult();
 
                this->NotifyObservers(AlgorithmComplete<AnyZeroDim>(*this));
            }
 
 
 

            const auto& SolutionsUserCoords() const
            {
                if (!solutions_user_coords_fresh_)
                {
                    solutions_user_coords_.clear();
                    solutions_user_coords_.reserve(solutions_post_endgame_.size());
                    for (const auto& s : solutions_post_endgame_)
                    {
                        // A path that failed before or during the endgame leaves its endpoint
                        // slot default-constructed (zero coordinates); keep index alignment with
                        // SolutionMetadata (output filters on it) by emitting an empty
                        // placeholder rather than trying to dehomogenize an unset point.
                        if (s.size() == 0)
                            solutions_user_coords_.emplace_back();
                        else
                            solutions_user_coords_.push_back(this->TargetSystem().DehomogenizePoint(s));
                    }
                    solutions_user_coords_fresh_ = true;
                }
                return solutions_user_coords_;
            }

            const auto& SolutionsInternalCoords() const
            {
                return solutions_post_endgame_;
            }

            template<typename Pred>
            SolnCont<Vec<BaseComplexT>> SolutionsWhere(Pred pred, bool user_coords = true,
                                                       bool merge_multiplicities = false) const
            {
                auto const& sols = user_coords ? SolutionsUserCoords() : SolutionsInternalCoords();
                auto const& md   = SolutionMetadata();
                SolnCont<Vec<BaseComplexT>> out;
                for (size_t i = 0; i < md.size() && i < sols.size(); ++i)
                {
                    // merge_multiplicities: keep only the chosen representative of each multiplicity
                    // cluster, so a multiplicity-m solution appears once rather than m times.  A no-op for
                    // simple / at-infinity / failed endpoints (each is its own representative).  Issue #299.
                    if (merge_multiplicities && !md[i].multiplicity_representative)
                        continue;
                    if (pred(md[i]))
                        out.push_back(sols[i]);
                }
                return out;
            }

            SolnCont<Vec<BaseComplexT>> FiniteSolutions(bool user_coords = true, bool merge_multiplicities = false) const
            {
                return SolutionsWhere([](auto const& m){
                    return m.endgame_success_code == SuccessCode::Success && m.is_finite && !m.is_nonsolution; }, user_coords, merge_multiplicities);
            }

            SolnCont<Vec<BaseComplexT>> RealSolutions(bool user_coords = true, bool merge_multiplicities = false) const
            {
                return SolutionsWhere([](auto const& m){
                    return m.endgame_success_code == SuccessCode::Success && m.is_finite && !m.is_nonsolution && m.is_real; }, user_coords, merge_multiplicities);
            }

            SolnCont<Vec<BaseComplexT>> NonsingularSolutions(bool user_coords = true, bool merge_multiplicities = false) const
            {
                return SolutionsWhere([](auto const& m){
                    return m.endgame_success_code == SuccessCode::Success && m.is_finite && !m.is_nonsolution && !m.is_singular; }, user_coords, merge_multiplicities);
            }

            SolnCont<Vec<BaseComplexT>> SingularSolutions(bool user_coords = true, bool merge_multiplicities = false) const
            {
                return SolutionsWhere([](auto const& m){
                    return m.endgame_success_code == SuccessCode::Success && m.is_finite && !m.is_nonsolution && m.is_singular; }, user_coords, merge_multiplicities);
            }

            SolnCont<Vec<BaseComplexT>> Nonsolutions(bool user_coords = true) const
            {
                return SolutionsWhere([](auto const& m){
                    return m.endgame_success_code == SuccessCode::Success && m.is_nonsolution; }, user_coords);
            }

            SolnCont<Vec<BaseComplexT>> InfiniteSolutions(bool user_coords = true) const
            {
                return SolutionsWhere([](auto const& m){ return !m.is_finite; }, user_coords);
            }

            const auto& SolutionMetadata() const
            {
                return solution_final_metadata_;
            }

            NumErrorT DefaultPointMatchTolerance() const
            {
                return endgame_.FinalTolerance();
            }

            NumErrorT SamePointTolerance() const
            {
                return endgame_.FinalTolerance() *
                       this->template Get<PostProcessing>().same_point_tolerance_multiplier;
            }

            SolutionMetaDataT MetadataForPoint(Vec<BaseComplexT> const& point, NumErrorT tol,
                                               bool user_coords = true, bool representatives_only = true) const
            {
                auto const& sols = user_coords ? SolutionsUserCoords() : SolutionsInternalCoords();
                auto const& md   = SolutionMetadata();
 
                std::vector<size_t> matches;
                for (size_t i = 0; i < md.size() && i < sols.size(); ++i)
                {
                    if (representatives_only && !md[i].multiplicity_representative)
                        continue;
                    if (!bertini::IsDistinct(sols[i], point, tol))
                        matches.push_back(i);
                }
 
                if (matches.empty())
                    throw std::runtime_error("metadata_for: no solution matches the given point within tol");
 
                // More than one representative among the matches => tol grouped two distinct clusters.
                // With representatives_only (the default) and the default tol this cannot happen -- distinct
                // representatives are at least SamePointTolerance apart and the window is smaller -- so it
                // only fires for a deliberately coarse user-supplied tol, where it is the honest answer.
                std::vector<size_t> reps;
                for (auto i : matches)
                    if (md[i].multiplicity_representative)
                        reps.push_back(i);
                if (reps.size() > 1)
                    throw std::runtime_error("metadata_for: the point matches more than one distinct solution cluster; reduce tol");
                if (reps.size() == 1)
                    return md[reps.front()];
                // No representative matched (only reachable with representatives_only=false, when the point
                // is a non-representative copy whose representative fell outside tol); return the matched
                // copy's own record (it still carries the cluster's multiplicity).
                return md[matches.front()];
            }

            std::vector<SolutionMetaDataT> CoincidentMetadataForPoint(Vec<BaseComplexT> const& point, NumErrorT tol, bool user_coords = true) const
            {
                auto const& sols = user_coords ? SolutionsUserCoords() : SolutionsInternalCoords();
                auto const& md   = SolutionMetadata();
 
                // Find the representative NEAREST the point (within tol).  Nearest, not "all within tol", so
                // the cluster is identified unambiguously even for a coarse tol.
                size_t rep = md.size();
                NumErrorT best = std::numeric_limits<NumErrorT>::infinity();
                for (size_t i = 0; i < md.size() && i < sols.size(); ++i)
                {
                    if (!md[i].multiplicity_representative)
                        continue;
                    NumErrorT d = static_cast<NumErrorT>((sols[i] - point).template lpNorm<Eigen::Infinity>());
                    if (d <= tol && d < best) { best = d; rep = i; }
                }
                if (rep == md.size())
                    throw std::runtime_error("metadata_for: no solution matches the given point within tol");
 
                // Gather the whole cluster: every endpoint the solver deemed the same solution as the
                // representative (within the same-point clustering tolerance).
                const NumErrorT same_tol = SamePointTolerance();
                std::vector<SolutionMetaDataT> out;
                for (size_t i = 0; i < md.size() && i < sols.size(); ++i)
                    if (!bertini::IsDistinct(sols[i], sols[rep], same_tol))
                        out.push_back(md[i]);
                return out;
            }

            const auto& EndgameBoundarySolutions() const
            {
                return solutions_at_endgame_boundary_;
            }

            const MidpathCheckReport& EndgameBoundaryMetadata() const
            {
                return midpath_report_;
            }

            SolveReport Report() const
            {
                return SummarizeSolve(SolutionMetadata(), EndgameBoundaryMetadata());
            }
 
        private:

            void PreSolveChecks() const
            {
                if (num_start_points_ > solutions_at_endgame_boundary_.max_size())
                    throw std::runtime_error("start system has more solutions than container for results.  I refuse to continue until this has been addressed.");
            }
 
 
            void PreSolveSetup()
            {
                // The midpath checker is a consumer of the solver's MidPathConfig, not a second
                // home for it, so whatever the caller last set on the solver takes effect here --
                // whether they named the config or reached the field through the flat surface.
                midpath_.Set(this->template Get<MidPathConfig>());
 
                // Fixed-multiple precision: FixedPrecisionConfig.precision (on the tracker) is the
                // authoritative precision for the whole solve.  Lift the tracker, the ambient/thread
                // precision, the start-point precision (via initial_ambient_precision), and the systems
                // all to that one value -- TrackerLoopInitialization requires them to agree.  (Double is
                // fixed at 16 and adaptive manages its own precision, so neither enters here.)
                if constexpr (!tracking::TrackerTraits<TrackerType>::IsAdaptivePrec)
                {
                    auto const& fp = GetTracker().template Get<PrecisionConfig>();
                    // Double validates (its PrecisionSetup throws on any precision other than 16);
                    // fixed-multiple adopts the value into precision_.
                    GetTracker().PrecisionSetup(fp);
                    if constexpr (!std::is_same<BaseComplexT, complex_dbl>::value)
                        if (fp.precision != 0)
                        {
                            ZeroDimConf zdc = this->template Get<ZeroDimConf>();
                            zdc.initial_ambient_precision = fp.precision;  // -> thread + start-point precision
                            this->template Set<ZeroDimConf>(zdc);
                        }
                }
 
                auto num_as_size_t = static_cast<SolnIndT>(num_start_points_);
 
                solution_final_metadata_.resize(num_as_size_t);
                solutions_at_endgame_boundary_.resize(num_as_size_t);
                solutions_post_endgame_.resize(num_as_size_t);
 
                SetMidpathRetrackTol(this->template Get<Tolerances>().newton_before_endgame);
 
                // Default the start-point precision to the initial ambient precision.  A caller can
                // override it via SetStartPointPrecision (e.g. to carry a higher precision forward
                // from a previous solve whose output feeds this one).
                if (!start_point_precision_set_by_user_)
                    start_point_precision_ = this->template Get<ZeroDimConf>().initial_ambient_precision;
            }

            struct BeforeEGContext
            {
                TrackerType& tracker;
                tracking::PathPrecisionRecorder<TrackerType>& path_prec;
            };
            struct DuringEGContext
            {
                // const ref: the residual/dehomogenize/precision calls are const-callable, and the
                // RefToGiven policy hands back a const target system.  A mutable thread-owned clone
                // binds here too.
                SystemType const& target_sys;
                TrackerType& tracker;
                EndgameType& endgame;
                tracking::PathPrecisionRecorder<TrackerType>& path_prec;
            };
 
            BeforeEGContext MemberBeforeEGContext()
            {
                return BeforeEGContext{ GetTracker(), path_prec_ };
            }
            DuringEGContext MemberDuringEGContext()
            {
                return DuringEGContext{ TargetSystem(), GetTracker(), GetEndgame(), path_prec_ };
            }

            static void HarvestPathPrecision(SolutionMetaData<BaseComplexT>& smd,
                                             tracking::PathPrecisionRecorder<TrackerType> const& rec)
            {
                if (rec.DidPrecisionIncrease() && !smd.precision_changed)
                {
                    smd.precision_changed = true;
                    smd.time_of_first_prec_increase = rec.TimeOfIncrease();
                }
                using std::max;
                smd.max_precision_used = max(smd.max_precision_used, rec.MaxPrecision());
            }

            Vec<BaseComplexT> ComputeStartPoint(SolnIndT soln_ind)
            {
                SetThreadPrecision(start_point_precision_);
                static std::mutex start_point_mutex;
                std::lock_guard<std::mutex> lock(start_point_mutex);
                return StartSystem().template StartPoint<BaseComplexT>(soln_ind);
            }

            void ExecuteBeforeEG(BeforeEGContext ctx, SolnIndT soln_ind, Vec<BaseComplexT> const& start_point)
            {
                ReseedThisThread(static_cast<uint64_t>(soln_ind));
 
                auto& smd = solution_final_metadata_[soln_ind];
                smd.path_index    = soln_ind;
                smd.solution_index = soln_ind;
 
                // The step tally spans the whole path -- this track and every sub-track the
                // endgame will issue -- so it starts here and is read after each stage.
                ctx.tracker.ResetCumulativeStepCounts();
 
                // The wall-clock budget, likewise per path and spanning both stages: armed here as
                // a deadline (so the endgame's many sub-tracks share it), cleared when the path is
                // over.  Recorded with the path whether or not it bites.
                auto const wall_clock_limit = this->template Get<ZeroDimConf>().max_path_wall_clock_duration;
                smd.wall_clock_limit_seconds = wall_clock_limit;
                ctx.tracker.ClearMaxWallClockTime();
                if (wall_clock_limit > 0)
                    ctx.tracker.SetMaxWallClockDuration(std::chrono::duration<double>(wall_clock_limit));
                // the whole-solve deadline, when there is one and it comes first, is the one that binds
                if (solve_deadline_ && (!ctx.tracker.MaxWallClockTime() || *solve_deadline_ < *ctx.tracker.MaxWallClockTime()))
                    ctx.tracker.SetMaxWallClockTime(*solve_deadline_);
 
                auto initial_prec = this->template Get<ZeroDimConf>().initial_ambient_precision;
                // SetThreadPrecision: writes thread-local only, safe from concurrent threads.
                SetThreadPrecision(initial_prec);
 
                // Draw the condition-number probe direction ONCE here, for the whole track of this
                // point (pre-endgame AND endgame).  We have just reseeded this thread's RNG to a value
                // determined solely by the path index, so the probe is deterministic per path and
                // identical whether the path is tracked in a serial loop or on a distributed worker --
                // regardless of how the paths were split across ranks.  It is NOT refreshed for the
                // endgame (which would churn the RNG across its sample-circle sub-tracks), so the same
                // direction persists start to finish.
                ctx.tracker.RefreshConditionDirection();
 
                // The times are stored precision-free (mpq_rational); materialize them as the tracking
                // complex type at the current working precision.
                BaseComplexT t_start           ( BaseRealT(this->template Get<ZeroDimConf>().start_time) );
                BaseComplexT t_endgame_boundary( BaseRealT(this->template Get<ZeroDimConf>().endgame_boundary) );
 
                // The start point is supplied by the caller (computed once, authoritatively, via
                // ComputeStartPoint) rather than regenerated here.  In a distributed solve rank 0
                // computes it and sends it with the task, so a worker never derives its own -- that,
                // plus authoritative pi (issue #156), is what makes the start point identical across
                // serial and distributed runs.
 
                // Reset to the configured initial step size for this fresh path.  In the
                // speculative-full-path model a path's endgame runs immediately before the next
                // path's pre-endgame tracking on the same tracker, and the endgame leaves
                // ReinitializeInitialStepSize(false) with a tiny step; without restoring it here the
                // next path would crawl from the start time with that tiny step (effectively a hang).
                ctx.tracker.ReinitializeInitialStepSize(true);
 
                // Begin every AMP path in hardware double precision; the AMP step criteria escalate to
                // multiprecision only where a path actually needs it (and the endgame, which clears
                // this override below, manages its own precision).  Total-degree start points are born
                // at LowestMultiplePrecision (the generator's arithmetic widens past the requested
                // digits, see issue #308), and initial_ambient_precision is itself a multiprecision
                // value, so without forcing double here the tracker starts -- and, in practice, stays
                // -- in multiprecision on every path, even trivially well-conditioned ones.
                if constexpr (tracking::TrackerTraits<TrackerType>::IsAdaptivePrec)
                    ctx.tracker.SetStartPrecision(DoublePrecision());
 
                Vec<BaseComplexT> result;
                auto tracking_success = ctx.tracker.TrackPath(result, t_start, t_endgame_boundary, start_point);
 
                solutions_at_endgame_boundary_[soln_ind] =
                    EGBoundaryMetaDataT({ result, tracking_success, ctx.tracker.CurrentStepsize(), ctx.tracker.CurrentPrecision() });
 
                // Clear the start-precision override so it does not leak into the endgame, which
                // shares this tracker instance for its sample circles.
                if constexpr (tracking::TrackerTraits<TrackerType>::IsAdaptivePrec)
                    ctx.tracker.SetStartPrecision(std::nullopt);
 
                smd.pre_endgame_success_code = tracking_success;
 
                if (tracking::TrackerTraits<TrackerType>::IsAdaptivePrec)
                    HarvestPathPrecision(smd, ctx.path_prec);
 
                StampWhereThePathGotTo(smd, ctx.tracker, tracking_success, ctx.tracker.CurrentTime());
            }

            void StampWhereThePathGotTo(SolutionMetaData<BaseComplexT>& smd, TrackerType const& tracker,
                                        SuccessCode code, BaseComplexT const& time_reached) const
            {
                using std::max;
                smd.final_time_used      = time_reached;
                smd.num_successful_steps = tracker.CumulativeSuccessfulSteps();
                smd.num_failed_steps     = tracker.CumulativeFailedSteps();
                smd.max_precision_used   = max(smd.max_precision_used,
                                               static_cast<decltype(smd.max_precision_used)>(tracker.CurrentPrecision()));
                if (code == SuccessCode::Success)
                    smd.latest_path_point.resize(0);
                else
                    smd.latest_path_point = tracker.CurrentPoint();
            }

            struct PathThreadState
            {
                System          sys;         // homotopy, tracked by `tracker`
                System          target_sys;  // for function residuals / dehomogenization
                TrackerType     tracker;
                EndgameType     endgame;
                tracking::PathPrecisionRecorder<TrackerType> path_prec;
 
                DuringEGContext Context()
                {
                    return DuringEGContext{ target_sys, tracker, endgame, path_prec };
                }
            };

            auto MakeThreadStateFactory()
            {
                return [this]() -> std::unique_ptr<PathThreadState>
                {
                    // Trackers/Endgames have no default ctor, so aggregate-init via new; {} defaults
                    // the path's precision recorder.
                    std::unique_ptr<PathThreadState> s(
                        new PathThreadState{ Clone(GetTracker().GetSystem()), Clone(TargetSystem()),
                                             GetTracker(), GetEndgame(), {} });
                    s->tracker.SetSystem(s->sys);       // re-point the copy at its own System
                    s->endgame.SetTracker(s->tracker);  // ... and the endgame at that tracker
                    SetThreadPrecision(this->template Get<ZeroDimConf>().initial_ambient_precision);
                    return s;
                };
            }

            auto MakeThreadTrackFn()
            {
                return [this](std::unique_ptr<PathThreadState>& state,
                              parallel::StartPointTask<BaseComplexT> const& task)
                           -> parallel::FullPathResult<BaseComplexT>
                {
                    auto idx = static_cast<SolnIndT>(task.path_index);
                    ExecuteOnePath(state->Context(), idx, task.start_point);
                    return PackFullPathResult(idx);
                };
            }

            void RunPathsThreaded(std::vector<SolnIndT> const& indices,
                                  std::vector<Vec<BaseComplexT>> const& start_points,
                                  unsigned n_threads)
            {
                using Task   = parallel::StartPointTask<BaseComplexT>;
                using Result = parallel::FullPathResult<BaseComplexT>;
 
                auto state_factory = MakeThreadStateFactory();
                auto track_fn      = MakeThreadTrackFn();
 
                parallel::WorkerThreadPool<Task, Result, decltype(state_factory), decltype(track_fn)>
                    pool(static_cast<int>(n_threads), state_factory, track_fn);
 
                for (auto idx : indices)
                    pool.submit(Task{ static_cast<std::size_t>(idx), start_points[idx] });
 
                for (std::size_t k = 0; k < indices.size(); ++k)
                    StoreFullPathResult(pool.collect());
 
                pool.shutdown();
            }

            void EscalateRetrackSettings()
            {
                midpath_retrack_tolerance_ *= this->template Get<AutoRetrack>().midpath_decrease_tolerance_factor;
                GetTracker().SetTrackingTolerance(midpath_retrack_tolerance_);
 
                const auto default_predictor = tracking::predict::DefaultPredictor();
                if (tracking::predict::Order(GetTracker().GetPredictor())
                    < tracking::predict::Order(default_predictor))
                    GetTracker().SetPredictor(default_predictor);
            }
 
 

            void ExecuteOnePath(DuringEGContext ctx, SolnIndT soln_ind, Vec<BaseComplexT> const& start_point)
            {
                // Carry the executing tracker on the path events: &ctx.tracker is the member tracker
                // in a serial solve and the thread-local clone in a threaded one, so a meta-observer
                // attaches its per-path sub-observer to the tracker that actually runs this path.
                Observable const* exec_tracker = &ctx.tracker;
 
                // Somebody asked us to stop before this path began, so do not begin it.  Its
                // metadata keeps the default NeverStarted, which is the honest description and
                // is already distinguishable from a path we started and abandoned.  Checked
                // here rather than only in the dispatch loops so the threaded pool, which has
                // every task submitted up front, drains without doing work.
                if (StopRequested() || SolveDeadlinePassed())
                    return;
 
                // wall-clock the whole path (pre-endgame + endgame).  steady_clock: monotonic, safe on
                // the worker thread; written to this path's own metadata slot, so no cross-path race.
                auto const path_start_clock = std::chrono::steady_clock::now();
                auto stamp_path_time = [&]{
                    solution_final_metadata_[soln_ind].path_time_seconds =
                        std::chrono::duration<double>(std::chrono::steady_clock::now() - path_start_clock).count();
                };
 
                this->NotifyObservers(PathStarted<AnyZeroDim>(*this, static_cast<std::size_t>(soln_ind), exec_tracker));
 
                // One precision record for the whole path, across the track to the boundary and
                // every endgame sub-track: started here, read after each stage, detached at the end.
                constexpr bool adaptive = tracking::TrackerTraits<TrackerType>::IsAdaptivePrec;
                if (adaptive)
                {
                    ctx.path_prec.Reset();
                    ctx.tracker.AddObserver(ctx.path_prec);
                }
 
                ctx.tracker.SetTrackingTolerance(midpath_retrack_tolerance_);
                ExecuteBeforeEG(BeforeEGContext{ ctx.tracker, ctx.path_prec }, soln_ind, start_point);
 
                ReadSolveTimeoutAsAnInterrupt(solution_final_metadata_[soln_ind].pre_endgame_success_code, path_start_clock);
                if (solution_final_metadata_[soln_ind].pre_endgame_success_code != SuccessCode::Success)
                {
                    if (adaptive)
                        ctx.tracker.RemoveObserver(ctx.path_prec);
                    ctx.tracker.ClearMaxWallClockTime();
                    stamp_path_time();
                    this->NotifyObservers(PathComplete<AnyZeroDim>(*this, static_cast<std::size_t>(soln_ind), exec_tracker,
                                                                   solution_final_metadata_[soln_ind].pre_endgame_success_code));
                    return;
                }
 
                ctx.tracker.SetTrackingTolerance(this->template Get<Tolerances>().newton_during_endgame);
                ExecuteDuringEG(ctx, soln_ind);
                if (adaptive)
                    ctx.tracker.RemoveObserver(ctx.path_prec);
 
                ReadSolveTimeoutAsAnInterrupt(solution_final_metadata_[soln_ind].endgame_success_code, path_start_clock);
                ctx.tracker.ClearMaxWallClockTime();   // the budget was this path's; do not carry it to the next
                stamp_path_time();
                this->NotifyObservers(PathComplete<AnyZeroDim>(*this, static_cast<std::size_t>(soln_ind), exec_tracker,
                                                               solution_final_metadata_[soln_ind].endgame_success_code));
            }
 

            template<typename RedoOne>
            void RunMidpathResolution(RedoOne&& redo_one)
            {
                auto passed = midpath_.Check(solutions_at_endgame_boundary_, StartSystem());
 
                midpath_report_ = MidpathCheckReport{};
                midpath_report_.passed = passed;
                for (auto const& v : midpath_.GetCrossedPaths())
                    midpath_report_.crossed_path_indices.push_back(v.index());
                midpath_report_.num_crossings_detected =
                    static_cast<unsigned>(midpath_report_.crossed_path_indices.size());
 
                const auto max_attempts = this->template Get<ZeroDimConf>().max_num_crossed_path_resolve_attempts;
                unsigned num_resolve_attempts = 0;
                while (!passed && num_resolve_attempts < max_attempts)
                {
                    EscalateRetrackSettings();
                    for (auto const& v : midpath_.GetCrossedPaths())
                        if (v.rerun())
                            redo_one(static_cast<SolnIndT>(v.index()));
                    passed = midpath_.Check(solutions_at_endgame_boundary_, StartSystem());
                    ++num_resolve_attempts;
                }
 
                midpath_report_.num_resolve_attempts = num_resolve_attempts;
                midpath_report_.passed = passed;
                FlagUnresolvedCrossings(passed);
                RecordMidpathOutcome();
 
                if (!passed)
                    std::cerr << "warning: " << midpath_report_.num_crossings_detected
                              << " path crossing(s) detected at the endgame boundary remained unresolved "
                              << "after " << num_resolve_attempts << " re-track attempt(s); "
                              << "the affected solutions may be wrong.  Consider a higher-order predictor "
                              << "or a tighter tracking tolerance.  See EndgameBoundaryMetadata()." << std::endl;
            }
 

            void FlagUnresolvedCrossings(bool passed)
            {
                for (auto& smd : solution_final_metadata_)
                    smd.crossing_unresolved = false;
                if (passed)
                    return;
                for (auto const& v : midpath_.GetCrossedPaths())
                    solution_final_metadata_[static_cast<SolnIndT>(v.index())].crossing_unresolved = true;
            }

            void RecordMidpathOutcome()
            {
                if (!records_ || recalling_ || records_run_id_.empty())
                    return;
 
                boost::json::array crossed;
                for (auto const i : midpath_report_.crossed_path_indices)
                    crossed.push_back(boost::json::value(static_cast<std::int64_t>(i)));
                records_->AppendResult(records_run_id_,
                    {{"kind", "midpath"},
                     {"run", records_run_id_},
                     {"passed", midpath_report_.passed},
                     {"num_crossings_detected",
                      static_cast<std::int64_t>(midpath_report_.num_crossings_detected)},
                     {"num_resolve_attempts",
                      static_cast<std::int64_t>(midpath_report_.num_resolve_attempts)},
                     {"crossed_path_indices", std::move(crossed)}});
 
                for (SolnIndT idx = 0; idx < solution_final_metadata_.size(); ++idx)
                    if (solution_final_metadata_[idx].crossing_unresolved)
                        RecordCompletedPath(PackFullPathResult(idx));
            }

            void ExecuteDuringEG(DuringEGContext ctx, SolnIndT soln_ind)
            {
                ReseedThisThread(static_cast<uint64_t>(soln_ind) + static_cast<uint64_t>(num_start_points_));
 
                auto& smd = solution_final_metadata_[soln_ind];
 
                const auto& bdry_point = solutions_at_endgame_boundary_[soln_ind].path_point;
 
                ctx.tracker.SetStepSize(solutions_at_endgame_boundary_[soln_ind].last_used_stepsize);
                ctx.tracker.ReinitializeInitialStepSize(false);
 
                // Resume the endgame at the precision the path was actually using at the boundary,
                // rather than inferring it from the (always-multiprecision) point's mantissa, which
                // is unreliable for paths that tracked in double.
                auto start_prec = solutions_at_endgame_boundary_[soln_ind].precision;
                SetThreadPrecision(start_prec);
 
                ctx.endgame.SetBoundaryTime(BaseComplexT(BaseRealT(this->template Get<ZeroDimConf>().endgame_boundary)));
                ctx.endgame.SetTargetTime  (BaseComplexT(BaseRealT(this->template Get<ZeroDimConf>().target_time)));
 
                auto eg_success = ctx.endgame.Run(bdry_point);
 
                // The endgame writes its approximation only when it converges, so on a failure
                // the accessor still holds whatever the PREVIOUS path on this endgame left there.
                // A failed path's solution slot is therefore emptied, never copied: its stamp
                // (latest_path_point, final_time_used) says where it got to.
                if (eg_success == SuccessCode::Success)
                    solutions_post_endgame_[soln_ind] = ctx.endgame.template FinalApproximation<BaseComplexT>();
                else
                    solutions_post_endgame_[soln_ind].resize(0);
 
                smd.endgame_success_code = eg_success;
                StampWhereThePathGotTo(smd, ctx.tracker, eg_success,
                                       eg_success == SuccessCode::Success ? ctx.endgame.LatestTime()
                                                                          : ctx.tracker.CurrentTime());
 
                // Harvest the path's precision record REGARDLESS of the endgame outcome: the
                // precision a path used is a fact about the tracking that happened, and it is
                // exactly the FAILED paths -- e.g. slow divergers escalating in the mp lane --
                // whose precision honesty matters most for diagnostics.
                if (tracking::TrackerTraits<TrackerType>::IsAdaptivePrec)
                    HarvestPathPrecision(smd, ctx.path_prec);
 
                // an unsuccessful endgame has no final approximation, so the final-point-dependent
                // metadata cannot be computed.
                if (eg_success != SuccessCode::Success)
                    return;
                if (tracking::TrackerTraits<TrackerType>::IsAdaptivePrec)
                {
                    assert(Precision(solutions_post_endgame_[soln_ind])==Precision(ctx.endgame.template FinalApproximation<BaseComplexT>()));
                    SetThreadPrecision(Precision(solutions_post_endgame_[soln_ind]));
                    // the system no longer needs aligning by hand -- the Eval below does it
                    // from the point it is handed (#377)
                }
                smd.function_residual = static_cast<NumErrorT>(ctx.target_sys.Eval(solutions_post_endgame_[soln_ind]).template lpNorm<Eigen::Infinity>());
                smd.singular_values  = EndpointSingularValues(ctx.target_sys, solutions_post_endgame_[soln_ind]);
                smd.condition_number = SpectralConditionNumberOf(smd.singular_values);
                smd.newton_residual = ctx.tracker.LatestNormOfStep();
 
                smd.accuracy_estimate_internal_coords = ctx.endgame.ApproximateError();
                smd.accuracy_estimate_user_coords =
                    static_cast<NumErrorT>( (ctx.target_sys.DehomogenizePoint(solutions_post_endgame_[soln_ind]) -
                    ctx.target_sys.DehomogenizePoint(ctx.endgame.template PreviousApproximation<BaseComplexT>())).template lpNorm<Eigen::Infinity>() );
                smd.cycle_num = ctx.endgame.CycleNumber();
 
                // Honest precision/accuracy of this solution, as digit counts.  precision_digits is the
                // precision the endgame actually finished in -- which is exactly the precision of the
                // stored solution point (DoublePrecision() for a fast-lane path, the mpfr precision for an
                // escalated one).  accuracy_digits is how many of those digits the convergence agreement
                // supports: floor(-log10(accuracy_estimate_internal_coords)), clamped to [0, precision_digits].
                smd.precision_digits = static_cast<unsigned>(Precision(solutions_post_endgame_[soln_ind]));
                {
                    using std::log10; using std::floor; using std::min; using std::max;
                    double err = static_cast<double>(smd.accuracy_estimate_internal_coords);
                    unsigned acc = (err > 0.0)
                        ? static_cast<unsigned>(max(0.0, floor(-log10(err))))
                        : smd.precision_digits;
                    smd.accuracy_digits = min(acc, smd.precision_digits);
                }
            }
 
 
        protected:
            // virtual so ZeroDimSolver can append its extraneous-solution filter after the engine's
            // classification (it overrides this to call the base, then filter against the original
            // over-determined system).  Everything from here down -- the classification helpers, the
            // pack/store helpers, and the data members -- is protected so the derived algorithm can
            // read the per-endpoint metadata and endpoints it needs for that filter.
            virtual void PostEGAction()
            {
                ComputePostTrackMetadata();
            }
 
 

            void ComputePostTrackMetadata()
            {
                ClassifyFiniteAndReal();
                ComputeMultiplicities();
                ClassifySingular();
            }

            void ClassifyFiniteAndReal()
            {
                const auto& post = this->template Get<PostProcessing>();
 
                for (decltype(num_start_points_) ii{0}; ii < num_start_points_; ++ii)
                {
                    auto& smd = solution_final_metadata_[ii];
 
                    if (smd.endgame_success_code==SuccessCode::GoingToInfinity ||
                        smd.endgame_success_code==SuccessCode::SecurityMaxNormReached)
                    {
                        smd.is_finite = false; // the endgame already decided this path diverges
                        continue;
                    }
                    if (smd.endgame_success_code!=SuccessCode::Success)
                        continue;              // failed otherwise: leave defaults (not finite/real/singular)
 
                    // the system renders both verdicts: it is what knows which of its coordinates
                    // the questions are about (b2#403)
                    smd.is_finite = this->TargetSystem().IsFinite(solutions_post_endgame_[ii],
                                                                  post.endpoint_finite_threshold);
                    if (smd.is_finite)
                        smd.is_real = this->TargetSystem().IsReal(solutions_post_endgame_[ii],
                                                                  post.real_threshold);
                }
            }

            void ComputeMultiplicities()
            {
                const NumErrorT same_tol = this->SamePointTolerance();
 
                std::vector<Vec<BaseComplexT>> user_pts(num_start_points_);
                std::vector<char> eligible(num_start_points_, 0);
                for (decltype(num_start_points_) ii{0}; ii < num_start_points_; ++ii)
                {
                    auto const& smd = solution_final_metadata_[ii];
                    if (smd.endgame_success_code==SuccessCode::Success && smd.is_finite)
                    {
                        user_pts[ii] = this->TargetSystem().DehomogenizePoint(solutions_post_endgame_[ii]);
                        eligible[ii] = 1;
                    }
                }
 
                for (decltype(num_start_points_) ii{0}; ii < num_start_points_; ++ii)
                {
                    if (!eligible[ii]) continue;
                    for (decltype(num_start_points_) jj{ii+1}; jj < num_start_points_; ++jj)
                    {
                        if (!eligible[jj]) continue;
                        if ( static_cast<NumErrorT>((user_pts[ii] - user_pts[jj]).template lpNorm<Eigen::Infinity>()) < same_tol )
                        {
                            ++solution_final_metadata_[ii].multiplicity;
                            ++solution_final_metadata_[jj].multiplicity;
                            // jj coincides with the earlier ii, so it is a duplicate, not the
                            // representative; the lowest-index member of a cluster is never marked
                            // here, so it remains the single representative.
                            solution_final_metadata_[jj].multiplicity_representative = false;
                        }
                    }
                }
            }

            void ClassifySingular()
            {
                const NumErrorT cond_threshold = this->template Get<PostProcessing>().condition_number_threshold;
                for (decltype(num_start_points_) ii{0}; ii < num_start_points_; ++ii)
                {
                    auto& smd = solution_final_metadata_[ii];
                    if (smd.endgame_success_code!=SuccessCode::Success)
                        continue;
                    smd.is_singular = (smd.multiplicity > 1) || (smd.condition_number > cond_threshold);
                }
            }

            static NumErrorT EndpointSpectralConditionNumber(SystemType const& sys, Vec<BaseComplexT> const& endpoint)
            {
                return SpectralConditionNumberOf(EndpointSingularValues(sys, endpoint));
            }

            static Vec<NumErrorT> EndpointSingularValues(SystemType const& sys, Vec<BaseComplexT> const& endpoint)
            {
                const auto J = sys.Jacobian(endpoint);
                const auto sv = Eigen::JacobiSVD<Mat<BaseComplexT>>(J).singularValues();
                Vec<NumErrorT> out(sv.size());
                for (Eigen::Index i = 0; i < sv.size(); ++i)
                    out(i) = static_cast<NumErrorT>(sv(i));
                return out;
            }

            static NumErrorT SpectralConditionNumberOf(Vec<NumErrorT> const& singular_values)
            {
                if (singular_values.size() == 0)
                    return std::numeric_limits<NumErrorT>::infinity();
                const NumErrorT smallest = singular_values(singular_values.size()-1);
                if (smallest <= 0)
                    return std::numeric_limits<NumErrorT>::infinity();
                return singular_values(0) / smallest;
            }
 
 

        //  MPI pack/store helpers (only compiled when BERTINI2_HAVE_MPI is defined)
 
            // Pack everything computed for one whole path into a single result value.  Used as the
            // thread-safe handoff in BOTH the MPI protocol (serialized over the wire) and the
            // standalone threaded solve (a worker thread packs, the main thread installs via
            // StoreFullPathResult) -- the compute/install split is what keeps shared result arrays
            // race-free when several threads run distinct paths.
            parallel::FullPathResult<BaseComplexT> PackFullPathResult(SolnIndT idx) const
            {
                parallel::FullPathResult<BaseComplexT> r;
                r.path_index          = idx;
                r.pre_endgame_success_code = solutions_at_endgame_boundary_[idx].success_code;
                r.boundary_point      = solutions_at_endgame_boundary_[idx].path_point;
                r.boundary_stepsize   = solutions_at_endgame_boundary_[idx].last_used_stepsize;
                r.boundary_precision  = solutions_at_endgame_boundary_[idx].precision;
 
                auto const& smd = solution_final_metadata_[idx];
                r.endgame_success_code   = smd.endgame_success_code;
                r.solution    = solutions_post_endgame_[idx];
                r.function_residual = smd.function_residual;
                r.condition_number  = smd.condition_number;
                r.singular_values   = smd.singular_values;
                r.newton_residual   = smd.newton_residual;
                r.final_time_used   = smd.final_time_used;
                r.accuracy_estimate_internal_coords = smd.accuracy_estimate_internal_coords;
                r.accuracy_estimate_user_coords = smd.accuracy_estimate_user_coords;
                r.cycle_num         = smd.cycle_num;
                r.precision_digits = smd.precision_digits;
                r.accuracy_digits  = smd.accuracy_digits;
                r.precision_changed = smd.precision_changed;
                r.time_of_first_prec_increase = smd.time_of_first_prec_increase;
                r.max_precision_used = smd.max_precision_used;
                r.path_time_seconds = smd.path_time_seconds;
                r.num_successful_steps = smd.num_successful_steps;
                r.num_failed_steps     = smd.num_failed_steps;
                r.latest_path_point           = smd.latest_path_point;
                r.wall_clock_limit_seconds = smd.wall_clock_limit_seconds;
                // Always false while paths are being tracked -- FlagUnresolvedCrossings runs after
                // the last resolve round -- so a re-track never carries an earlier round's verdict
                // back out.  It is set here for the amending record that reports a flagged path
                // once the check is finally over.
                r.crossing_unresolved = smd.crossing_unresolved;
                return r;
            }
 
            // Install a whole-path result on the manager.  Overwrites cleanly, so re-dispatching a
            // crossed path in a later resolve round simply replaces its earlier (crossed) result.
            void StoreFullPathResult(parallel::FullPathResult<BaseComplexT> const& r)
            {
                auto idx = static_cast<SolnIndT>(r.path_index);
                solutions_at_endgame_boundary_[idx] =
                    EGBoundaryMetaDataT{ r.boundary_point, r.pre_endgame_success_code, r.boundary_stepsize, r.boundary_precision };
                solutions_post_endgame_[idx] = r.solution;
 
                auto& smd = solution_final_metadata_[idx];
                smd.path_index          = idx;
                smd.solution_index      = idx;
                smd.pre_endgame_success_code = r.pre_endgame_success_code;
                smd.endgame_success_code     = r.endgame_success_code;
                smd.function_residual   = r.function_residual;
                smd.condition_number    = r.condition_number;
                smd.singular_values     = r.singular_values;
                smd.newton_residual     = r.newton_residual;
                smd.final_time_used     = r.final_time_used;
                smd.accuracy_estimate_internal_coords = r.accuracy_estimate_internal_coords;
                smd.accuracy_estimate_user_coords = r.accuracy_estimate_user_coords;
                smd.cycle_num           = r.cycle_num;
                smd.precision_digits  = r.precision_digits;
                smd.accuracy_digits   = r.accuracy_digits;
                smd.precision_changed   = r.precision_changed;
                smd.time_of_first_prec_increase = r.time_of_first_prec_increase;
                smd.max_precision_used  = r.max_precision_used;
                smd.path_time_seconds   = r.path_time_seconds;
                smd.num_successful_steps = r.num_successful_steps;
                smd.num_failed_steps     = r.num_failed_steps;
                smd.latest_path_point           = r.latest_path_point;
                smd.wall_clock_limit_seconds = r.wall_clock_limit_seconds;
                // a recalled path brings its archived crossing verdict back with it; a freshly
                // tracked one carries false, which is what a path mid-solve is
                smd.crossing_unresolved = r.crossing_unresolved;
 
                // the records seam (ADR-0046): every topology installs completed paths here on
                // the main/manager thread, so emission is single-writer by construction
                RecordCompletedPath(r);
            }
 

        //
        //  the structured output directory seam (records; ADR-0046)
        //
 
        public:

            void RecordTo(std::shared_ptr<records::OutputDirectory> directory)
            {
                records_ = std::move(directory);
            }

            void RecordToPath(std::string const& path) override
            {
                RecordTo(records::OutputDirectory::Shared(path));
            }

            std::string RecordsRunIdentity() const override { return records_run_id_; }

            virtual SystemType const& RecordsTargetSystem() const { return TargetSystem(); }

            void AppendRecordJson(std::string const& record_json) override
            {
                if (records_)
                    records_->Append(boost::json::parse(record_json).as_object());
            }

            std::string PutRecordsDefinition(std::string const& content,
                                             std::string const& kind,
                                             std::string const& label = {}) override
            {
                return records_ ? records_->PutDefinition(content, kind, std::nullopt, label)
                                : std::string();
            }

            void DeclareFiniteSolutionsResult()
            {
                if (!records_ || records_run_id_.empty())
                    return;
                boost::json::array points;
                for (auto const& smd : solution_final_metadata_)
                    if (smd.is_finite)
                        points.push_back(boost::json::object{
                            {"run", records_run_id_},
                            {"index", static_cast<std::int64_t>(smd.path_index)}});
                boost::json::object result;
                result["kind"] = "result";
                result["name"] = "finite solutions [run " + records_run_id_ + "]";
                result["when"] = records::TimeStampNow();
                result["description"] = "auto-declared by the solver after post-processing";
                result["points"] = points;
                records_->Append(result);
                records_->RefreshIndex();
            }

            std::shared_ptr<records::OutputDirectory> const& Records() const { return records_; }

            unsigned long long NumPaths() const { return num_start_points_; }

            unsigned long long NumPathsRecalled() const { return num_recalled_; }

            bool WasStoppedEarly() const { return stopped_early_; }

            unsigned long long NumPathsNeverStarted() const { return num_never_started_; }

            std::string const& RecordsRunId() const { return records_run_id_; }

            void MaybeAttachAmbientRecords()
            {
                if (records_ || !parallel::IsManager())
                    return;
                if (auto const dir = records::AmbientRecordsPath())
                    records_ = records::OutputDirectory::Shared(*dir);
            }

            void SetRecordedStartProvenance(std::vector<boost::json::object> refs,
                                            std::string start_identity)
            {
                records_start_refs_ = std::move(refs);
                records_start_identity_ = std::move(start_identity);
            }

            std::string CanonicalSettingsText() const
            {
                std::ostringstream text;
                text << records::ConfigEncodingVersion << "\n"
                     << records::CanonicalEncoding(this->template Get<ZeroDimConf>()) << "\n"
                     << records::CanonicalEncoding(this->template Get<Tolerances>()) << "\n"
                     << records::CanonicalEncoding(this->template Get<AutoRetrack>()) << "\n"
                     << records::CanonicalEncoding(this->template Get<PostProcessing>()) << "\n"
                     << records::CanonicalEncoding(GetTracker().template Get<tracking::SteppingConfig>()) << "\n"
                     << records::CanonicalEncoding(GetTracker().template Get<tracking::NewtonConfig>()) << "\n"
                     << records::CanonicalEncoding(GetTracker().GetPredictor()) << "\n"
                     << records::CanonicalEncoding(GetTracker().template Get<PrecisionConfig>()) << "\n";
                AppendEndgameConfigEncodings(text,
                    typename endgame::AlgoTraits<EndgameType>::NeededConfigs{});
                text << records::CanonicalEncoding(this->template Get<MidPathConfig>()) << "\n";
                return text.str();
            }

            boost::json::object RecordsAsk() const
            {
                auto const settings = detail::Sha256(CanonicalSettingsText());
                boost::json::object ask;
                ask["op"] = "zerodim";
                ask["tracker"] = tracking::TrackerTraits<TrackerType>::kRecordName;
                ask["endgame"] = endgame::AlgoTraits<EndgameType>::kRecordName;
                ask["target"] = RecordsTargetSystem().ContentDigest().Hex();
                // the homotopy actually tracked: for seed-derived homotopies this is
                // redundant with (target, seed) but harmless; for USER homotopies (the
                // engine driven directly) it is the only thing telling asks apart
                ask["homotopy"] = Homotopy().ContentDigest().Hex();
                if (!records_start_identity_.empty())
                    ask["start"] = records_start_identity_;   // external start data is identity
                ask["config"] = settings.Hex();
                ask["seed"] = static_cast<std::int64_t>(GetGlobalSeed());
                return ask;
            }

            template <typename... EGConfTs>
            void AppendEndgameConfigEncodings(std::ostringstream& text,
                                              detail::TypeList<EGConfTs...>) const
            {
                ((text << records::CanonicalEncoding(
                    GetEndgame().template Get<EGConfTs>()) << "\n"), ...);
            }

            void EnsureRunRecorded()
            {
                auto const ask = RecordsAsk();
                records_run_id_ = detail::Sha256(boost::json::serialize(ask)).Hex().substr(0, 16);
 
                // the run's results file (the payload store) is created BEFORE anything
                // refers to it: line 1 is its self-description, so a wandering results
                // file says what run and ask it answers.  Idempotent on resume.
                records_->EnsureResultsFile(records_run_id_,
                    {{"kind", "results_header"},
                     {"schema", records::RecordSchemaVersion},
                     {"run", records_run_id_},
                     {"ask", ask}});
 
                records_run_resumed_ = false;
                for (auto const& rec : records_->Scan())
                    if (auto const* k = rec.if_contains("kind");
                        k && k->is_string() && k->get_string() == "run"
                        && rec.if_contains("run") && rec.at("run").is_string()
                        && rec.at("run").as_string() == records_run_id_)
                    {
                        // header already on record: this ask was asked before.  The
                        // recall event narrating THIS session lands after recall runs
                        // (RecallRecordedPaths knows the counts).
                        records_run_resumed_ = true;
                        return;
                    }
 
                // the archived system is a JSON document (like the configs): a
                // structured "system" parts view (variable groups, functions, ...) and
                // the CANONICAL ENCODING (b2sysenc) -- exact and complete: every block
                // (slices, randomization, blends), patch, and gamma survives.  Classic
                // syntax appears nowhere (an input/compat format, not an output one).
                // The encoding is the digest PREIMAGE, so the definition id EQUALS
                // target_digest and a copied-out file verifies itself (hash .encoding).
                auto const target_digest_hex = RecordsTargetSystem().ContentDigest().Hex();
                records_->PutDefinition(
                    records::SystemEncodingAsJson(RecordsTargetSystem().CanonicalEncodingText(),
                                                  target_digest_hex,
                                                  io::SystemPartsJson(RecordsTargetSystem())),
                    "systems", target_digest_hex);
                // the homotopy ACTUALLY TRACKED is archived too: its exact coefficients
                // (gamma, start-system constants, blend structure) are what the paths
                // followed, and for user homotopies (chained solves, the engine driven
                // directly) the archived encoding is the ONLY complete record -- no
                // seed can rebuild a homotopy constructed before this solve began.
                auto const homotopy_digest_hex = Homotopy().ContentDigest().Hex();
                records_->PutDefinition(
                    records::SystemEncodingAsJson(Homotopy().CanonicalEncodingText(),
                                                  homotopy_digest_hex,
                                                  io::SystemPartsJson(Homotopy())),
                    "systems", homotopy_digest_hex);
                // the settings, reconstructible: the SAME canonical text the digest is over
                std::string const config_text = CanonicalSettingsText();
                records_->PutDefinition(
                    records::ConfigTextAsJson(config_text, std::string(ask.at("config").as_string())),
                    "configs",
                    std::string(ask.at("config").as_string()));
                boost::json::object header;
                header["kind"] = "run";
                header["schema"] = records::RecordSchemaVersion;
                header["when"] = records::TimeStampNow();
                header["run"] = records_run_id_;
                header["op"] = "zerodim";
                // which software wrote this -- descriptive only, NEVER part of the ask
                // identity (a newer build answering the same ask must still recall)
                header["producer"] = records::ProducerInfo();
                header["ask"] = ask;
                // the field holds a DIGEST -- dereferencing it to the definitions/
                // file is the reader's job (its id equals this digest)
                header["target_digest"] = target_digest_hex;
                // the INTERNAL variable ordering, so views can label endpoint coordinates
                // truthfully (homogenized points have more coordinates than the user's
                // rendering declares)
                {
                    boost::json::array variable_names, user_variable_names;
                    auto const& homogenizers = TargetSystem().HomogenizingVariables();
                    for (auto const& v : TargetSystem().Variables())
                    {
                        variable_names.push_back(boost::json::value(v->name()));
                        bool const is_homogenizer = std::find(homogenizers.begin(),
                            homogenizers.end(), v) != homogenizers.end();
                        if (!is_homogenizer)
                            user_variable_names.push_back(boost::json::value(v->name()));
                    }
                    header["variables"] = variable_names;
                    // the USER's variables, in dehomogenized-point order: labels for
                    // the endpoint_user coordinates on track records
                    header["variables_user"] = user_variable_names;
                }
                // (the settings digest lives in ask.config; deref is the reader's job)
                header["num_paths"] = static_cast<std::int64_t>(num_start_points_);
                // where this run's computed paths live: history refers, results/ holds
                header["results_file"] = "results/" + records_run_id_.substr(0, 2)
                                         + "/" + records_run_id_ + ".jsonl";
                records_->Append(header);
            }

            static bool AbandonmentIsReusable(parallel::FullPathResult<BaseComplexT> const& r, double current_limit)
            {
                auto const is = [&](SuccessCode c){
                    return r.pre_endgame_success_code == c || r.endgame_success_code == c; };
                if (is(SuccessCode::ExternallyTerminated))
                    return false;
                if (is(SuccessCode::WallClockLimitReached))
                    return current_limit > 0 && r.wall_clock_limit_seconds > 0
                           && current_limit <= r.wall_clock_limit_seconds;
                return true;
            }

            std::vector<SolnIndT> RecallRecordedPaths(std::vector<SolnIndT> const& all_indices)
            {
                std::map<std::size_t, boost::json::object> recorded;   // last record per index wins
                for (auto const& rec : records_->ResultsOf(records_run_id_))
                {
                    auto const* k = rec.if_contains("kind");
                    if (!k || !k->is_string() || k->get_string() != "path")
                        continue;
                    if (auto const* idx = rec.if_contains("index"); idx && idx->is_int64())
                        recorded[static_cast<std::size_t>(idx->as_int64())] = rec;
                }
 
                auto const policy = this->template Get<RecordsConfig>().recall;
                auto const current_limit = this->template Get<ZeroDimConf>().max_path_wall_clock_duration;
 
                num_recalled_ = 0;
                recalling_ = true;
                std::vector<SolnIndT> missing;
                for (auto const idx : all_indices)
                {
                    auto const found = recorded.find(static_cast<std::size_t>(idx));
                    if (found == recorded.end() || policy == RecallPolicy::Nothing)
                    {
                        missing.push_back(idx);
                        continue;
                    }
 
                    auto decoded = records::DecodeFullPathResult<BaseComplexT>(found->second, idx);
 
                    // The rule for reusing a record: a path record is reusable when the run
                    // that would reuse it asks no more of that path than the run that produced
                    // it.  A completed path -- success, divergence, a failure under settings
                    // that are part of the ask -- always qualifies, since nothing more is being
                    // asked.  An ABANDONED path was ended for a reason deliberately kept out of
                    // the ask, so the record says "we gave up here", which is true and is not
                    // an answer; whether to reuse it is the user's call (RecallPolicy).  Under
                    // the default: a path somebody stopped is always re-tracked, since an
                    // interrupt has no budget to compare; one a wall-clock limit cut off is
                    // reused only while the current budget asks no more patience than the one
                    // that abandoned it, and re-tracked once the budget goes up -- which is
                    // exactly when paying again is the point.  Either way its record stays, so
                    // the store still says what happened; it is simply not mistaken for a
                    // conclusion.
                    if (policy == RecallPolicy::Completed && !AbandonmentIsReusable(decoded, current_limit))
                    {
                        missing.push_back(idx);
                        continue;
                    }
 
                    StoreFullPathResult(std::move(decoded));
                    ++num_recalled_;
                }
                recalling_ = false;
 
                // the narrative stays complete: a recalled ask WAS asked.  One small
                // history line per re-ask -- when, how much came from the store, how
                // much this session still computed (the kill-and-rerun story, in the
                // record).  Fresh runs (no prior header) narrate via their run header.
                if (records_run_resumed_)
                    records_->Append({{"kind", "recall"},
                                      {"run", records_run_id_},
                                      {"when", records::TimeStampNow()},
                                      {"num_recalled", static_cast<std::int64_t>(num_recalled_)},
                                      {"num_computed", static_cast<std::int64_t>(missing.size())},
                                      {"producer", records::ProducerInfo()}});
                return missing;
            }

            void TallyStoppedSolve(std::vector<SolnIndT> const& dispatched)
            {
                num_never_started_ = 0;
                bool abandoned_any = false;
 
                for (auto idx : dispatched)
                {
                    auto const code = solution_final_metadata_[idx].pre_endgame_success_code;
                    if (code == SuccessCode::NeverStarted)
                        ++num_never_started_;
                    else if (code == SuccessCode::ExternallyTerminated)
                        abandoned_any = true;
 
                    if (solution_final_metadata_[idx].endgame_success_code
                            == SuccessCode::ExternallyTerminated)
                        abandoned_any = true;
                }
 
                stopped_early_ = abandoned_any || num_never_started_ > 0;
            }

            void RecordCompletedPath(parallel::FullPathResult<BaseComplexT> const& r)
            {
                if (!records_ || recalling_ || records_run_id_.empty())
                    return;
 
                // A path we never began has nothing to report.  This happens when a solve is
                // stopped partway: every path after the stop is left NeverStarted.  Their
                // absence is explained once, at the run level, by the record that says the
                // solve was cut short -- which is better than thousands of identical rows
                // saying we did not get to it.  A path we DID begin and then abandoned is a
                // different thing entirely and is recorded, with how far it got.
                if (r.pre_endgame_success_code == SuccessCode::NeverStarted)
                    return;
                auto record = records::EncodeFullPathResult(r);
                record["kind"] = "path";
                record["run"] = records_run_id_;
                record["index"] = static_cast<std::int64_t>(r.path_index);
                // success / diverged / failed -- a truncation near infinity is a verdict,
                // not a failure (the SummarizeSolve bucketing, PR #46, now in the records)
                record["status"] = records::CoarsePathStatus(r.endgame_success_code);
                // the endpoint in USER coordinates too: audits read results in the
                // variables the user wrote, not the internal homogenized ones
                if (r.solution.size() != 0)
                    record["endpoint_user"] = records::EncodePoint(
                        this->TargetSystem().DehomogenizePoint(r.solution));
                // a chained/given start carries its true provenance; otherwise the
                // canonical start-system label (provenance bottoms out here)
                record["start"] = (r.path_index < records_start_refs_.size())
                    ? records_start_refs_[r.path_index]
                    : boost::json::object{{"kind", "start_label"},
                                          {"index", static_cast<std::int64_t>(r.path_index)}};
                records_->AppendResult(records_run_id_, record);
            }
 
 
        protected:
 
        protected:

        //  private data members
 
            std::shared_ptr<records::OutputDirectory> records_;  
            std::string records_run_id_;      
            std::vector<boost::json::object> records_start_refs_;  
            std::string records_start_identity_;  
            bool recalling_ = false;          
            bool records_run_resumed_ = false; 
            unsigned long long num_recalled_ = 0;  
 
            bool stopped_early_ = false;          
            std::optional<std::chrono::steady_clock::time_point> solve_deadline_;  

            bool SolveDeadlinePassed() const
            {
                return solve_deadline_ && std::chrono::steady_clock::now() >= *solve_deadline_;
            }

            void ReadSolveTimeoutAsAnInterrupt(SuccessCode& code, std::chrono::steady_clock::time_point path_start) const
            {
                if (code != SuccessCode::WallClockLimitReached || !SolveDeadlinePassed())
                    return;
                auto const path_limit = this->template Get<ZeroDimConf>().max_path_wall_clock_duration;
                if (path_limit > 0)
                {
                    auto const path_deadline = path_start
                        + std::chrono::duration_cast<std::chrono::steady_clock::duration>(std::chrono::duration<double>(path_limit));
                    if (path_deadline <= *solve_deadline_)
                        return;     // the path's own budget ran out first: that is what stopped it
                }
                code = SuccessCode::ExternallyTerminated;
            }
            unsigned long long num_never_started_ = 0;  
 
            unsigned long long num_start_points_;  
            NumErrorT midpath_retrack_tolerance_;  
            MidpathCheckReport midpath_report_; 
            unsigned start_point_precision_ = DoublePrecision(); 
            bool start_point_precision_set_by_user_ = false;  
 

            // i feel like these should be factored out into some policy class which prescribes how they are used, so that the actions taken are customizable.
            tracking::PathPrecisionRecorder<TrackerType> path_prec_;   
 

            TrackerType tracker_;    
            EndgameType endgame_;    
            MidpathType midpath_;    
 
 

            SolnCont< EGBoundaryMetaDataT > solutions_at_endgame_boundary_; 
            SolnCont<Vec<BaseComplexT> > solutions_post_endgame_;          
            mutable SolnCont<Vec<BaseComplexT> > solutions_user_coords_; 
            mutable bool solutions_user_coords_fresh_ = false;            
            SolnCont<SolutionMetaDataT> solution_final_metadata_;        
 
 
        }; // struct HomotopySolver
 
 
        // impl namespace deliberately NOT named `detail`: bertini::detail (TypeList, Configured)
        // is referenced unqualified throughout this file, and a bertini::algorithm::detail would
        // shadow it for any code defined after this point (e.g. ApplyParsedConfigs).
        namespace zero_dim_detail {

        template<typename SystemType>
        struct OwnedHomotopy
        {
            using StartSystemBaseT = bertini::start_system::StartSystem;                  
            using FactoryT = bertini::start_system::StartSystemFactory<SystemType>;       

            OwnedHomotopy(SystemType const& target, FactoryT factory, std::string const& path_variable_name)
             : user_target_(Clone(target)), owned_target_(Clone(target)), owned_factory_(std::move(factory))
            {
                ConsistencyCheck();              // feasibility (no path var; not under-constrained; polynomial)
                SquareUp();                      // randomize an over-determined system down to square
                RankCheck();                     // square, but can isolated solutions even exist?
                PrepareTarget(owned_target_);    // homogenize + auto-patch the (now square) target
                owned_start_    = owned_factory_(owned_target_);                         // start system over the prepared target
                // Treat the configured name as a BASE: mangle it only if it would collide with a
                // user variable, so the injected path variable can never overlap a user symbol.
                owned_homotopy_ = MakeHomotopy(owned_target_, *owned_start_,
                                               UniquePathVariableName(owned_target_, path_variable_name));
            }

            SystemType const&       BuiltTarget()   const { return owned_target_;   }
            SystemType const&       UserTarget()    const { return user_target_;    }
            StartSystemBaseT const& BuiltStart()    const { return *owned_start_;   }
            SystemType const&       BuiltHomotopy() const { return owned_homotopy_; }
 
        protected:
            SystemType                        user_target_;     
            SystemType                        owned_target_;    
            std::shared_ptr<StartSystemBaseT> owned_start_;     
            SystemType                        owned_homotopy_;  
            FactoryT                          owned_factory_;   
 
            // When the user's system was over-determined, SquareUp randomizes it down to square for
            // tracking and keeps the ORIGINAL (natural, un-homogenized) system here so ZeroDimSolver
            // can discard the extraneous solutions the squaring introduces.  was_randomized_ gates
            // that filter; randomization_matrix_ records the (exact) coefficient matrix used.
            bool                  was_randomized_ = false;       
            SystemType            original_natural_target_;      
            Mat<complex_mp>       randomization_matrix_;         

            void ConsistencyCheck() const
            {
                if (owned_target_.HavePathVariable())
                    throw std::runtime_error("unable to perform zero dim solve on target system -- has path variable, use a HomotopySolver instead.");
 
                // A square zero-dim system needs one equation per dimension.  Each projective
                // (homogeneous) variable group of size k spans P^{k-1}: its k coordinates carry only
                // k-1 dimensions because scale is free, so it needs one fewer equation than it has
                // variables.  Subtract that free scale per projective group before comparing.  An
                // UNDER-determined system has a positive-dimensional solution set, so no isolated
                // solutions to compute -- raise a helpful error rather than tracking garbage.
                if (owned_target_.NumVariables() - owned_target_.NumHomVariableGroups() > owned_target_.NumTotalFunctions())
                    throw std::runtime_error("unable to perform zero dim solve on target system -- it is under-determined (fewer equations than variables), so its solution set is positive-dimensional, not zero-dimensional.  ZeroDimSolver computes isolated solutions only; add equations, or use a positive-dimensional method.");
 
                if (!owned_target_.IsPolynomial())
                    throw std::runtime_error("unable to perform zero dim solve on target system -- system is non-polynomial, use a HomotopySolver instead.");
            }

            void SquareUp()
            {
                auto const dimension = owned_target_.NumVariables() - owned_target_.NumHomVariableGroups();
                if (owned_target_.NumTotalFunctions() <= dimension)
                    return; // already square (or under-determined, which ConsistencyCheck already rejected)
 
                original_natural_target_ = Clone(owned_target_);   // keep the original N-function system
                owned_target_            = owned_target_.Randomize();
                randomization_matrix_    = owned_target_.RandomizationMatrix();
                was_randomized_          = true;
            }

            void RankCheck() const
            {
                if (owned_target_.HavePathVariable())
                    return;
                auto const n = static_cast<Eigen::Index>(owned_target_.NumVariables());
                if (static_cast<Eigen::Index>(owned_target_.NumTotalFunctions()) != n || n == 0)
                    return; // not a square map: ConsistencyCheck governs feasibility here
 
                // a generic complex sample point: a zero-dimensional variety misses it almost surely,
                // where the Jacobian attains its generic (full) rank.
                Vec<complex_dbl> pt = Vec<complex_dbl>::Random(n);
                Mat<complex_dbl> J  = owned_target_.template Jacobian<complex_dbl>(pt);
                Eigen::FullPivLU<Mat<complex_dbl>> lu(J);
                lu.setThreshold(1e-10);
                if (lu.rank() < n)
                    throw std::runtime_error("unable to perform zero dim solve on target system -- the Jacobian is rank-deficient at a generic point, so the solution set is positive-dimensional, not zero-dimensional.  ZeroDimSolver computes isolated solutions only.");
            }

            static void PrepareTarget(SystemType& target)
            {
                target.Homogenize(); // work over projective coordinates
                target.AutoPatch();  // then patch if needed
            }
        };
 
        } // ns zero_dim_detail
 

        template<typename TrackerType, typename EndgameType, typename SystemType>
        struct ZeroDimSolver :
            private zero_dim_detail::OwnedHomotopy<SystemType>,
            public  HomotopySolver<TrackerType, EndgameType, SystemType>
        {
            using OwnedT   = zero_dim_detail::OwnedHomotopy<SystemType>;            
            using EngineT  = HomotopySolver<TrackerType, EndgameType, SystemType>;  
            using FactoryT = typename OwnedT::FactoryT;  

            ZeroDimSolver(SystemType const& target,
                          FactoryT factory = bertini::start_system::MakeStartFactory<bertini::start_system::TotalDegreeBinomial, SystemType>())
             : OwnedT(target, std::move(factory), ZeroDimConfig{}.path_variable_name),
               EngineT(OwnedT::BuiltTarget(), OwnedT::BuiltStart(), OwnedT::BuiltHomotopy())
            {}

            SystemType const& RecordsTargetSystem() const override { return OwnedT::UserTarget(); }

            bool WasRandomized() const { return this->was_randomized_; }
            Mat<complex_mp> const& RandomizationMatrix() const { return this->randomization_matrix_; }
 
        protected:
            void PostEGAction() override
            {
                EngineT::PostEGAction();   // the engine's finite/real/multiplicity/singular classification
 
                if (!this->was_randomized_)
                    return;
 
                // genuine roots satisfy the original system to ~solve accuracy; spurious ones miss it
                // by O(1).  A generous multiple of the endgame tolerance separates the two cleanly.
                const double threshold =
                    1e3 * static_cast<double>(this->template Get<TolerancesConfig>().newton_during_endgame);
 
                for (decltype(this->num_start_points_) ii{0}; ii < this->num_start_points_; ++ii)
                {
                    auto& smd = this->solution_final_metadata_[ii];
                    if (smd.endgame_success_code != SuccessCode::Success || !smd.is_finite)
                        continue;
 
                    auto user_pt = this->TargetSystem().DehomogenizePoint(this->solutions_post_endgame_[ii]);
                    // Evaluate the original system in DOUBLE precision: the filter only needs to tell a
                    // genuine root (tiny residual) from an extraneous one (O(1)), and a double residual
                    // does that robustly without a precision mismatch between the (possibly AMP/high
                    // precision) solution point and the original system's working precision.
                    Vec<complex_dbl> pt_d(user_pt.size());
                    for (Eigen::Index k = 0; k < user_pt.size(); ++k)
                        pt_d(k) = complex_dbl(user_pt(k));
                    auto residual = static_cast<NumErrorT>(
                        this->original_natural_target_.template Eval<complex_dbl>(pt_d).template lpNorm<Eigen::Infinity>());
                    smd.function_residual = residual;            // report residual against the ORIGINAL system
                    if (static_cast<double>(residual) > threshold)
                        smd.is_nonsolution = true;               // extraneous: finite, but not a solution of the original system
                }
            }
 
#ifdef BERTINI2_HAVE_MPI
            // Broadcast rank 0's authoritative owned systems to every rank, then re-seat the engine's
            // references and let the caller (RunParallel) re-point the tracker at the homotopy.
            void DistributeSystems(MPI_Comm comm) override
            {
                parallel::mpi_broadcast_serialized(comm, this->owned_target_,   0);
                // the OWNING shared_ptr<StartSystem> carries the concrete derived type polymorphically
                // (BOOST_CLASS_EXPORT); a base reference would serialize only the System slice.
                parallel::mpi_broadcast_serialized(comm, this->owned_start_,    0);
                parallel::mpi_broadcast_serialized(comm, this->owned_homotopy_, 0);
                this->ResetSystems(this->owned_homotopy_, *this->owned_start_, this->owned_target_);
            }
#endif
        };
 
    } // ns algo
 
} // ns bertini
 
// Include output formatters after ZeroDim is fully defined.
// output.hpp includes zero_dim_solve.hpp, so #pragma once prevents re-inclusion
// and the circular dependency is resolved.  Any TU that gets zero_dim_solve.hpp
// therefore also gets output.hpp, making WriteMainData/WriteRawData instantiable.
#include "bertini2/nag_algorithms/output.hpp"
 
namespace bertini {
namespace algorithm {
 
template<typename TrackerType, typename EndgameType, typename SystemType>
inline void
HomotopySolver<TrackerType,EndgameType,SystemType>::WriteMainData(std::ostream& out) const
{
    output::Classic<HomotopySolver>::MainData(out, *this);
}
 
template<typename TrackerType, typename EndgameType, typename SystemType>
inline void
HomotopySolver<TrackerType,EndgameType,SystemType>::WriteRawData(std::ostream& out) const
{
    output::Classic<HomotopySolver>::RawData(out, *this);
}
 
template<typename TrackerType, typename EndgameType, typename SystemType>
inline void
HomotopySolver<TrackerType,EndgameType,SystemType>::WriteFiniteSolutions(std::ostream& out) const
{
    output::Classic<HomotopySolver>::FiniteSolutions(out, *this);
}
 
template<typename TrackerType, typename EndgameType, typename SystemType>
inline void
HomotopySolver<TrackerType,EndgameType,SystemType>::WriteRealFiniteSolutions(std::ostream& out) const
{
    output::Classic<HomotopySolver>::RealFiniteSolutions(out, *this);
}
 
template<typename TrackerType, typename EndgameType, typename SystemType>
inline void
HomotopySolver<TrackerType,EndgameType,SystemType>::WriteNonsingularSolutions(std::ostream& out) const
{
    output::Classic<HomotopySolver>::NonsingularSolutions(out, *this);
}
 
template<typename TrackerType, typename EndgameType, typename SystemType>
inline void
HomotopySolver<TrackerType,EndgameType,SystemType>::WriteSingularSolutions(std::ostream& out) const
{
    output::Classic<HomotopySolver>::SingularSolutions(out, *this);
}
 
template<typename TrackerType, typename EndgameType, typename SystemType>
inline void
HomotopySolver<TrackerType,EndgameType,SystemType>::WriteRawSolutions(std::ostream& out) const
{
    output::Classic<HomotopySolver>::RawSolutions(out, *this);
}
 
} // ns algorithm
} // ns bertini
 
// Include config parsers after ZeroDim is fully defined, then provide the
// out-of-line definition of ApplyParsedConfigs.  Parsing headers do not include
// zero_dim_solve.hpp, so there is no circular dependency here.
#include "bertini2/io/parsing/settings_parsers.hpp"
 
namespace bertini {
namespace algorithm {

// Injects every element of a std::tuple<Ts...> into `target` via Set<T>.
// Works for any Configured<>-derived target whose typelist contains all Ts.
template<typename Target, typename... Ts>
void InjectParsedTuple(Target& target, std::tuple<Ts...> const& t) {
    (target.template Set<Ts>(std::get<Ts>(t)), ...);
}
 
template<typename TrackerType, typename EndgameType, typename SystemType>
void
HomotopySolver<TrackerType,EndgameType,SystemType>
    ::ApplyParsedConfigs(std::string const& config_str)
{
    using namespace parsing::classic;
 
    // 1. solver-owned configs (Tolerances, PostProcessing, ZeroDimConf, AutoRetrack)
    using ZDConfs = typename Config::UsedConfigs;
    auto zd = ConfigParser<ZDConfs>::Parse(config_str);
    InjectParsedTuple(*this, zd);
    // The path variable name comes from ZeroDimConf, so the homotopy must be
    // re-formed if the parsed name differs from the one used at construction.
    // Only re-run setup in that case: re-running unconditionally re-randomizes
    // gamma and the start system for no reason, and (before Homogenize() was made
    // idempotent) re-corrupted the prepared target system.
    // Compare against the collision-free name actually used at construction (the configured
    // name is a base that UniquePathVariableName may have mangled), so an unchanged config
    // does not spuriously look like a rename and force a rebuild.
    auto const effective_path_variable_name =
        UniquePathVariableName(this->TargetSystem(), this->template Get<ZeroDimConf>().path_variable_name);
    if (!Homotopy().HavePathVariable()
        || Homotopy().GetPathVariable()->name() != effective_path_variable_name)
    {
        DefaultSystemSetup();
    }
 
    // 2. Tracker — the classic file's own tracker settings, then the phase tolerance the solver
    //    drives.  Setup() still takes stepping and Newton positionally.
    using TkConfs = detail::TypeList<
        tracking::SteppingConfig,
        tracking::NewtonConfig,
        tracking::TrackerConfig>;
    auto tk = ConfigParser<TkConfs>::Parse(config_str);
    tracker_.template Set<tracking::TrackerConfig>(std::get<tracking::TrackerConfig>(tk));
    tracker_.Setup(
        std::get<tracking::TrackerConfig>(tk).predictor,
        this->template Get<Tolerances>().newton_before_endgame,
        std::get<tracking::SteppingConfig>(tk),
        std::get<tracking::NewtonConfig>(tk));
    tracker_.PrecisionSetup(PrecisionConfig(Homotopy()));
    endgame_.SetTracker(tracker_); // keep endgame's tracker ref consistent after reconfigure
 
    // 3. Endgame-owned configs — endgame_ is the concrete EndgameType deriving from
    //    Configured<its_configs>, so Set<T> is available directly.
    //    AlgoTraits<EndgameType>::NeededConfigs is public (unlike EndgameType::Configs).
    using EGConfs = typename endgame::AlgoTraits<EndgameType>::NeededConfigs;
    auto eg = ConfigParser<EGConfs>::Parse(config_str);
    InjectParsedTuple(endgame_, eg);
 
    // 4. Midpath config
    SetMidpath(ConfigParser<MidPathConfig>::Parse(config_str));
}
 
} // ns algorithm
} // ns bertini
 
 
// Explicit instantiation declarations — suppress re-instantiation of the production solver types in
// every including TU.  Six ZeroDimSolver combos (Tracker x Endgame) cover EVERY clone-owned start
// system (the start system is held polymorphically), plus six HomotopySolver combos for user
// homotopies.  Definitions in core/src/eti/zero_dim_eti.cpp + zero_dim_blackbox_eti.cpp; see ADR-0014.
#include "bertini2/endgames.hpp"
#include "bertini2/system/start_systems.hpp"
 
namespace bertini{ namespace algorithm{
 
extern template struct ZeroDimSolver<tracking::DoublePrecisionTracker,   typename endgame::EndgameSelector<tracking::DoublePrecisionTracker>::PSEG,     System>;
extern template struct ZeroDimSolver<tracking::DoublePrecisionTracker,   typename endgame::EndgameSelector<tracking::DoublePrecisionTracker>::Cauchy,   System>;
extern template struct ZeroDimSolver<tracking::MultiplePrecisionTracker, typename endgame::EndgameSelector<tracking::MultiplePrecisionTracker>::PSEG,   System>;
extern template struct ZeroDimSolver<tracking::MultiplePrecisionTracker, typename endgame::EndgameSelector<tracking::MultiplePrecisionTracker>::Cauchy, System>;
extern template struct ZeroDimSolver<tracking::AMPTracker,               typename endgame::EndgameSelector<tracking::AMPTracker>::PSEG,                 System>;
extern template struct ZeroDimSolver<tracking::AMPTracker,               typename endgame::EndgameSelector<tracking::AMPTracker>::Cauchy,               System>;
 
// user homotopies: the engine over caller-owned systems.  definitions in zero_dim_blackbox_eti.cpp
extern template struct HomotopySolver<tracking::DoublePrecisionTracker,   typename endgame::EndgameSelector<tracking::DoublePrecisionTracker>::PSEG,     System>;
extern template struct HomotopySolver<tracking::DoublePrecisionTracker,   typename endgame::EndgameSelector<tracking::DoublePrecisionTracker>::Cauchy,   System>;
extern template struct HomotopySolver<tracking::MultiplePrecisionTracker, typename endgame::EndgameSelector<tracking::MultiplePrecisionTracker>::PSEG,   System>;
extern template struct HomotopySolver<tracking::MultiplePrecisionTracker, typename endgame::EndgameSelector<tracking::MultiplePrecisionTracker>::Cauchy, System>;
extern template struct HomotopySolver<tracking::AMPTracker,               typename endgame::EndgameSelector<tracking::AMPTracker>::PSEG,                 System>;
extern template struct HomotopySolver<tracking::AMPTracker,               typename endgame::EndgameSelector<tracking::AMPTracker>::Cauchy,               System>;
 
}} // namespaces