minimum_tracer.h

// Copyright (C) 2024 Lisa Biermann, Margarete Mühlleitner, Rui Santos, João
// Viana
//
// SPDX-FileCopyrightText: 2024 Lisa Biermann, Margarete Mühlleitner, Rui
// Santos, João Viana
//
// SPDX-License-Identifier: GPL-3.0-or-later
 
#pragma once
 
#include "Eigen/Eigenvalues"                   // Eigenvalues utility
#include <BSMPT/minimizer/Minimizer.h>         // for Minimizer
#include <BSMPT/models/ClassPotentialOrigin.h> // for Class_Potential_Origin
#include <BSMPT/utility/Logger.h>              // for Logger Class
#include <BSMPT/utility/asciiplotter/asciiplotter.h>
#include <BSMPT/utility/utility.h>
#include <Eigen/Dense> // Eigenvalues matrix
#include <chrono>
#include <cmath>    // std::pow
#include <memory>   // for shared_ptr
#include <optional> // std::optional
#include <stdlib.h> // std::strtol
 
namespace BSMPT
{
 
const double EmptyValue = NAN;
enum class StatusNLOStability
{
  NotSet,
  Off,
  Success,
  NoNLOStability
};
const std::unordered_map<StatusNLOStability, std::string>
    StatusNLOStabilityToString{
        {StatusNLOStability::NotSet, "not_set"},
        {StatusNLOStability::Off, "off"},
        {StatusNLOStability::Success, "success"},
        {StatusNLOStability::NoNLOStability, "no_nlo_stability"}};
enum class StatusEWSR
{
  NotSet,
  Off,
  Failure,
  NotBFB,
  FlatRegion,
  EWSymNonRes,
  EWSymRes
};
const std::unordered_map<StatusEWSR, std::string> StatusEWSRToString{
    {StatusEWSR::NotSet, "not_set"},
    {StatusEWSR::Off, "off"},
    {StatusEWSR::Failure, "failure"},
    {StatusEWSR::NotBFB, "non_bfb"},
    {StatusEWSR::FlatRegion, "flat_region"},
    {StatusEWSR::EWSymNonRes, "ew_syum_non_res"},
    {StatusEWSR::EWSymRes, "ew_sym_res"}};
enum class StatusTracing
{
  NotSet,
  Success,
  NoCoverage,
  NoMinsAtBoundaries,
  NoGlobMinCoverage,
  Failure
};
const std::unordered_map<StatusTracing, std::string> StatusTracingToString{
    {StatusTracing::NotSet, "not_set"},
    {StatusTracing::Success, "success"},
    {StatusTracing::NoCoverage, "no_coverage"},
    {StatusTracing::NoMinsAtBoundaries, "no_mins_at_boundaries"},
    {StatusTracing::NoGlobMinCoverage, "no_glob_min_coverage"},
    {StatusTracing::Failure, "failure"}};
enum class StatusCoexPair
{
  NotSet,
  Success,
  NoCoexPairs
};
const std::unordered_map<StatusCoexPair, std::string> StatusCoexPairToString{
    {StatusCoexPair::NotSet, "not_set"},
    {StatusCoexPair::Success, "success"},
    {StatusCoexPair::NoCoexPairs, "no_coex_pair"}};
enum class StatusCrit
{
  NotSet,
  Success,
  FalseLower,
  TrueLower,
  Failure
};
const std::unordered_map<StatusCrit, std::string> StatusCritToString{
    {StatusCrit::NotSet, "not_set"},
    {StatusCrit::Success, "success"},
    {StatusCrit::FalseLower, "false_lower"},
    {StatusCrit::TrueLower, "true_lower"},
    {StatusCrit::Failure, "failure"}};
enum class StatusTemperature
{
  NotSet,
  Success,
  NotMet,
  NaN
};
const std::unordered_map<StatusTemperature, std::string>
    StatusTemperatureToString{{StatusTemperature::NotSet, "not_set"},
                              {StatusTemperature::Success, "success"},
                              {StatusTemperature::NotMet, "not_met"},
                              {StatusTemperature::NaN, "nan"}};
enum class StatusGW
{
  NotSet,
  Success,
  Failure
};
const std::unordered_map<StatusGW, std::string> StatusGWToString{
    {StatusGW::NotSet, "not_set"},
    {StatusGW::Success, "success"},
    {StatusGW::Failure, "failure"}};
 
std::ostream &operator<<(std::ostream &os, const StatusNLOStability &status);
std::ostream &operator<<(std::ostream &os, const StatusEWSR &status);
std::ostream &operator<<(std::ostream &os, const StatusTracing &status);
std::ostream &operator<<(std::ostream &os, const StatusCoexPair &status);
std::ostream &operator<<(std::ostream &os, const StatusCrit &status);
std::ostream &operator<<(std::ostream &os, const StatusGW &status);
std::ostream &operator<<(std::ostream &os, const StatusTemperature &status);
 
struct Minimum
{
  std::vector<double> point;
  double temp;
  double potential;
  bool is_glob_min = false;
  int EdgeOfPhase  = 0;
 
  bool operator<(const Minimum &a) const { return temp < a.temp; }
};
 
class MinimumTracer
{
private:
  int WhichMinimizer;
  bool UseMultithreading;
 
protected:
  std::shared_ptr<Class_Potential_Origin> modelPointer;
 
public:
  double GradientThreshold = 1e-3;
 
  double HessianDiagonalShift = 1e-3;
 
  std::vector<double> HighTemperatureVEV;
 
  std::vector<Minimum> SavedMinimaFromVEVSplitting;
 
  MinimumTracer();
 
  MinimumTracer(const std::shared_ptr<Class_Potential_Origin> &pointer_in,
                const int &WhichMinimizer_in,
                const bool &UseMultithreading_in);
 
  void FindFlatDirections();
 
  void ConvertToNonFlatDirections(std::vector<double> &point);
 
  void FindDiscreteSymmetries();
 
  std::vector<double> ConvertToVEVDim(const std::vector<double> &point);
 
  std::vector<double> GetGlobalMinimum(const double &Temp,
                                       std::vector<double> &check,
                                       const std::vector<double> &start);
 
  std::vector<double> GetGlobalMinimum(const double &Temp,
                                       const std::vector<double> &start);
 
  std::vector<double> GetGlobalMinimum(const double &Temp);
 
  void IsGlobMin(Minimum &min);
 
  StatusNLOStability GetStatusNLOVEV(const bool &out);
 
  StatusEWSR GetStatusEWSR(const int &out);
 
  int IsThereEWSymmetryRestoration();
 
  double SmallestEigenvalue(
      const std::vector<double> &point,
      const std::function<std::vector<std::vector<double>>(std::vector<double>)>
          &Hessian);
 
  std::vector<double> FindZeroSmallestEigenvalue(std::vector<double> point_1,
                                                 double T_1,
                                                 std::vector<double> point_2,
                                                 double T_2);
 
  std::vector<Minimum> TrackPhase(double &globMinEndT,
                                  const std::vector<double> &point_In,
                                  const double &currentT_In,
                                  const double &finalT,
                                  const double &dT_In     = 1,
                                  const bool &output      = true,
                                  const bool &unprotected = false);
 
  std::vector<Minimum> TrackPhase(const std::vector<double> &point_In,
                                  const double &currentT_In,
                                  const double &finalT,
                                  const double &dT_In     = 1,
                                  const bool &output      = true,
                                  const bool &unprotected = false);
 
  bool flat_dirs_found = false;
 
  std::vector<int> NonFlatDirections;
 
  std::vector<std::size_t> flat_1D_dirs;
 
  std::vector<std::vector<std::size_t>> flat_2D_dirs;
 
  std::vector<std::vector<std::size_t>> flat_3D_dirs;
 
  std::vector<Eigen::MatrixXd> GroupElements;
 
  void ReduceVEV(std::vector<double> &vev);
 
  void ReduceVEV(Minimum &min);
 
  const std::vector<std::vector<double>>
  WarpPath(const std::vector<std::vector<double>> &path,
           const std::vector<double> &T1,
           const std::vector<double> &F1,
           const std::vector<double> &T2,
           const std::vector<double> &F2);
 
  std::vector<double> LocateMinimum(
      const std::vector<double> &guess_In,
      std::function<std::vector<double>(std::vector<double>)> &df,
      std::function<std::vector<std::vector<double>>(std::vector<double>)>
          &Hessian,
      const double &error            = 1e-4,
      const double &const_multiplier = 1e-2,
      const int &maxiter             = 100);
 
  std::vector<std::string> GetLegend(const int &num_coex_phases,
                                     const bool &do_gw_calc);
};
 
std::vector<std::vector<double>>
Create1DimGrid(const std::vector<double> &point,
               const int k,
               const double low_value,
               const double high_value,
               const int nsteps = 100);
 
std::vector<std::vector<double>>
Create1DimGrid(const std::vector<double> &min_start,
               const std::vector<double> &min_end,
               const int npoints = 100);
 
bool almost_the_same(const double &a,
                     const double &b,
                     const double &rel_precision = 0.01,
                     const double &num_zero      = 1e-10);
 
bool almost_the_same(const std::vector<double> &a,
                     const std::vector<double> &b,
                     const bool &allow_for_sign_flip = false,
                     const double &rel_precision     = 0.01,
                     const double &num_zero          = 1e-10);
 
struct Phase
{
  int id = 0;
 
  double T_low = 0;
 
  double T_high = 0;
 
  std::vector<Minimum> MinimumPhaseVector;
 
  std::shared_ptr<MinimumTracer> MinTracer;
 
  void Add(Minimum min);
 
  Minimum Get(double T);
 
  Phase();
 
  Phase(const std::vector<double> &phase_start,
        const double &initialT,
        const double &finalT,
        double &globMinEndT,
        std::shared_ptr<MinimumTracer> &MinTracerIn);
 
  Phase(const double &initialT,
        const std::vector<double> &phase_start,
        const double &LowT,
        const double &HighT,
        double &globMinEndT,
        std::shared_ptr<MinimumTracer> &MinTracerIn);
 
  Phase(const std::vector<double> &phase_start,
        const double &initialT,
        const double &finalT,
        std::shared_ptr<MinimumTracer> &MinTracerIn);
 
  Phase(const double &initialT,
        const std::vector<double> &phase_start,
        const double &LowT,
        const double &HighT,
        std::shared_ptr<MinimumTracer> &MinTracerIn);
 
  Phase(const double &initialT,
        const double &LowT,
        const double &HighT,
        std::shared_ptr<MinimumTracer> &MinTracerIn);
};
 
struct CoexPhases
{
  int coex_pair_id;
  double T_high;
  double T_low;
  Phase false_phase;
  Phase true_phase;
  double crit_temp       = -1;
  StatusCrit crit_status = StatusCrit::NotSet;
 
  CoexPhases();
 
  CoexPhases(const int &pair_id,
             const Phase &false_phase,
             const Phase &true_phase,
             const double &Tlow_in,
             const double &Thigh_in);
 
  void CalculateTc();
};
 
struct Vacuum
{
  double T_low = 0;
 
  double T_high = 300;
 
  double T_low_highTempPhase = -1;
 
  double T_high_lowTempPhase = -1;
 
  int num_points = 0;
 
  StatusTracing status_vacuum = StatusTracing::NotSet;
 
  StatusCoexPair status_coex_pairs = StatusCoexPair::NotSet;
 
  std::shared_ptr<MinimumTracer> MinTracer;
 
  std::shared_ptr<Class_Potential_Origin> modelPointer;
 
  std::vector<Phase> PhasesList;
 
  std::vector<CoexPhases> CoexPhasesList;
 
  Vacuum(const double &T_lowIn,
         const double &T_highIn,
         std::shared_ptr<MinimumTracer> &MinTracerIn,
         std::shared_ptr<Class_Potential_Origin> &modelPointerIn,
         const int &UseMultiStepPTModeIn,
         const int &num_pointsIn     = 10,
         const bool &do_only_tracing = false);
 
  void MultiStepPTTracer(const double &Temp, const double &deltaT = 0);
 
  void print(const Phase &phase);
 
  void setCoexPhases();
 
  void setCoexRegion(const int &UseMultiStepPTMode);
 
  void addPhase(Phase &phase);
 
  int MinimumFoundAlready(const Minimum &minimum);
 
  void MultiStepPTMode0(const std::vector<double> &LowTempPoint,
                        const std::vector<double> &HighTempPoint);
 
  void MultiStepPTMode1(const std::vector<double> &LowTempPoint,
                        const std::vector<double> &HighTempPoint);
 
  bool DoPhasesOverlap(Phase &new_phase, Phase &old_phase);
 
  bool DoGlobMinOverlap(const Phase &new_phase, const Phase &old_phase);
 
  void MultiStepPTMode2(const std::vector<double> &LowTempPoint,
                        const std::vector<double> &HighTempPoint);
 
  void PrintPhasesDiagram(int size = 100);
};
 
} // namespace BSMPT