classic_writer.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
//Bertini 2 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.
//
//Bertini 2 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 Bertini 2. If not, see <http://www.gnu.org/licenses/>.
//
// Copyright(C) 2015 - 2026 by Bertini2 developers
//
// 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/system/system.hpp"
#include "bertini2/function_tree/find.hpp"
#include "bertini2/trackers/config.hpp"
 
#include <ostream>
#include <sstream>
#include <iomanip>
#include <string>
#include <vector>
 
namespace bertini{
    namespace classic{

        inline void EmitSystem(std::ostream& out, System const& sys)
        {
            auto emit_groups = [&out](char const* keyword, auto const& groups){
                for (auto const& grp : groups)
                {
                    out << keyword << " ";
                    for (size_t i = 0; i < grp.size(); ++i)
                        out << (i ? ", " : "") << *grp[i];
                    out << ";\n";
                }
            };
            emit_groups("variable_group", sys.VariableGroups());
            emit_groups("hom_variable_group", sys.HomVariableGroups());
 
            // a homotopy is a system too: its path variable must be declared, or the functions that
            // use it do not parse back and the file does not round-trip (b2#366)
            if (sys.HavePathVariable())
                out << "pathvariable " << *sys.GetPathVariable() << ";\n";
 
            auto functions = sys.NaturalFunctionsAsNodes();
 
            // Named subexpressions are not stored separately; they are discovered inside the function
            // trees (nested ones included).  Emit each as `name = expr;` before the functions, the way
            // the classic parser expects a subfunction to be defined ahead of its use.
            {
                std::vector<std::shared_ptr<const node::Node>> roots(functions.begin(), functions.end());
                for (auto const& ne : node::Find<node::NamedExpression>(roots))
                    out << ne->name() << " = " << ne->EntryNode() << ";\n";
            }
 
            out << "function ";
            for (size_t i = 0; i < functions.size(); ++i)
                out << (i ? ", " : "") << "f" << i;
            out << ";\n";
            for (size_t i = 0; i < functions.size(); ++i)
                out << "f" << i << " = " << functions[i] << ";\n";
        }

        inline std::string SystemToClassic(System const& sys)
        {
            std::ostringstream ss;
            EmitSystem(ss, sys);
            return ss.str();
        }
 

        struct ClassicWriteOptions
        {
            int           tracktype              = 0;      
            int           mptype                 = 2;      
            int           odepredictor           = 5;      
            double        tracktolbeforeeg        = 1e-5;  
            double        tracktolduringeg        = 1e-6;  
            double        finaltol                = 1e-11; 
            // step-size cadence -- governs how aggressively a path may take a big step (i.e. jump),
            // so it must be controlled for a fair crossing-rate comparison against Bertini 1.
            double        maxstepsize             = 0.1;   
            double        stepsuccessfactor       = 2.0;   
            double        stepfailfactor          = 0.5;   
            unsigned      stepsforincrease         = 5;     
            unsigned long maxnumbersteps          = 100000;
            unsigned      maxnewtonits             = 2;     
            unsigned      maxcrossedpathresolves   = 2;     
        };

        inline int PredictorToClassic(tracking::Predictor p)
        {
            using P = tracking::Predictor;
            switch (p)
            {
                case P::Constant:          return -1;
                case P::Euler:             return 0;
                case P::Heun:              return 1;
                case P::RK4:               return 2;
                case P::HeunEuler:         return 1;
                case P::RKNorsett34:       return 4;
                case P::RKF45:             return 5;
                case P::RKCashKarp45:      return 6;
                case P::RKDormandPrince56: return 7;
                case P::RKVerner67:        return 8;
            }
            return 5; // RKF45 -- the Bertini 2 default, a safe fallback
        }

        inline void EmitConfig(std::ostream& out, System const& sys, ClassicWriteOptions const& opt)
        {
            auto num = [](double v){ std::ostringstream s; s << std::setprecision(15) << v; return s.str(); };
            out << "tracktype: "              << opt.tracktype              << ";\n";
            out << "mptype: "                 << opt.mptype                 << ";\n";
            out << "odepredictor: "           << opt.odepredictor           << ";\n";
            out << "tracktolbeforeeg: "       << num(opt.tracktolbeforeeg)  << ";\n";
            out << "tracktolduringeg: "       << num(opt.tracktolduringeg)  << ";\n";
            out << "finaltol: "               << num(opt.finaltol)          << ";\n";
            out << "maxstepsize: "            << num(opt.maxstepsize)       << ";\n";
            out << "stepsuccessfactor: "      << num(opt.stepsuccessfactor) << ";\n";
            out << "stepfailfactor: "         << num(opt.stepfailfactor)    << ";\n";
            out << "stepsforincrease: "       << opt.stepsforincrease       << ";\n";
            out << "maxnumbersteps: "         << opt.maxnumbersteps         << ";\n";
            out << "maxnewtonits: "           << opt.maxnewtonits           << ";\n";
            out << "maxcrossedpathresolves: " << opt.maxcrossedpathresolves << ";\n";
            out << "coefficientbound: "       << num(static_cast<double>(sys.CoefficientBound<complex_dbl>())) << ";\n";
            // A system that is not a polynomial one has no degree bound, and asking for one now
            // throws.  Write the file anyway -- "this is what Bertini 1 would eat, if it could
            // eat it" is a useful thing to look at -- but say why the setting is absent instead
            // of inventing a number.  `%` is a comment in this format, and our own parser strips
            // comments, so the file still reads back into b2.
            if (sys.IsPolynomial())
                out << "degreebound: "        << sys.DegreeBound()          << ";\n";
            else
                out << "% degreebound: omitted -- this system is not polynomial, so it has none\n";
        }

        inline void WriteClassicInput(std::ostream& out, System const& sys,
                                      ClassicWriteOptions const& opt = ClassicWriteOptions{})
        {
            out << "CONFIG\n";
            EmitConfig(out, sys, opt);
            out << "END;\n\nINPUT\n";
            EmitSystem(out, sys);
            out << "END;\n";
        }

        inline std::string SystemToClassicFile(System const& sys,
                                               ClassicWriteOptions const& opt = ClassicWriteOptions{})
        {
            std::ostringstream ss;
            WriteClassicInput(ss, sys, opt);
            return ss.str();
        }
 
    } // namespace classic
} // namespace bertini