include/Pythia8/MathTools.h

0001 // MathTools.h is a part of the PYTHIA event generator.
0002 // Copyright (C) 2024 Torbjorn Sjostrand.
0003 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
0004 // Please respect the MCnet Guidelines, see GUIDELINES for details.
0005
0006 // Header file for some mathematics tools, like special functions.
0007 #ifndef Pythia8_MathTools_H
0008 #define Pythia8_MathTools_H
0009
0010 // Header file for the MathTools methods.
0011 #include "Pythia8/Basics.h"
0012 #include "Pythia8/PythiaStdlib.h"
0013
0014 namespace Pythia8 {
0015
0016 //==========================================================================
0017
0018 // The Gamma function for real argument.
0019 double gammaReal(double x);
0020
0021 // Modified Bessel functions of the first and second kinds.
0022 double besselI0(double x);
0023 double besselI1(double x);
0024 double besselK0(double x);
0025 double besselK1(double x);
0026
0027 // Integrate f(x) dx over the specified range.
0028 bool integrateGauss(double& resultOut, function<double(double)> f,
0029   double xLo, double xHi, double tol=1e-6);
0030
0031 // Solve f(x) = target for x in the specified range.
0032 bool brent(double& solutionOut, function<double(double)> f,
0033   double target, double xLo, double xHi, double tol=1e-6, int maxIter = 10000);
0034
0035 // Gram determinant, invariants used in the argument = 2*pi*pj.
0036 double gramDet(double s01tilde, double s12tilde, double s02tilde,
0037   double m0, double m1, double m2);
0038 double gramDet(Vec4 p0, Vec4 p1, Vec4 p2);
0039
0040 // Dilogarithm.
0041 double Li2 (const double, const double kmax = 100.0, const double xerr = 1e-9);
0042
0043 // Standard factorial.
0044 double factorial(const int);
0045
0046 // Binomial coefficient.
0047 int binomial (const int,int);
0048
0049 // Lambert W function.
0050 double lambertW (const double x);
0051
0052 // Kallen function.
0053 double kallenFunction(const double x, const double y, const double z);
0054
0055 // Generate linearly or logarithmically spaced points.
0056 vector<double> linSpace(int nPts, double xMin, double xMax);
0057 vector<double> logSpace(int nPts, double xMin, double xMax);
0058
0059 //==========================================================================
0060
0061 // LinearInterpolator class.
0062 // Used to interpolate between values in linearly spaced data.
0063 // The interpolation uses the arithmetic mean.
0064
0065 class LinearInterpolator {
0066
0067 public:
0068
0069   LinearInterpolator() = default;
0070
0071   // Constructor.
0072   LinearInterpolator(double leftIn, double rightIn, vector<double> ysIn)
0073     : leftSave(leftIn), rightSave(rightIn), ysSave(ysIn) { }
0074
0075   // Function to get y-values of interpolation data.
0076   const vector<double>& data() const { return ysSave; }
0077
0078   // x-values are linearly spaced on the interpolation region.
0079   double left()  const { return leftSave; }
0080   double right() const { return rightSave; }
0081   double dx()    const { return (rightSave - leftSave) / (ysSave.size() - 1); }
0082   double xi(int i) const { return leftSave + i * dx(); }
0083
0084   // Get interpolated value at the specified point.
0085   double at(double x) const;
0086   double operator()(double x) const { return at(x); }
0087
0088   // Plot the data points of this LinearInterpolator in a histogram.
0089   Hist plot(string title) const;
0090   Hist plot(string title, double xMin, double xMax) const;
0091
0092   // Sample a random number distributed as given by this LinearInterpolator.
0093   double sample(Rndm& rndm) const;
0094
0095 private:
0096
0097   // Data members
0098   double leftSave, rightSave;
0099   vector<double> ysSave;
0100
0101 };
0102
0103 //==========================================================================
0104
0105 // LogInterpolator class.
0106 // Used to interpolate between values in logarithmically spaced data.
0107 // The interpolation uses the geometric mean.
0108
0109 class LogInterpolator {
0110
0111 public:
0112
0113   // Default constructor.
0114   LogInterpolator() = default;
0115
0116   // Constructor.
0117   LogInterpolator(double leftIn, double rightIn, vector<double> ysIn)
0118     : leftSave(leftIn), rightSave(rightIn), ysSave(ysIn) {
0119       if (ysIn.size() <= 1)
0120         rxSave = numeric_limits<double>::quiet_NaN();
0121       else
0122         rxSave = pow(rightSave / leftSave, 1. / (ysSave.size() - 1));
0123     }
0124
0125   // Function to get y-values of interpolation data.
0126   const vector<double>& data() const { return ysSave; }
0127
0128   // x-values are logarithmically spaced on the interpolation region.
0129   double left()  const { return leftSave; }
0130   double right() const { return rightSave; }
0131   double rx()    const { return rxSave; }
0132
0133   // Get interpolated value at the specified point.
0134   double at(double x) const;
0135   double operator()(double x) const { return at(x); }
0136
0137   // Plot the data points of this LogInterpolator in a histogram.
0138   Hist plot(string title, int nBins, double xMin, double xMax) const;
0139
0140 private:
0141
0142   // Data members.
0143   double leftSave, rightSave, rxSave;
0144   vector<double> ysSave;
0145
0146 };
0147
0148 //==========================================================================
0149
0150 // Class for the "Hungarian" pairing algorithm. Adapted for Vincia
0151 // from an implementation by M. Buehren and C. Ma, see notices below.
0152
0153 // This is a C++ wrapper with slight modification of a hungarian algorithm
0154 // implementation by Markus Buehren. The original implementation is a few
0155 // mex-functions for use in MATLAB, found here:
0156 // http://www.mathworks.com/matlabcentral/fileexchange/
0157 //    6543-functions-for-the-rectangular-assignment-problem
0158 //
0159 // Both this code and the orignal code are published under the BSD license.
0160 // by Cong Ma, 2016.
0161 //
0162 // Copyright (c) 2014, Markus Buehren
0163 // All rights reserved.
0164 //
0165 // Redistribution and use in source and binary forms, with or without
0166 // modification, are permitted provided that the following conditions are
0167 // met:
0168 //
0169 // * Redistributions of source code must retain the above copyright
0170 // notice, this list of conditions and the following disclaimer.
0171 // * Redistributions in binary form must reproduce the above copyright
0172 // notice, this list of conditions and the following disclaimer in
0173 // the documentation and/or other materials provided with the distribution
0174 //
0175 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
0176 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
0177 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
0178 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
0179 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
0180 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
0181 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
0182 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
0183 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
0184 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
0185 // POSSIBILITY OF SUCH DAMAGE.
0186
0187 class HungarianAlgorithm {
0188
0189 public:
0190
0191   // Function wrapper for solving assignment problem.
0192   double solve(vector <vector<double> >& distMatrix, vector<int>& assignment);
0193
0194  private:
0195
0196   // Solve optimal solution for assignment problem using Munkres algorithm,
0197   // also known as the Hungarian algorithm.
0198   void optimal(vector<int>& assignment, double& cost,
0199     vector<double>& distMatrix, int nOfRows, int nOfColumns);
0200   // Build the assignment vector.
0201   void vect(vector<int>& assignment, vector<bool>& starMatrix, int nOfRows,
0202     int nOfColumns);
0203   // Calculate the assignment cost.
0204   void calcCost(vector<int>& assignment, double& cost,
0205     vector<double>& distMatrix, int nOfRows);
0206   // Factorized step 2a of the algorithm.
0207   void step2a(vector<int>& assignment, vector<double>& distMatrix,
0208     vector<bool>& starMatrix, vector<bool>& newStarMatrix,
0209     vector<bool>& primeMatrix, vector<bool>& coveredColumns,
0210     vector<bool>& coveredRows, int nOfRows, int nOfColumns, int minDim);
0211   // Factorized step 2b of the algorithm.
0212   void step2b(vector<int>& assignment, vector<double>& distMatrix,
0213     vector<bool>& starMatrix, vector<bool>& newStarMatrix,
0214     vector<bool>& primeMatrix, vector<bool>& coveredColumns,
0215     vector<bool>& coveredRows, int nOfRows, int nOfColumns, int minDim);
0216   // Factorized step 3 of the algorithm.
0217   void step3(vector<int>& assignment, vector<double>& distMatrix,
0218     vector<bool>& starMatrix, vector<bool>& newStarMatrix,
0219     vector<bool>& primeMatrix, vector<bool>& coveredColumns,
0220     vector<bool>& coveredRows, int nOfRows, int nOfColumns, int minDim);
0221   // Factorized step 4 of the algorithm.
0222   void step4(vector<int>& assignment, vector<double>& distMatrix,
0223     vector<bool>& starMatrix, vector<bool>& newStarMatrix,
0224     vector<bool>& primeMatrix, vector<bool>& coveredColumns,
0225     vector<bool>& coveredRows, int nOfRows, int nOfColumns, int minDim,
0226     int row, int col);
0227   // Factorized step 5 of the algorithm.
0228   void step5(vector<int>& assignment, vector<double>& distMatrix,
0229     vector<bool>& starMatrix, vector<bool>& newStarMatrix,
0230     vector<bool>& primeMatrix, vector<bool>& coveredColumns,
0231     vector<bool>& coveredRows, int nOfRows, int nOfColumns, int minDim);
0232
0233 };
0234
0235 //==========================================================================
0236
0237 } // end namespace Pythia8
0238
0239 #endif // Pythia8_MathTools_H