interface to NLopt minimisers (tested and working)

configure.ac

@@ -45,6 +45,12 @@ AC_ARG_WITH([Minuit2],
     [AM_CFLAGS="$AM_CFLAGS -I$with_Minuit2/include"]
     [AM_CXXFLAGS="$AM_CXXFLAGS -I$with_Minuit2/include"]
 	[AM_LDFLAGS="$AM_LDFLAGS -L$with_Minuit2/lib"])
+AC_ARG_WITH([nlopt],
+    [AS_HELP_STRING([--with-nlopt=prefix],
+        [try this for a non-standard install prefix of the nlopt library])],
+    [AM_CFLAGS="$AM_CFLAGS -I$with_nlopt/include"]
+    [AM_CXXFLAGS="$AM_CXXFLAGS -I$with_nlopt/include"]
+    [AM_LDFLAGS="$AM_LDFLAGS -L$with_nlopt/lib"])
 AC_ARG_WITH([LatCore],
     [AS_HELP_STRING([--with-LatCore=prefix],
     [use this option for a non-standard install prefix of the LatCore library])],
@@ -84,6 +90,13 @@ AC_CHECK_LIB([m],[cos],[],[AC_MSG_ERROR([libm library not found])])
 AC_CHECK_LIB([gslcblas],[cblas_dgemm],[],
              [AC_MSG_ERROR([GSL CBLAS library not found])])
 AC_CHECK_LIB([gsl],[gsl_blas_dgemm],[],[AC_MSG_ERROR([GSL library not found])])
+AC_CHECK_LIB([nlopt_cxx],[nlopt_create],
+    [AC_DEFINE([HAVE_NLOPT],
+    [1],
+    [Define to 1 if you have the `NLopt' library (-lnlopt_cxx).])]
+    [have_nlopt=true]
+    [LIBS="$LIBS -lnlopt_cxx"],[])
+AM_CONDITIONAL([HAVE_NLOPT], [test x$have_nlopt = xtrue])
 AC_CHECK_LIB([LatCore],[_ZN7LatCore12testFunctionEv],[],
     [AC_MSG_ERROR([LatCore library not found])])
 SAVED_LDFLAGS=$LDFLAGS

examples/exFit.cpp

@@ -3,6 +3,7 @@
 #include <LatAnalyze/CompiledModel.hpp>
 #include <LatAnalyze/Io.hpp>
 #include <LatAnalyze/MinuitMinimizer.hpp>
+#include <LatAnalyze/NloptMinimizer.hpp>
 #include <LatAnalyze/RandGen.hpp>
 #include <LatAnalyze/XYStatData.hpp>
 
@@ -37,25 +38,36 @@ int main(void)
             data.x(i2, 1)                  = xBuf[1];
             data.y(data.dataIndex(i1, i2)) = rg.gaussian(f(xBuf, exactPar),
                                                          yErr);
-            printf("% 8e % 8e % 8e % 8e % 8e\n", data.x(i1, 0), xErr,
-                   data.x(i2, 1), data.y(i1), yErr);
         }
     }
-    cout << endl;
     data.setXError(0, DVec::Constant(data.getXSize(0), xErr));
     data.assumeXExact(true, 1);
     data.setYError(0, DVec::Constant(data.getYSize(), yErr));
-    cout << data << endl;
+    
+    // set minimizers
+    DVec                init = DVec::Constant(2, 0.1);
+    FitResult           p;
+    NloptMinimizer      globalMin(NloptMinimizer::Algorithm::GN_CRS2_LM);
+    MinuitMinimizer     localMin;
+    vector<Minimizer *> min{&globalMin, &localMin};
+    
+    globalMin.setVerbosity(Minimizer::Verbosity::Normal);
+    globalMin.setPrecision(0.01);
+    globalMin.setMaxIteration(10000);
+    globalMin.useLowLimit(0);
+    globalMin.setLowLimit(0, 0.);
+    globalMin.useHighLimit(0);
+    globalMin.setHighLimit(0, 20.);
+    globalMin.useLowLimit(1);
+    globalMin.setLowLimit(1, 0.);
+    globalMin.useHighLimit(1);
+    globalMin.setHighLimit(1, 20.);
+    localMin.setVerbosity(Minimizer::Verbosity::Normal);
     
     // fit
-    DVec init = DVec::Constant(2, 0.1);
-    FitResult p;
-    MinuitMinimizer minimizer;
-    
     f.parName().setName(0, "m");
     f.parName().setName(1, "A");
-    minimizer.setVerbosity(Minimizer::Verbosity::Normal);
+    p = data.fit(min, init, f);
-    p = data.fit(minimizer, init, f);
     p.print();
     
     return EXIT_SUCCESS;

examples/exFitSample.cpp

@@ -2,6 +2,7 @@
 #include <cmath>
 #include <LatAnalyze/CompiledModel.hpp>
 #include <LatAnalyze/MinuitMinimizer.hpp>
+#include <LatAnalyze/NloptMinimizer.hpp>
 #include <LatAnalyze/RandGen.hpp>
 #include <LatAnalyze/XYSampleData.hpp>
 
@@ -43,16 +44,29 @@ int main(void)
         }
     }
     data.assumeXExact(true, 1);
-    cout << data << endl;
+    
+    // set minimizers
+    DVec                init = DVec::Constant(2, 0.1);
+    SampleFitResult     p;
+    NloptMinimizer      globalMin(NloptMinimizer::Algorithm::GN_CRS2_LM);
+    MinuitMinimizer     localMin;
+    vector<Minimizer *> min{&globalMin, &localMin};
+    
+    globalMin.setPrecision(0.01);
+    globalMin.setMaxIteration(10000);
+    globalMin.useLowLimit(0);
+    globalMin.setLowLimit(0, 0.);
+    globalMin.useHighLimit(0);
+    globalMin.setHighLimit(0, 20.);
+    globalMin.useLowLimit(1);
+    globalMin.setLowLimit(1, 0.);
+    globalMin.useHighLimit(1);
+    globalMin.setHighLimit(1, 20.);
     
     // fit
-    DVec init = DVec::Constant(2, 0.1);
-    SampleFitResult p;
-    MinuitMinimizer minimizer;
-
     f.parName().setName(0, "m");
     f.parName().setName(1, "A");
-    p = data.fit(minimizer, init, f);
+    p = data.fit(min, init, f);
     p.print();
     
     return EXIT_SUCCESS;

lib/Makefile.am

@@ -12,7 +12,7 @@ AM_YFLAGS = -d
 
 BUILT_SOURCES = AsciiParser.hpp MathParser.hpp
 
-lib_LTLIBRARIES    = libLatAnalyze.la
+lib_LTLIBRARIES = libLatAnalyze.la
 
 libLatAnalyze_la_SOURCES = \
     AsciiFile.cpp          \
@@ -86,6 +86,10 @@ if HAVE_MINUIT
     libLatAnalyze_la_SOURCES += MinuitMinimizer.cpp
     libLatAnalyze_la_HEADERS += MinuitMinimizer.hpp
 endif
+if HAVE_NLOPT
+    libLatAnalyze_la_SOURCES += NloptMinimizer.cpp
+    libLatAnalyze_la_HEADERS += NloptMinimizer.hpp
+endif
 libLatAnalyze_la_CXXFLAGS = $(COM_CXXFLAGS)
 
 ACLOCAL_AMFLAGS = -I .buildutils/m4

lib/NloptMinimizer.cpp

@@ -0,0 +1,175 @@
+/*
+ * NloptMinimizer.cpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2016 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/NloptMinimizer.hpp>
+#include <LatAnalyze/includes.hpp>
+
+using namespace std;
+using namespace Latan;
+
+/******************************************************************************
+ *                      NloptMinimizer implementation                         *
+ ******************************************************************************/
+// constructors ////////////////////////////////////////////////////////////////
+NloptMinimizer::NloptMinimizer(const Algorithm algorithm)
+{
+    setAlgorithm(algorithm);
+}
+
+NloptMinimizer::NloptMinimizer(const Index dim, const Algorithm algorithm)
+: Minimizer(dim)
+{
+    setAlgorithm(algorithm);
+}
+
+// access //////////////////////////////////////////////////////////////////////
+NloptMinimizer::Algorithm NloptMinimizer::getAlgorithm(void) const
+{
+    return algorithm_;
+}
+
+void NloptMinimizer::setAlgorithm(const Algorithm algorithm)
+{
+    algorithm_ = algorithm;
+}
+
+// minimization ////////////////////////////////////////////////////////////////
+const DVec & NloptMinimizer::operator()(const DoubleFunction &f)
+{
+    DVec &x = getState();
+    
+    // resize minimizer state to match function number of arguments
+    if (f.getNArg() != x.size())
+    {
+        resize(f.getNArg());
+    }
+    
+    // create and set minimizer
+    nlopt::opt     min(getAlgorithm(), x.size());
+    NloptFuncData  data;
+    vector<double> lb(x.size()), hb(x.size());
+    
+    min.set_maxeval(getMaxIteration());
+    min.set_xtol_rel(getPrecision());
+    data.f = &f;
+    min.set_min_objective(&funcWrapper, &data);
+    for (Index i = 0; i < x.size(); ++i)
+    {
+        lb[i] = hasLowLimit(i)  ? getLowLimit(i)  : -HUGE_VAL;
+        hb[i] = hasHighLimit(i) ? getHighLimit(i) :  HUGE_VAL;
+    }
+    min.set_lower_bounds(lb);
+    min.set_upper_bounds(hb);
+    
+    // minimise
+    double         res;
+    vector<double> vx(x.size());
+    nlopt::result  status;
+    unsigned int   n = 0;
+    
+    for (Index i = 0; i < x.size(); ++i)
+    {
+        vx[i] = x(i);
+    }
+    do
+    {
+        if (getVerbosity() >= Verbosity::Normal)
+        {
+            cout << "========== NLopt minimization, pass #" << n + 1;
+            cout << " ==========" << endl;
+            cout << "Algorithm: " << min.get_algorithm_name() << endl;
+            cout << "Max eval.= " << min.get_maxeval();
+            cout << " -- Precision= " << min.get_xtol_rel() << endl;
+            cout << "starting f(x)= " << f(x) << endl;
+        }
+        try
+        {
+            status = min.optimize(vx, res);
+        }
+        catch (invalid_argument &e)
+        {
+            LATAN_ERROR(Runtime, "NLopt has reported receving invalid "
+                        "arguments (if you are using a global minimizer, did "
+                        "you specify limits for all variables?)");
+        }
+        if (getVerbosity() >= Verbosity::Normal)
+        {
+            cout << "Found minimum " << res << " at:" << endl;
+            for (Index i = 0; i < x.size(); ++i)
+            {
+                printf("%8s= %e\n", f.varName().getName(i).c_str(), vx[i]);
+            }
+            cout << "after " << data.evalCount << " evaluations" << endl;
+            cout << "Minimization ended with code " << status;
+            cout << " (" << returnMessage(status) << ")";
+            cout << endl;
+        }
+        data.evalCount = 0;
+        for (Index i = 0; i < x.size(); ++i)
+        {
+            x(i) = vx[i];
+        }
+        n++;
+    } while ((status != nlopt::XTOL_REACHED) and (status != nlopt::SUCCESS)
+             and (n < getMaxPass()));
+    if (getVerbosity() >= Verbosity::Normal)
+    {
+        cout << "=================================================" << endl;
+    }
+    if ((status != nlopt::XTOL_REACHED) and (status != nlopt::SUCCESS))
+    {
+        LATAN_WARNING("invalid minimum: " + returnMessage(status));
+    }
+    
+    return x;
+}
+
+// NLopt return code parser ////////////////////////////////////////////////////
+string NloptMinimizer::returnMessage(const nlopt::result status)
+{
+    switch (status)
+    {
+        case nlopt::SUCCESS:
+            return "success";
+        case nlopt::STOPVAL_REACHED:
+            return "stopping value reached";
+        case nlopt::FTOL_REACHED:
+            return "tolerance on function reached";
+        case nlopt::XTOL_REACHED:
+            return "tolerance on variable reached";
+        case nlopt::MAXEVAL_REACHED:
+            return "maximum function evaluation reached";
+        case nlopt::MAXTIME_REACHED:
+            return "maximum time reached";
+        default:
+            return "";
+    }
+}
+
+// NLopt function wrapper //////////////////////////////////////////////////////
+double NloptMinimizer::funcWrapper(unsigned int n __dumb, const double *arg,
+                                   double *grad , void *vdata)
+{
+    NloptFuncData &data = *static_cast<NloptFuncData *>(vdata);
+    
+    assert(grad == nullptr);
+    data.evalCount++;
+    
+    return (*data.f)(arg);
+}

lib/NloptMinimizer.hpp

@@ -0,0 +1,73 @@
+/*
+ * NloptMinimizer.hpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2016 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/>.
+ */
+
+#ifndef Latan_NloptMinimizer_hpp_
+#define Latan_NloptMinimizer_hpp_
+
+#include <LatAnalyze/Global.hpp>
+#include <LatAnalyze/Derivative.hpp>
+#include <LatAnalyze/Function.hpp>
+#include <LatAnalyze/Minimizer.hpp>
+#include <nlopt.hpp>
+
+BEGIN_LATAN_NAMESPACE
+
+/******************************************************************************
+ *                      interface to NLOpt minimizers                         *
+ * ( http://ab-initio.mit.edu/wiki/index.php/NLopt )                          *
+ * -------------------------------------------------------------------------- *
+ * cf. http://ab-initio.mit.edu/wiki/index.php/NLopt_Algorithms for algorithm *
+ * references and naming conventions                                          *
+ ******************************************************************************/
+class NloptMinimizer: public Minimizer
+{
+public:
+    typedef nlopt::algorithm Algorithm;
+private:
+    struct NloptFuncData
+    {
+        const DoubleFunction *f{nullptr};
+        unsigned int         evalCount{0};
+    };
+public:
+    // constructor
+    NloptMinimizer(const Algorithm algorithm = defaultAlg_);
+    explicit NloptMinimizer(const Index dim,
+                            const Algorithm algorithm = defaultAlg_);
+    // destructor
+    virtual ~NloptMinimizer(void) = default;
+    // access
+    Algorithm getAlgorithm(void) const;
+    void      setAlgorithm(const Algorithm algorithm);
+    // minimization
+    virtual const DVec & operator()(const DoubleFunction &f);
+private:
+    // NLopt return code parser
+    static std::string returnMessage(const nlopt::result status);
+    // NLopt function wrapper
+    static double funcWrapper(unsigned int n, const double *arg,
+                              double *grad , void *vdata);
+private:
+    Algorithm                  algorithm_;
+    static constexpr Algorithm defaultAlg_ = Algorithm::LN_NELDERMEAD;
+};
+
+END_LATAN_NAMESPACE
+
+#endif // Latan_NloptMinimizer_hpp_