solver_recording.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
//solver_recording.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.
//
//solver_recording.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 solver_recording.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 <cstdio>
#include <cstdlib>
#include <string>
 
#include <boost/json.hpp>
 
#include "bertini2/num_traits.hpp"
#include "bertini2/parallel/path_result.hpp"
 
namespace bertini {
namespace records {

inline std::string ExactDoubleText(double v)
{
    char buffer[32];
    std::snprintf(buffer, sizeof(buffer), "%.17g", v);
    return buffer;
}

inline double DoubleFromText(std::string const& s)
{
    return std::strtod(s.c_str(), nullptr);
}

inline std::string EncodeRealScalar(double v) { return ExactDoubleText(v); }
inline std::string EncodeRealScalar(real_mp const& v) { return v.str(0, std::ios::scientific); }

inline boost::json::array EncodeComplexScalar(complex_dbl const& z)
{
    return {ExactDoubleText(z.real()), ExactDoubleText(z.imag()),
            static_cast<std::int64_t>(DoublePrecision())};
}

inline boost::json::array EncodeComplexScalar(complex_mp const& z)
{
    return {z.real().str(0, std::ios::scientific), z.imag().str(0, std::ios::scientific),
            static_cast<std::int64_t>(z.precision())};
}

template <typename ComplexT>
ComplexT DecodeComplexScalar(boost::json::array const& triple);

template <>
inline complex_dbl DecodeComplexScalar<complex_dbl>(boost::json::array const& triple)
{
    return complex_dbl(DoubleFromText(std::string(triple.at(0).as_string())),
                       DoubleFromText(std::string(triple.at(1).as_string())));
}

template <>
inline complex_mp DecodeComplexScalar<complex_mp>(boost::json::array const& triple)
{
    auto const digits = static_cast<unsigned>(triple.at(2).as_int64());
    complex_mp z(real_mp(std::string(triple.at(0).as_string()), digits),
                 real_mp(std::string(triple.at(1).as_string()), digits));
    z.precision(digits);
    return z;
}

template <typename ComplexT>
boost::json::array EncodePoint(Vec<ComplexT> const& point)
{
    boost::json::array out;
    for (Eigen::Index ii = 0; ii < point.size(); ++ii)
        out.push_back(EncodeComplexScalar(point(ii)));
    return out;
}

template <typename ComplexT>
Vec<ComplexT> DecodePoint(boost::json::array const& coords)
{
    Vec<ComplexT> out(static_cast<Eigen::Index>(coords.size()));
    for (std::size_t ii = 0; ii < coords.size(); ++ii)
        out(static_cast<Eigen::Index>(ii)) = DecodeComplexScalar<ComplexT>(coords[ii].as_array());
    return out;
}

inline std::string CanonicalName(SuccessCode code)
{
    switch (code)
    {
        case SuccessCode::NeverStarted: return "NeverStarted";
        case SuccessCode::Success: return "Success";
        case SuccessCode::HigherPrecisionNecessary: return "HigherPrecisionNecessary";
        case SuccessCode::ReduceStepSize: return "ReduceStepSize";
        case SuccessCode::GoingToInfinity: return "GoingToInfinity";
        case SuccessCode::FailedToConverge: return "FailedToConverge";
        case SuccessCode::MatrixSolveFailure: return "MatrixSolveFailure";
        case SuccessCode::MatrixSolveFailureFirstPartOfPrediction: return "MatrixSolveFailureFirstPartOfPrediction";
        case SuccessCode::MaxNumStepsTaken: return "MaxNumStepsTaken";
        case SuccessCode::MaxPrecisionReached: return "MaxPrecisionReached";
        case SuccessCode::MinStepSizeReached: return "MinStepSizeReached";
        case SuccessCode::Failure: return "Failure";
        case SuccessCode::SingularStartPoint: return "SingularStartPoint";
        case SuccessCode::ExternallyTerminated: return "ExternallyTerminated";
        case SuccessCode::MinTrackTimeReached: return "MinTrackTimeReached";
        case SuccessCode::SecurityMaxNormReached: return "SecurityMaxNormReached";
        case SuccessCode::CycleNumTooHigh: return "CycleNumTooHigh";
        case SuccessCode::FailedToSelectPrecisionAndStepsize: return "FailedToSelectPrecisionAndStepsize";
        case SuccessCode::WallClockLimitReached: return "WallClockLimitReached";
    }
    return "UnknownSuccessCode";
}

inline std::string CoarsePathStatus(SuccessCode code)
{
    if (code == SuccessCode::Success)
        return "success";
    if (code == SuccessCode::GoingToInfinity || code == SuccessCode::SecurityMaxNormReached)
        return "diverged";
    return "failed";
}

template <typename ComplexT>
boost::json::object EncodeFullPathResult(parallel::FullPathResult<ComplexT> const& r)
{
    boost::json::object out;
    out["pre_endgame_success_code"] = static_cast<std::int64_t>(r.pre_endgame_success_code);
    out["pre_endgame_success_code_name"] = CanonicalName(r.pre_endgame_success_code);
    out["boundary_point"] = EncodePoint(r.boundary_point);
    out["boundary_stepsize"] = EncodeRealScalar(r.boundary_stepsize);
    out["boundary_precision"] = static_cast<std::int64_t>(r.boundary_precision);
    out["endgame_success_code"] = static_cast<std::int64_t>(r.endgame_success_code);
    out["endgame_success_code_name"] = CanonicalName(r.endgame_success_code);
    out["endpoint"] = EncodePoint(r.solution);
    out["function_residual"] = ExactDoubleText(r.function_residual);
    out["condition_number"] = ExactDoubleText(r.condition_number);
    {
        boost::json::array sv;
        for (Eigen::Index i = 0; i < r.singular_values.size(); ++i)
            sv.push_back(boost::json::value(ExactDoubleText(r.singular_values(i))));
        out["singular_values"] = std::move(sv);
    }
    out["newton_residual"] = ExactDoubleText(r.newton_residual);
    out["final_time_used"] = EncodeComplexScalar(r.final_time_used);
    out["accuracy_estimate_internal_coords"] = ExactDoubleText(r.accuracy_estimate_internal_coords);
    out["accuracy_estimate_user_coords"] = ExactDoubleText(r.accuracy_estimate_user_coords);
    out["cycle_num"] = static_cast<std::int64_t>(r.cycle_num);
    out["precision_digits"] = static_cast<std::int64_t>(r.precision_digits);
    out["accuracy_digits"] = static_cast<std::int64_t>(r.accuracy_digits);
    out["precision_changed"] = r.precision_changed;
    out["time_of_first_prec_increase"] = EncodeComplexScalar(r.time_of_first_prec_increase);
    out["max_precision_used"] = static_cast<std::int64_t>(r.max_precision_used);
    out["path_time_seconds"] = ExactDoubleText(r.path_time_seconds);
    // where the path got to: the step tally always, the point only when the path did not
    // succeed (a successful path's point is `endpoint`)
    out["num_successful_steps"] = static_cast<std::int64_t>(r.num_successful_steps);
    out["num_failed_steps"] = static_cast<std::int64_t>(r.num_failed_steps);
    if (r.latest_path_point.size() > 0)
        out["latest_path_point"] = EncodePoint(r.latest_path_point);
    // the wall-clock budget the path ran under, when there was one: what a later run compares
    // its own patience against before reusing an abandonment
    if (r.wall_clock_limit_seconds > 0)
        out["wall_clock_limit_seconds"] = ExactDoubleText(r.wall_clock_limit_seconds);
    // the crossing verdict, written only when the library could NOT vouch for the endpoint.
    // Absent means clean -- which is also what every record written before this field existed
    // meant, since a store from then has no way to say otherwise.
    if (r.crossing_unresolved)
        out["crossing_unresolved"] = true;
    return out;
}

template <typename ComplexT>
parallel::FullPathResult<ComplexT> DecodeFullPathResult(boost::json::object const& rec,
                                                        std::size_t path_index)
{
    using RealT = typename parallel::FullPathResult<ComplexT>::RealT;
    parallel::FullPathResult<ComplexT> r;
    r.path_index = path_index;
    r.pre_endgame_success_code =
        static_cast<SuccessCode>(rec.at("pre_endgame_success_code").as_int64());
    r.boundary_point = DecodePoint<ComplexT>(rec.at("boundary_point").as_array());
    if constexpr (std::is_same<RealT, double>::value)
        r.boundary_stepsize = DoubleFromText(std::string(rec.at("boundary_stepsize").as_string()));
    else
        r.boundary_stepsize = RealT(std::string(rec.at("boundary_stepsize").as_string()));
    r.boundary_precision = static_cast<unsigned>(rec.at("boundary_precision").as_int64());
    r.endgame_success_code = static_cast<SuccessCode>(rec.at("endgame_success_code").as_int64());
    r.solution = DecodePoint<ComplexT>(rec.at("endpoint").as_array());
    r.function_residual = DoubleFromText(std::string(rec.at("function_residual").as_string()));
    r.condition_number = DoubleFromText(std::string(rec.at("condition_number").as_string()));
    // records written before the spectrum was archived simply lack it
    if (auto const* sv = rec.if_contains("singular_values"))
    {
        auto const& arr = sv->as_array();
        r.singular_values.resize(static_cast<Eigen::Index>(arr.size()));
        for (std::size_t i = 0; i < arr.size(); ++i)
            r.singular_values(static_cast<Eigen::Index>(i)) = DoubleFromText(std::string(arr[i].as_string()));
    }
    r.newton_residual = DoubleFromText(std::string(rec.at("newton_residual").as_string()));
    r.final_time_used = DecodeComplexScalar<ComplexT>(rec.at("final_time_used").as_array());
    // the internal estimate was recorded under the bare key `accuracy_estimate` until the
    // field was named for its coordinates (ADR-0069); a record written before that reads the
    // same, under the key it was written with
    r.accuracy_estimate_internal_coords = DoubleFromText(std::string(
        (rec.if_contains("accuracy_estimate_internal_coords")
            ? rec.at("accuracy_estimate_internal_coords")
            : rec.at("accuracy_estimate")).as_string()));
    r.accuracy_estimate_user_coords =
        DoubleFromText(std::string(rec.at("accuracy_estimate_user_coords").as_string()));
    r.cycle_num = static_cast<unsigned>(rec.at("cycle_num").as_int64());
    r.precision_digits = static_cast<unsigned>(rec.at("precision_digits").as_int64());
    r.accuracy_digits = static_cast<unsigned>(rec.at("accuracy_digits").as_int64());
    r.precision_changed = rec.at("precision_changed").as_bool();
    r.time_of_first_prec_increase =
        DecodeComplexScalar<ComplexT>(rec.at("time_of_first_prec_increase").as_array());
    r.max_precision_used = static_cast<unsigned>(rec.at("max_precision_used").as_int64());
    r.path_time_seconds = DoubleFromText(std::string(rec.at("path_time_seconds").as_string()));
    // records written before paths were stamped with where they got to simply lack these
    if (auto const* v = rec.if_contains("num_successful_steps"))
        r.num_successful_steps = static_cast<unsigned>(v->as_int64());
    if (auto const* v = rec.if_contains("num_failed_steps"))
        r.num_failed_steps = static_cast<unsigned>(v->as_int64());
    if (auto const* v = rec.if_contains("latest_path_point"))
        r.latest_path_point = DecodePoint<ComplexT>(v->as_array());
    if (auto const* v = rec.if_contains("wall_clock_limit_seconds"))
        r.wall_clock_limit_seconds = DoubleFromText(std::string(v->as_string()));
    if (auto const* v = rec.if_contains("crossing_unresolved"))
        r.crossing_unresolved = v->as_bool();
    return r;
}
 
} // namespace records
} // namespace bertini