code cleaning + derivative support for NLopt

lib/Derivative.cpp

@@ -31,11 +31,12 @@ using namespace Math;
 // constructor /////////////////////////////////////////////////////////////////
 Derivative::Derivative(const DoubleFunction &f, const Index dir,
                        const double step)
-: f_(f)
-, dir_(dir)
-, step_(step)
-, buffer_(new DVec(f.getNArg()))
-{}
+: buffer_(new DVec(f.getNArg()))
+{
+    setFunction(f);
+    setDir(dir);
+    setStep(step);
+}
 
 Derivative::Derivative(const DoubleFunction &f, const Index dir,
                        const Index order, const DVec &point, const double step)
@@ -45,6 +46,11 @@ Derivative::Derivative(const DoubleFunction &f, const Index dir,
 }
 
 // access //////////////////////////////////////////////////////////////////////
+Index Derivative::getDir(void) const
+{
+    return dir_;
+}
+
 Index Derivative::getOrder(void) const
 {
     return order_;
@@ -60,6 +66,11 @@ double Derivative::getStep(void) const
     return step_;
 }
 
+void Derivative::setDir(const Index dir)
+{
+    dir_ = dir;
+}
+
 void Derivative::setFunction(const DoubleFunction &f)
 {
     f_ = f;

lib/Derivative.hpp

@@ -39,9 +39,11 @@ public:
     // destructor
     virtual ~Derivative(void) = default;
     // access
+    Index  getDir(void) const;
     Index  getNPoint(void) const;
     Index  getOrder(void) const;
     double getStep(void) const;
+    void   setDir(const Index dir);
     void   setFunction(const DoubleFunction &f);
     void   setOrderAndPoint(const Index order, const DVec &point);
     void   setStep(const double step);
@@ -73,16 +75,13 @@ public:
     static const Index defaultPrecOrder = 2;
 public:
     // constructor
-    CentralDerivative(const DoubleFunction &f, const Index dir = 0,
+    CentralDerivative(const DoubleFunction &f = DoubleFunction(),
+                      const Index dir = 0,
                       const Index order = 1,
                       const Index precOrder = defaultPrecOrder);
     // destructor
     virtual ~CentralDerivative(void) = default;
     // access
-    using Derivative::getNPoint;
-    using Derivative::getStep;
-    using Derivative::getOrder;
-    using Derivative::setStep;
     Index getPrecOrder(void) const;
     void  setOrder(const Index order, const Index precOrder = defaultPrecOrder);
     // function call

lib/Minimizer.cpp

@@ -23,12 +23,6 @@
 using namespace std;
 using namespace Latan;
 
-// constructor /////////////////////////////////////////////////////////////////
-Minimizer::Minimizer(const Index dim)
-{
-    resize(dim);
-}
-
 // access //////////////////////////////////////////////////////////////////////
 void Minimizer::resize(const Index dim)
 {

lib/Minimizer.hpp

@@ -36,7 +36,6 @@ class Minimizer: public Solver
 public:
     // constructor
     Minimizer(void) = default;
-    explicit Minimizer(const Index dim);
     // destructor
     virtual ~Minimizer(void) = default;
     // access

lib/MinuitMinimizer.cpp

@@ -36,12 +36,6 @@ MinuitMinimizer::MinuitMinimizer(const Algorithm algorithm)
     setAlgorithm(algorithm);
 }
 
-MinuitMinimizer::MinuitMinimizer(const Index dim, const Algorithm algorithm)
-: Minimizer(dim)
-{
-    setAlgorithm(algorithm);
-}
-
 // access //////////////////////////////////////////////////////////////////////
 MinuitMinimizer::Algorithm MinuitMinimizer::getAlgorithm(void) const
 {

lib/MinuitMinimizer.hpp

@@ -41,9 +41,7 @@ public:
     };
 public:
     // constructor
-    MinuitMinimizer(const Algorithm algorithm = defaultAlg_);
-    explicit MinuitMinimizer(const Index dim,
-                             const Algorithm algorithm = defaultAlg_);
+    explicit MinuitMinimizer(const Algorithm algorithm = defaultAlg_);
     // destructor
     virtual ~MinuitMinimizer(void) = default;
     // access

lib/NloptMinimizer.cpp

@@ -30,12 +30,7 @@ using namespace Latan;
 NloptMinimizer::NloptMinimizer(const Algorithm algorithm)
 {
     setAlgorithm(algorithm);
-}
-
-NloptMinimizer::NloptMinimizer(const Index dim, const Algorithm algorithm)
-: Minimizer(dim)
-{
-    setAlgorithm(algorithm);
+    der_.setOrder(1, 1);
 }
 
 // access //////////////////////////////////////////////////////////////////////
@@ -67,7 +62,10 @@ const DVec & NloptMinimizer::operator()(const DoubleFunction &f)
     
     min.set_maxeval(getMaxIteration());
     min.set_xtol_rel(getPrecision());
+    min.set_ftol_rel(-1.);
+    der_.setFunction(f);
     data.f = &f;
+    data.d = &der_;
     min.set_min_objective(&funcWrapper, &data);
     for (Index i = 0; i < x.size(); ++i)
     {
@@ -126,13 +124,12 @@ const DVec & NloptMinimizer::operator()(const DoubleFunction &f)
             x(i) = vx[i];
         }
         n++;
-    } while ((status != nlopt::XTOL_REACHED) and (status != nlopt::SUCCESS)
-             and (n < getMaxPass()));
+    } while (!minSuccess(status) and (n < getMaxPass()));
     if (getVerbosity() >= Verbosity::Normal)
     {
         cout << "=================================================" << endl;
     }
-    if ((status != nlopt::XTOL_REACHED) and (status != nlopt::SUCCESS))
+    if (!minSuccess(status))
     {
         LATAN_WARNING("invalid minimum: " + returnMessage(status));
     }
@@ -163,13 +160,37 @@ string NloptMinimizer::returnMessage(const nlopt::result status)
 }
 
 // NLopt function wrapper //////////////////////////////////////////////////////
-double NloptMinimizer::funcWrapper(unsigned int n __dumb, const double *arg,
+double NloptMinimizer::funcWrapper(unsigned int n, const double *arg,
                                    double *grad , void *vdata)
 {
     NloptFuncData &data = *static_cast<NloptFuncData *>(vdata);
     
-    assert(grad == nullptr);
+    if (grad)
+    {
+        for (unsigned int i = 0; i < n; ++i)
+        {
+            data.d->setDir(i);
+            grad[i] = (*(data.d))(arg);
+        }
+        data.evalCount += data.d->getNPoint()*n;
+    }
     data.evalCount++;
     
     return (*data.f)(arg);
 }
+
+// NLopt return status parser //////////////////////////////////////////////////
+bool NloptMinimizer::minSuccess(const nlopt::result status)
+{
+    switch (status)
+    {
+        case nlopt::SUCCESS:
+        case nlopt::FTOL_REACHED:
+        case nlopt::XTOL_REACHED:
+            return true;
+            break;
+        default:
+            return false;
+            break;
+    }
+}

lib/NloptMinimizer.hpp

@@ -43,13 +43,12 @@ private:
     struct NloptFuncData
     {
         const DoubleFunction *f{nullptr};
+        Derivative           *d{nullptr};
         unsigned int         evalCount{0};
     };
 public:
     // constructor
-    NloptMinimizer(const Algorithm algorithm = defaultAlg_);
-    explicit NloptMinimizer(const Index dim,
-                            const Algorithm algorithm = defaultAlg_);
+    explicit NloptMinimizer(const Algorithm algorithm = defaultAlg_);
     // destructor
     virtual ~NloptMinimizer(void) = default;
     // access
@@ -63,9 +62,12 @@ private:
     // NLopt function wrapper
     static double funcWrapper(unsigned int n, const double *arg,
                               double *grad , void *vdata);
+    // NLopt return status parser
+    static bool minSuccess(const nlopt::result status);
 private:
     Algorithm                  algorithm_;
     static constexpr Algorithm defaultAlg_ = Algorithm::LN_NELDERMEAD;
+    CentralDerivative          der_;
 };
 
 END_LATAN_NAMESPACE