remove static derivatives

.gitignore

@@ -16,3 +16,5 @@ analyze*.py
 *.synctex.gz
 
 CMakeUserPresets.json
+__pycache__/Build.cpython-310.pyc
+__pycache__/Setup.cpython-310.pyc

include/BSMPT/bounce_solution/action_calculation.h

@@ -334,59 +334,6 @@ class BounceActionInt
    */
   void PrintVector(std::vector<double> vec);
 
-  /**
-   * @brief Numerical method to calculate the potential's (or other functions's)
-   * gradient using finite differences method.
-   *
-   * This method is used while BSMPT is not able to
-   * calculate the potential derivative analytically. We used the 4th order
-   * method
-   *
-   * \f$\frac{\partial V}{\partial \phi_i} = \frac{1}{12
-   * \epsilon}\left(-V(\dots ,\vec{\phi}_i + 2  \epsilon ) + 8 V(\dots
-   * ,\vec{\phi}_i + \epsilon )- 8 V(\dots ,\vec{\phi}_i - \epsilon ) +
-   * V(\dots ,\vec{\phi}_i - 2  \epsilon )\right)\f$
-   *
-   * where \f$ \epsilon \f$ is a small step.
-   *
-   * @param phi Where we want to calculate the gradient
-   * @param V Potential (or other function)
-   * @param eps Size of finite differences step
-   * @param dim Dimensions of the VEV space (or dimensions of V argument)
-   * @return std::vector<double> The \f$ dim \times 1 \f$ gradient of V taken at
-   * phi
-   */
-  static std::vector<double>
-  NablaNumerical(const std::vector<double> &phi,
-                 const std::function<double(std::vector<double>)> &V,
-                 const double &eps,
-                 const int &dim);
-  /**
-   * @brief Numerical method to calculate the potential's (or other functions's)
-   * hessian matrix using finite differences method.
-   *
-   * \f$\frac{\partial^2 V}{\partial \phi_i \phi_j} = \frac{1}{4
-   * \epsilon^2}\left(V(\dots, \vec{\phi}_i + \epsilon , \vec{\phi}_j +
-   * \epsilon) - V(\dots, \vec{\phi}_i - \epsilon , \vec{\phi}_j +
-   * \epsilon) - V(\dots, \vec{\phi}_i + \epsilon , \vec{\phi}_j -
-   * \epsilon) + V(\dots, \vec{\phi}_i - \epsilon , \vec{\phi}_j -
-   * \epsilon) \right)\f$
-   *
-   * where \f$ \epsilon \f$ is a small step.
-   *
-   * @param phi Where we want to calculate the Hessian matrix
-   * @param V Potential (or other function)
-   * @param eps Size of finite differences step
-   * @param dim Dimensions of the VEV space (or dimensions of V argument)
-   * @return std::vector<std::vector<double>> The \f$ dim \times \dim \f$
-   *  hessian matrix of V taken at phi
-   */
-  static std::vector<std::vector<double>>
-  HessianNumerical(const std::vector<double> &phi,
-                   const std::function<double(std::vector<double>)> &V,
-                   double eps,
-                   const int &dim);
-
   /**
    * @brief Calculates \f$ \frac{dV}{dl} \f$ using the spline and potential
    * derivatives

include/BSMPT/minimum_tracer/minimum_tracer.h

@@ -279,6 +279,13 @@ class MinimumTracer
    */
   double GradientThreshold = 1e-3;
 
+  /**
+   * @brief Add a constant to the diagonals of the hessian matrix in the
+   * LocateMinimum function. Helps with convergence.
+   *
+   */
+  double HessianDiagonalShift = 1e-3;
+
   /**
    * @brief Minimum found in IsThereEWSymmetryRestoration()
    *
@@ -544,60 +551,6 @@ class MinimumTracer
                                      const bool &do_gw_calc);
 };
 
-/**
- * @brief Numerical method to calculate the
- * gradient of a function f using finite differences method.
- *
- * This method is used while BSMPT is not able to
- * calculate the potential derivative analytically. We used the 4th order
- * method
- *
- * \f$\frac{\partial f}{\partial \phi_i} = \frac{1}{12
- * \epsilon}\left(-f(\dots ,\vec{\phi}_i + 2  \epsilon ) + 8 f(\dots
- * ,\vec{\phi}_i + \epsilon )- 8 f(\dots ,\vec{\phi}_i - \epsilon ) +
- * f(\dots ,\vec{\phi}_i - 2  \epsilon )\right)\f$
- *
- * where \f$ \epsilon \f$ is a small step.
- *
- * @param phi Where we want to calculate the gradient
- * @param f function
- * @param eps Size of finite differences step
- * @param dim Dimensions of the VEV space (or dimensions of V argument)
- * @return std::vector<double> The \f$ dim \times 1 \f$ gradient of V taken at
- * phi
- */
-std::vector<double>
-NablaNumerical(const std::vector<double> &phi,
-               const std::function<double(std::vector<double>)> &f,
-               const double &eps,
-               const int &dim);
-
-/**
- * @brief Numerical method to calculate the
- * hessian matrix of a function f using finite differences method.
- *
- * \f$\frac{\partial^2 f}{\partial \phi_i \phi_j} = \frac{1}{4
- * \epsilon^2}\left(V(\dots, \vec{\phi}_i + \epsilon , \vec{\phi}_j +
- * \epsilon) - f(\dots, \vec{\phi}_i - \epsilon , \vec{\phi}_j +
- * \epsilon) - f(\dots, \vec{\phi}_i + \epsilon , \vec{\phi}_j -
- * \epsilon) + f(\dots, \vec{\phi}_i - \epsilon , \vec{\phi}_j -
- * \epsilon) \right)\f$
- *
- * where \f$ \epsilon \f$ is a small step.
- *
- * @param phi Where we want to calculate the Hessian matrix
- * @param f Potential (or other function)
- * @param eps Size of finite differences step
- * @param dim Dimensions of the VEV space (or dimensions of f argument)
- * @return std::vector<std::vector<double>> The \f$ dim \times \dim \f$
- *  hessian matrix of f taken at phi
- */
-std::vector<std::vector<double>>
-HessianNumerical(const std::vector<double> &phi,
-                 const std::function<double(std::vector<double>)> &f,
-                 const double &eps,
-                 const int &dim);
-
 /**
  * @brief Create1DimGrid creates a 1-dim grid of given size in index-direction
  * @param point

include/BSMPT/utility/NumericalDerivatives.h

@@ -0,0 +1,56 @@
+#pragma once
+#include <functional>
+#include <vector>
+
+namespace BSMPT
+{
+/**
+ * @brief Numerical method to calculate the
+ * gradient of a function f using finite differences method.
+ *
+ * This method is used while BSMPT is not able to
+ * calculate the potential derivative analytically. We used the 4th order
+ * method
+ *
+ * \f$\frac{\partial f}{\partial \phi_i} = \frac{1}{12
+ * \epsilon}\left(-f(\dots ,\vec{\phi}_i + 2  \epsilon ) + 8 f(\dots
+ * ,\vec{\phi}_i + \epsilon )- 8 f(\dots ,\vec{\phi}_i - \epsilon ) +
+ * f(\dots ,\vec{\phi}_i - 2  \epsilon )\right)\f$
+ *
+ * where \f$ \epsilon \f$ is a small step.
+ *
+ * @param phi Where we want to calculate the gradient
+ * @param f function
+ * @param eps Size of finite differences step
+ * @return std::vector<double> The \f$ dim \times 1 \f$ gradient of V taken at
+ * phi
+ */
+std::vector<double>
+NablaNumerical(const std::vector<double> &phi,
+               const std::function<double(std::vector<double>)> &f,
+               const double &eps);
+
+/**
+ * @brief Numerical method to calculate the potential's (or other functions's)
+ * hessian matrix using finite differences method.
+ *
+ * \f$\frac{\partial^2 V}{\partial \phi_i \phi_j} = \frac{1}{4
+ * \epsilon^2}\left(V(\dots, \vec{\phi}_i + \epsilon , \vec{\phi}_j +
+ * \epsilon) - V(\dots, \vec{\phi}_i - \epsilon , \vec{\phi}_j +
+ * \epsilon) - V(\dots, \vec{\phi}_i + \epsilon , \vec{\phi}_j -
+ * \epsilon) + V(\dots, \vec{\phi}_i - \epsilon , \vec{\phi}_j -
+ * \epsilon) \right)\f$
+ *
+ * where \f$ \epsilon \f$ is a small step.
+ *
+ * @param phi Where we want to calculate the Hessian matrix
+ * @param V Potential (or other function)
+ * @param eps Size of finite differences step
+ * @return std::vector<std::vector<double>> The \f$ dim \times \dim \f$
+ *  hessian matrix of V taken at phi
+ */
+std::vector<std::vector<double>>
+HessianNumerical(const std::vector<double> &phi,
+                 const std::function<double(std::vector<double>)> &V,
+                 double eps);
+} // namespace BSMPT

src/bounce_solution/action_calculation.cpp

@@ -8,6 +8,7 @@
  */
 
 #include <BSMPT/bounce_solution/action_calculation.h>
+#include <BSMPT/utility/NumericalDerivatives.h>
 
 namespace BSMPT
 {
@@ -29,7 +30,7 @@ BounceActionInt::BounceActionInt(
   this->V  = [&](std::vector<double> vev) { return V_In(vev) - this->Vfalse; };
   this->dV = dV_In;
   this->Hessian = [=](auto const &arg)
-  { return HessianNumerical(arg, V_In, this->eps, this->dim); };
+  { return HessianNumerical(arg, V_In, this->eps); };
   this->TrueVacuum          = TrueVacuum_In;
   this->FalseVacuum         = FalseVacuum_In;
   this->InitPath            = InitPath_In;
@@ -53,9 +54,9 @@ BounceActionInt::BounceActionInt(
   this->V = [&](std::vector<double> vev) { return V_In(vev) - this->Vfalse; };
   // Use numerical derivative
   this->dV = [=](auto const &arg)
-  { return NablaNumerical(arg, this->V, this->eps, this->dim); };
+  { return NablaNumerical(arg, this->V, this->eps); };
   this->Hessian = [=](auto const &arg)
-  { return HessianNumerical(arg, V_In, this->eps, this->dim); };
+  { return HessianNumerical(arg, V_In, this->eps); };
   this->TrueVacuum          = TrueVacuum_In;
   this->FalseVacuum         = FalseVacuum_In;
   this->InitPath            = InitPath_In;
@@ -118,75 +119,6 @@ void BounceActionInt::PrintVector(std::vector<double> vec)
   BSMPT::Logger::Write(BSMPT::LoggingLevel::BounceDetailed, ss.str());
 }
 
-std::vector<double> BounceActionInt::NablaNumerical(
-    const std::vector<double> &phi,
-    const std::function<double(std::vector<double>)> &V,
-    const double &eps,
-    const int &dim)
-{
-  // Numerical gradient of the potential up to second order
-  std::vector<double> result(dim);
-
-  for (int i = 0; i < dim; i++)
-  {
-    std::vector<double> lp2 = phi;
-    lp2[i] += 2 * eps;
-    std::vector<double> lp1 = phi;
-    lp1[i] += eps;
-    std::vector<double> lm1 = phi;
-    lm1[i] -= eps;
-    std::vector<double> lm2 = phi;
-    lm2[i] -= 2 * eps;
-    result[i] = (-V(lp2) + 8 * V(lp1) - 8 * V(lm1) + V(lm2)) / (12 * eps);
-  }
-  return result;
-}
-
-std::vector<std::vector<double>> BounceActionInt::HessianNumerical(
-    const std::vector<double> &phi,
-    const std::function<double(std::vector<double>)> &V,
-    double eps,
-    const int &dim)
-{
-  std::vector<std::vector<double>> result(dim, std::vector<double>(dim));
-  for (int i = 0; i < dim; i++)
-  {
-    // https://en.wikipedia.org/wiki/Finite_difference
-    for (int j = i; j < dim; j++)
-    {
-      double r = 0;
-
-      if (i == j) eps /= 2;
-
-      std::vector<double> xp = phi; // F(x+h, y+h)
-      xp[i] += eps;
-      xp[j] += eps;
-      r += V(xp);
-
-      xp = phi; //-F(x+h, y-h)
-      xp[i] += eps;
-      xp[j] -= eps;
-      r -= V(xp);
-
-      xp = phi; //-F(x-h, y+h)
-      xp[i] -= eps;
-      xp[j] += eps;
-      r -= V(xp);
-
-      xp = phi; // F(x-h, y-h)
-      xp[i] -= eps;
-      xp[j] -= eps;
-      r += V(xp);
-
-      result[i][j] = r / (4 * eps * eps);
-      result[j][i] = r / (4 * eps * eps);
-
-      if (i == j) eps *= 2;
-    }
-  }
-  return result;
-}
-
 std::vector<double>
 BounceActionInt::NormalForce(const double &l,
                              const double &dldrho,