Re #4158. Improvements to new Leveberg-Marquardt.

Code/Mantid/Framework/API/inc/MantidAPI/FunctionFactory.h

@@ -134,8 +134,8 @@ namespace API
     for( std::vector<std::string>::const_iterator it = names.begin(); 
          it != names.end(); ++it )
     {
-      IFunction *func = this->createFitFunction(*it);
-      if( dynamic_cast<FunctionType*>(func) )
+      IFunction_sptr func = this->createFitFunction(*it);
+      if ( func && dynamic_cast<FunctionType*>(func.get()) )
       {
         typeNames.push_back(*it);
       }

Code/Mantid/Framework/API/inc/MantidAPI/IConstraint.h

@@ -60,6 +60,9 @@ class MANTID_API_DLL IConstraint: public ParameterReference
   /// Returns the derivative of the penalty for each active parameter
   virtual double checkDeriv() = 0;
 
+  /// Returns the derivative of the penalty for each active parameter
+  virtual double checkDeriv2() = 0;
+
   /// Set the parameters of IFitFunction to satisfy constraint. For example
   /// for a BoundaryConstraint this if param value less than lower boundary
   /// it is set to that value and vice versa for if the param value is larger

Code/Mantid/Framework/CurveFitting/CMakeLists.txt

@@ -139,15 +139,12 @@ set ( TEST_FILES
 	test/CompositeFunctionTest.h
 #	test/ConvolutionTest.h
 	test/ExpDecayTest.h
-#	test/ExpDecayMuonTest.h
-#	test/ExpDecayOscTest.h
-#	test/FitTest.h
+	test/ExpDecayMuonTest.h
+	test/ExpDecayOscTest.h
 	test/FitMWTest.h
 	test/FRConjugateGradientTest.h
-#	test/NewFitTest.h
 #	test/FuncMinimizerFactoryTest.h
-#	test/FunctionTest.h
-#	test/FunctionFactoryTest.h
+	test/FunctionFactoryTest.h
 #	test/Gaussian1DTest.h
 	test/GausDecayTest.h
 	test/GaussianTest.h
@@ -169,9 +166,9 @@ set ( TEST_FILES
 	test/SimplexTest.h
 #	test/SpecialFunctionSupportTest.h
 #	test/SplineBackgroundTest.h
-#	test/StretchExpTest.h
+	test/StretchExpTest.h
 #	test/UserFunction1DTest.h
-#	test/UserFunctionTest.h
+	test/UserFunctionTest.h
 )
 
 

Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/BoundaryConstraint.h

@@ -116,6 +116,7 @@ namespace Mantid
       /// overwrite IConstraint base class methods
       virtual double check();
       virtual double checkDeriv();
+      virtual double checkDeriv2();
       virtual void setParamToSatisfyConstraint();
       virtual std::string asString()const;
 

Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/LevenbergMarquardtMinimizer.h

@@ -75,6 +75,7 @@ class DLLExport LevenbergMarquardtMinimizer : public IFuncMinimizer
   GSLVector m_der;
   /// To keep second derivatives
   GSLMatrix m_hessian;
+  double m_oldDder;
 	/// Static reference to the logger class
 	static Kernel::Logger& g_log;
 };

Code/Mantid/Framework/CurveFitting/src/BoundaryConstraint.cpp

@@ -197,10 +197,16 @@ double BoundaryConstraint::check()
 
   if (m_hasLowerBound)
     if ( paramValue < m_lowerBound )
-      penalty = (m_lowerBound-paramValue)* m_penaltyFactor;
+    {
+      double dp = m_lowerBound-paramValue;
+      penalty = m_penaltyFactor * dp * dp;
+    }
   if (m_hasUpperBound)
     if ( paramValue > m_upperBound )
-      penalty = (paramValue-m_upperBound)* m_penaltyFactor;
+    {
+      double dp = paramValue-m_upperBound;
+      penalty = m_penaltyFactor * dp * dp;
+    }
 
   return penalty;
 }
@@ -220,10 +226,39 @@ double BoundaryConstraint::checkDeriv()
 
   if (m_hasLowerBound)
     if ( paramValue < m_lowerBound )
-      penalty = -m_penaltyFactor;
+    {
+      double dp = m_lowerBound-paramValue;
+      penalty = 2 * m_penaltyFactor * dp;
+    }
+  if (m_hasUpperBound)
+    if ( paramValue > m_upperBound )
+    {
+      double dp = paramValue-m_upperBound;
+      penalty = 2 * m_penaltyFactor * dp;
+    }
+
+  return penalty;
+}
+
+double BoundaryConstraint::checkDeriv2()
+{
+  double penalty = 0.0;
+
+  if (/*m_activeParameterIndex < 0 ||*/ !(m_hasLowerBound || m_hasUpperBound))
+  {
+    // no point in logging any warning here since checkDeriv() will always be called after
+    // check() is called 
+    return penalty;
+  } 
+
+  double paramValue = getFunction()->getParameter(getIndex());
+
+  if (m_hasLowerBound)
+    if ( paramValue < m_lowerBound )
+      penalty = 2 * m_penaltyFactor;
   if (m_hasUpperBound)
     if ( paramValue > m_upperBound )
-      penalty = m_penaltyFactor;
+      penalty = 2 * m_penaltyFactor;
 
   return penalty;
 }

Code/Mantid/Framework/CurveFitting/src/CostFuncLeastSquares.cpp

@@ -233,6 +233,14 @@ double CostFuncLeastSquares::valDerivHessian(GSLVector& der, GSLMatrix& hessian,
       {
         d += jacobian.get(k,i) * jacobian.get(k,j);
       }
+      if (i == j)
+      {
+        API::IConstraint* c = m_function->getConstraint(i);
+        if (c)
+        {
+          d += c->checkDeriv2();
+        }
+      }
       hessian.set(i1,i2,d);
       //std::cerr << "hess " << i1 << ' ' << i2 << std::endl;
       if (i1 != i2)

Code/Mantid/Framework/CurveFitting/src/LevenbergMarquardtMinimizer.cpp

@@ -43,6 +43,7 @@ void LevenbergMarquardtMinimizer::initialize(API::ICostFunction_sptr function)
   m_mu = 0;
   m_nu = 2.0;
   m_rho = 0;
+  m_oldDder = 0.0;
 }
 
 /// Do one iteration.
@@ -75,20 +76,7 @@ bool LevenbergMarquardtMinimizer::iterate()
 
   // Calculate damping to hessian
   if (m_mu == 0) // first iteration
-  {// initial m_mu = m_tau * max( hessian[i,i] )
-    //double maxH = 0.0;
-    //for(size_t i = 0; i < n; ++i)
-    //{
-    //  double tmp = m_hessian.get(i,i);
-    //  if (tmp > maxH)
-    //  {
-    //    maxH = tmp;
-    //  }
-    //}
-    //if (maxH <= 0)
-    //{
-    //  throw std::runtime_error("Hessian in Levenberg-Marquardt must be positively defined");
-    //}
+  {
     m_mu = m_tau;// * maxH;
     m_nu = 2.0;
   }
@@ -99,35 +87,14 @@ bool LevenbergMarquardtMinimizer::iterate()
   }
   // copy the hessian
   GSLMatrix H(m_hessian);
-  GSLMatrix I(n,n);
-  I.identity();
-  I *= m_mu;
-  H += I;
-
-  if (debug)
+  for(size_t i = 0; i < n; ++i)
   {
-    std::cerr << "hessian:" << std::endl;
-    for(size_t i = 0; i < n; ++i)
-    {
-      for(size_t j = 0; j < n; ++j)
-      {
-        std::cerr << m_hessian.get(i,j) << ' ';
-      }
-      std::cerr << std::endl;
-    }
-    std::cerr << "\n";
-
-    std::cerr << "I:" << std::endl;
-    for(size_t i = 0; i < n; ++i)
-    {
-      for(size_t j = 0; j < n; ++j)
-      {
-        std::cerr << I.get(i,j) << ' ';
-      }
-      std::cerr << std::endl;
-    }
-    std::cerr << "\n";
+    double d = H.get(i,i) + m_mu * fabs(m_der.get(i));
+    H.set(i,i,d);
+  }
 
+  if (debug && m_rho > 0)
+  {
     std::cerr << "H:" << std::endl;
     for(size_t i = 0; i < n; ++i)
     {
@@ -155,15 +122,18 @@ bool LevenbergMarquardtMinimizer::iterate()
   {
     for(size_t j = 0; j < n; ++j)  {std::cerr << m_der.get(j) << ' '; } std::cerr << std::endl;
     for(size_t j = 0; j < n; ++j)  {std::cerr << dx.get(j) << ' '; } std::cerr << std::endl << std::endl;
-    system("pause");
+    //system("pause");
   }
 
   // Update the parameters of the cost function.
   for(size_t i = 0; i < n; ++i)
   {
     double d = m_leastSquares->getParameter(i) + dx.get(i);
     m_leastSquares->setParameter(i,d);
-    //std::cerr << "P" << i << ' ' << d << std::endl;
+    if (debug)
+    {
+      std::cerr << "P" << i << ' ' << d << std::endl;
+    }
   }
   m_leastSquares->getFittingFunction()->applyTies();
 
@@ -188,27 +158,42 @@ bool LevenbergMarquardtMinimizer::iterate()
     std::cerr << "F " << m_F << ' ' << F1 << ' ' << dL << std::endl;
   }
   // Try the stop condition
-  if (fabs(dL) < m_relTol)
+  //if (fabs(dL) < m_relTol)
+  if (m_oldDder > 0.0 && fabs(m_oldDder - dder)/m_oldDder < 0.001 || dder < 0.001)
   {
+    if (debug)
+    std::cerr << "stopped at " << dder << ' ' << fabs(m_oldDder - dder)/m_oldDder << std::endl;
     return false;
   }
+
+  m_oldDder = dder;
+
+  if (fabs(dL) == 0.0) m_rho = 0;
+  else
+    m_rho = (m_F - F1) / dL;
+  if (debug)
+  {
+    std::cerr << "rho=" << m_rho << std::endl;
+  }
 
-  m_rho = (m_F - F1) / dL;
   if (m_rho > 0)
   {// good progress, decrease m_mu but no more than by 1/3
     // rho = 1 - (2*rho - 1)^3
     m_rho = 2.0 * m_rho - 1.0;
     m_rho = 1.0 - m_rho * m_rho * m_rho;
     const double I3 = 1.0 / 3.0;
-    if (I3 > m_rho) m_rho = I3;
+    if (m_rho > I3) m_rho = I3;
+    if (m_rho < 0.0001) m_rho = 0.1;
     m_mu *= m_rho;
     m_nu = 2.0;
     m_F = F1;
+    if (debug)
+    std::cerr << "times " << m_rho << std::endl;
   }
   else
   {// bad iteration. increase m_mu and revert changes to parameters
     m_mu *= m_nu;
-    m_nu *= 2.0;
+    //m_nu *= 2.0;
     // undo parameter update
     for(size_t i = 0; i < n; ++i)
     {
@@ -220,10 +205,6 @@ bool LevenbergMarquardtMinimizer::iterate()
     m_F = m_leastSquares->val();
   }
 
-  if (debug)
-  {
-    std::cerr << "rho=" << m_rho << std::endl;
-  }
   return true;
 }