first cmake draft, source relocation, not working

lib/LatAnalyze/Numerical/MinuitMinimizer.cpp

@@ -0,0 +1,198 @@
+/*
+ * MinuitMinimizer.cpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2020 Antonin Portelli
+ *
+ * LatAnalyze 3 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.
+ *
+ * LatAnalyze 3 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 LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <LatAnalyze/Numerical/MinuitMinimizer.hpp>
+#include <LatAnalyze/includes.hpp>
+
+// forward declaration necessary in the ROOT-based version of Minuit2
+namespace ROOT
+{
+    namespace Fit
+    {
+        class ParameterSettings;
+    };
+};
+
+// macros necessary in the ROOT-based version of Minuit2
+#define ROOT_Math_VecTypes
+#define MATHCORE_STANDALONE
+
+#include <Minuit2/Minuit2Minimizer.h>
+#include <Math/Functor.h>
+
+using namespace std;
+using namespace Latan;
+
+static constexpr double initErr = 0.1;
+
+/******************************************************************************
+ *                    MinuitMinimizer implementation                          *
+ ******************************************************************************/
+// constructors ////////////////////////////////////////////////////////////////
+MinuitMinimizer::MinuitMinimizer(const Algorithm algorithm)
+{
+    setAlgorithm(algorithm);
+}
+
+// access //////////////////////////////////////////////////////////////////////
+MinuitMinimizer::Algorithm MinuitMinimizer::getAlgorithm(void) const
+{
+    return algorithm_;
+}
+
+void MinuitMinimizer::setAlgorithm(const Algorithm algorithm)
+{
+    algorithm_ = algorithm;
+}
+
+bool MinuitMinimizer::supportLimits(void) const
+{
+    return true;
+}
+
+// minimization ////////////////////////////////////////////////////////////////
+const DVec & MinuitMinimizer::operator()(const DoubleFunction &f)
+{
+    using namespace ROOT;
+    using namespace Minuit2;
+    
+    DVec           &x = getState();
+    int            printLevel = 0;
+    EMinimizerType minuitAlg  = kCombined;
+    double         prec = getPrecision();
+    
+    // convert Latan parameters to Minuit parameters
+    switch (getVerbosity())
+    {
+        case Verbosity::Silent:
+            printLevel = 0;
+            break;
+        case Verbosity::Normal:
+            printLevel = 2;
+            break;
+        case Verbosity::Debug:
+            printLevel = 3;
+            break;
+    }
+    // The factor of 0.002 here is to compensate the dirty hack in Minuit
+    // source used to match the C++ and F77 versions
+    // (cf. VariableMetricBuilder.cxx)
+    switch (getAlgorithm())
+    {
+        case Algorithm::migrad:
+            minuitAlg  = kMigrad;
+            prec      /= 0.002;
+            break;
+        case Algorithm::simplex:
+            minuitAlg = kSimplex;
+            break;
+        case Algorithm::combined:
+            minuitAlg  = kCombined;
+            prec      /= 0.002;
+            break;
+    }
+    
+    // resize minimizer state to match function number of arguments
+    if (f.getNArg() != x.size())
+    {
+        resize(f.getNArg());
+    }
+    
+    // create and set minimizer
+    Minuit2Minimizer min(minuitAlg);
+    
+    min.SetStrategy(2);
+    min.SetMaxFunctionCalls(getMaxIteration());
+    min.SetTolerance(prec);
+    min.SetPrintLevel(printLevel);
+    
+    // set function and variables
+    Math::Functor minuitF(f, x.size());
+    string        name;
+    double        val, step;
+    
+    min.SetFunction(minuitF);
+    for (Index i = 0; i < x.size(); ++i)
+    {
+        name = f.varName().getName(i);
+        val  = x(i);
+        step = (fabs(x(i)) != 0.) ? initErr*fabs(x(i)) : 1.;
+        if (hasHighLimit(i) and !hasLowLimit(i))
+        {
+            min.SetUpperLimitedVariable(i, name, val, step, getHighLimit(i));
+        }
+        else if (!hasHighLimit(i) and hasLowLimit(i))
+        {
+            min.SetLowerLimitedVariable(i, name, val, step, getLowLimit(i));
+        }
+        else if (hasHighLimit(i) and hasLowLimit(i))
+        {
+            min.SetLimitedVariable(i, name, val, step, getLowLimit(i),
+                                   getHighLimit(i));
+        }
+        else
+        {
+            min.SetVariable(i, name, val, step);
+        }
+    }
+    
+    // minimize
+    int          status;
+    unsigned int n = 0;
+    
+    do
+    {
+        if (getVerbosity() >= Verbosity::Normal)
+        {
+            cout << "========== Minuit minimization, pass #" << n + 1;
+            cout << " =========" << endl;
+        }
+        min.Minimize();
+        status = min.Status();
+        n++;
+    } while ((status >= 2) and (n < getMaxPass()));
+    if (getVerbosity() >= Verbosity::Normal)
+    {
+        cout << "=================================================" << endl;
+    }
+    switch (status)
+    {
+        case 1:
+            // covariance matrix was made positive, the minimum is still good
+            // it just means that Minuit error analysis is inaccurate
+            break;
+        case 2:
+            LATAN_WARNING("invalid minimum: Hesse analysis is not valid");
+            break;
+        case 3:
+            LATAN_WARNING("invalid minimum: requested precision not reached");
+            break;
+        case 4:
+            LATAN_WARNING("invalid minimum: iteration limit reached");
+            break;
+    }
+    
+    // save and return result
+    for (Index i = 0; i < x.size(); ++i)
+    {
+        x(i) = min.X()[i];
+    }
+    
+    return x;
+}