Prepare modification of interface to minimizer. Test this on SimplexMinimizer first. Refs #1254.

Code/Mantid/CurveFitting/inc/MantidCurveFitting/Fit.h

@@ -6,6 +6,8 @@
 //----------------------------------------------------------------------
 #include "MantidAPI/Algorithm.h"
 #include "MantidAPI/IFunction.h"
+#include "MantidCurveFitting/CostFuncLeastSquares.h"
+#include "MantidCurveFitting/CostFuncIgnorePosPeaks.h"
 
 namespace Mantid
 {
@@ -74,7 +76,7 @@ namespace Mantid
       /// Set a function for fitting
       void setFunction(API::IFunction* fun);
       /// Get the function for fitting
-      const API::IFunction* getFunction()const{return m_function;}
+      API::IFunction* getFunction()const{return m_function;}
 
     protected:
       // Overridden Algorithm methods

Code/Mantid/CurveFitting/inc/MantidCurveFitting/GSLFunctions.h

@@ -80,9 +80,9 @@ namespace Mantid
   };
 
   /// Structure to contain least squares data and used by GSL
-  struct FitData1 {
+  struct GSL_FitData {
     /// Constructor
-    FitData1(Fit* f);
+    GSL_FitData(Fit* f);
     /// number of points to be fitted (size of X, Y and sqrtWeightData arrays)
     size_t n;
     /// number of (active) fit parameters

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

@@ -7,6 +7,8 @@
 #include "MantidCurveFitting/IFuncMinimizer.h"
 #include <gsl/gsl_multifit_nlin.h>
 #include "MantidAPI/IFunction.h"
+#include "MantidCurveFitting/Fit.h"
+#include "MantidCurveFitting/GSLFunctions.h"
 
 namespace Mantid
 {
@@ -52,16 +54,27 @@ public:
   int hasConverged();
   double costFunctionVal();
   void calCovarianceMatrix(double epsrel, gsl_matrix * covar);
+  void initialize(double* X, const double* Y, double *sqrtWeight, const int& nData, const int& nParam, 
+    gsl_vector* startGuess, Fit* fit, const std::string& costFunction);
 
 private:
   /// name of this minimizer
   const std::string m_name;
 
+  /// GSL data container
+  GSL_FitData *m_data;
+
+  /// GSL minimizer container
+  gsl_multifit_function_fdf gslContainer;
+
   /// pointer to the GSL solver doing the work
   gsl_multifit_fdfsolver *m_gslSolver;
 
   /// Stored to access IFunction interface in iterate()
   API::IFunction* m_function;
+
+	/// Static reference to the logger class
+	static Kernel::Logger& g_log;
 };
 
 

Code/Mantid/CurveFitting/inc/MantidCurveFitting/SimplexMinimizer.h

@@ -5,6 +5,8 @@
 // Includes
 //----------------------------------------------------------------------
 #include "MantidCurveFitting/IFuncMinimizer.h"
+#include "MantidCurveFitting/Fit.h"
+#include "MantidCurveFitting/GSLFunctions.h"
 #include <gsl/gsl_multimin.h>
 
 namespace Mantid
@@ -42,7 +44,10 @@ class DLLExport SimplexMinimizer : public IFuncMinimizer
 public:
   /// constructor and destructor
   ~SimplexMinimizer();
-  SimplexMinimizer(gsl_multimin_function& gslContainer, gsl_vector* startGuess, const double& size);
+  SimplexMinimizer(): m_name("Simplex"), m_size(1.0) {}
+
+  void resetSize(double* X, const double* Y, double *sqrtWeight, const int& nData, const int& nParam, 
+    gsl_vector* startGuess, const double& size, Fit* fit, const std::string& costFunction);
 
   /// Overloading base class methods
   std::string name()const;
@@ -50,6 +55,8 @@ public:
   int hasConverged();
   double costFunctionVal();
   void calCovarianceMatrix(double epsrel, gsl_matrix * covar);
+  void initialize(double* X, const double* Y, double *sqrtWeight, const int& nData, const int& nParam, 
+    gsl_vector* startGuess, Fit* fit, const std::string& costFunction);
 
 private:
   /// name of this minimizer
@@ -58,8 +65,23 @@ private:
   /// pointer to the GSL solver doing the work
   gsl_multimin_fminimizer *m_gslSolver;
 
+  /// clear memory
+  void clearMemory();
+
+  /// size of simplex
+  double m_size;
+
   /// used by GSL
   gsl_vector *m_simplexStepSize;
+
+  /// GSL data container
+  GSL_FitData *m_data;
+
+  /// GSL simplex minimizer container
+  gsl_multimin_function gslContainer;
+
+	/// Static reference to the logger class
+	static Kernel::Logger& g_log;
 };
 
 

Code/Mantid/CurveFitting/src/Fit.cpp

@@ -12,8 +12,6 @@
 #include "MantidDataObjects/Workspace2D.h"
 #include "MantidKernel/UnitFactory.h"
 #include "MantidCurveFitting/ICostFunction.h"
-#include "MantidCurveFitting/CostFuncLeastSquares.h"
-#include "MantidCurveFitting/CostFuncIgnorePosPeaks.h"
 #include "MantidCurveFitting/BoundaryConstraint.h"
 #include "MantidCurveFitting/LevenbergMarquardtMinimizer.h"
 #include "MantidCurveFitting/SimplexMinimizer.h"
@@ -211,7 +209,7 @@ namespace CurveFitting
     // since as a rule of thumb this is required as a minimum to obtained 'accurate'
     // fitting parameter values.
 
-    FitData1 l_data(this);
+    GSL_FitData l_data(this);
 
     l_data.p = m_function->nActive();
     l_data.n = m_maxX - m_minX; // m_minX and m_maxX are array index markers. I.e. e.g. 0 & 19.
@@ -282,10 +280,10 @@ namespace CurveFitting
 
     // set-up GSL container to be used with GSL simplex algorithm
 
-    gsl_multimin_function gslSimplexContainer;
+   /* gsl_multimin_function gslSimplexContainer;
     gslSimplexContainer.n = l_data.p;  // n here refers to number of parameters
     gslSimplexContainer.f = &gsl_costFunction;
-    gslSimplexContainer.params = &l_data;
+    gslSimplexContainer.params = &l_data;*/
 
 
     // set-up GSL container to be used with none-least squares GSL routines using derivatives
@@ -341,7 +339,11 @@ namespace CurveFitting
 
     if ( methodUsed.compare("Simplex") == 0 )
     {
-      minimizer = new SimplexMinimizer(gslSimplexContainer, initFuncArg, 1.0);
+      //minimizer = new SimplexMinimizer(gslSimplexContainer, initFuncArg, 1.0);
+      SimplexMinimizer* sm = new SimplexMinimizer;
+      sm->initialize(l_data.X, l_data.Y, l_data.sqrtWeightData, l_data.n, l_data.p, 
+                     initFuncArg, this, costFunction);
+      minimizer = sm;
     }
     else
     {
@@ -392,7 +394,10 @@ namespace CurveFitting
             methodUsed = "Simplex";
             //simplexFallBack = true;
             delete minimizer;
-            minimizer = new SimplexMinimizer(gslSimplexContainer, initFuncArg, 1.0);
+            SimplexMinimizer* sm = new SimplexMinimizer;
+            sm->initialize(l_data.X, l_data.Y, l_data.sqrtWeightData, l_data.n, l_data.p, 
+                           initFuncArg, this, costFunction);
+            minimizer = sm;
             iter = 0;
             g_log.warning() << "Fit algorithm using Levenberg-Marquardt failed "
               << "reporting the following: " << gsl_strerror(status) << "\n"
@@ -423,8 +428,14 @@ namespace CurveFitting
           { 
             g_log.information() << "Simplex step size reduced to 0.1\n";
             delete minimizer;
-            minimizer = new SimplexMinimizer(gslSimplexContainer, initFuncArg, 0.1);
-            //iter = 0;
+            SimplexMinimizer* sm = new SimplexMinimizer;
+            sm->initialize(l_data.X, l_data.Y, l_data.sqrtWeightData, l_data.n, l_data.p, 
+                           initFuncArg, this, costFunction);
+            sm->resetSize(l_data.X, l_data.Y, l_data.sqrtWeightData, l_data.n, l_data.p, 
+                           initFuncArg, 0.1, this, costFunction);
+            minimizer = sm;
+                  //sm = static_cast<SimplexMinimizer*>(minimizer);
+
             status = GSL_CONTINUE;
             continue;
           }
@@ -781,24 +792,5 @@ namespace CurveFitting
     return deriv;
   }
 
-  /**
-   * Constructor. Creates declared -> active index map
-   * @param f Pointer to the Fit algorithm
-   */
-  FitData1::FitData1(Fit* f):fit(f)
-  {
-    int j = 0;
-    for(int i=0;i<f->m_function->nParams();++i)
-    {
-      if (f->m_function->isActive(i))
-      {
-        J.m_index.push_back(j);
-        j++;
-      }
-      else
-        J.m_index.push_back(-1);
-    }
-  }
-
 } // namespace Algorithm
 } // namespace Mantid

Code/Mantid/CurveFitting/src/GSLFunctions.cpp

@@ -4,6 +4,7 @@
 #include "MantidCurveFitting/GSLFunctions.h"
 #include "MantidCurveFitting/Fit.h"
 #include "MantidCurveFitting/ICostFunction.h"
+#include "MantidAPI/CompositeFunction.h"
 
 
 namespace Mantid
@@ -21,7 +22,7 @@ namespace CurveFitting
   */
   int gsl_f(const gsl_vector * x, void *params, gsl_vector * f) {
 
-    struct FitData1 *p = (struct FitData1 *)params;
+    struct GSL_FitData *p = (struct GSL_FitData *)params;
     p->fit->function (x->data, f->data, p->X, p->n);
 
     // function() return calculated data values. Need to convert this values into
@@ -42,7 +43,7 @@ namespace CurveFitting
   */
   int gsl_df(const gsl_vector * x, void *params, gsl_matrix * J) {
 
-    struct FitData1 *p = (struct FitData1 *)params;
+    struct GSL_FitData *p = (struct GSL_FitData *)params;
 
     p->J.setJ(J);
 
@@ -83,7 +84,7 @@ namespace CurveFitting
   double gsl_costFunction(const gsl_vector * x, void *params)
   {
 
-    struct FitData1 *p = (struct FitData1 *)params;
+    struct GSL_FitData *p = (struct GSL_FitData *)params;
     double * l_holdCalculatedData = p->holdCalculatedData;
 
     // calculate yCal and store in l_holdCalculatedData
@@ -101,7 +102,7 @@ namespace CurveFitting
   void gsl_costFunction_df(const gsl_vector * x, void *params, gsl_vector *df)
   {
 
-    struct FitData1 *p = (struct FitData1 *)params;
+    struct GSL_FitData *p = (struct GSL_FitData *)params;
     double * l_holdCalculatedData = p->holdCalculatedData;
 
     p->fit->function (x->data, l_holdCalculatedData, p->X, p->n);
@@ -127,6 +128,24 @@ namespace CurveFitting
     gsl_costFunction_df(x, params, df); 
   }
 
+  /**
+   * Constructor. Creates declared -> active index map
+   * @param f Pointer to the Fit algorithm
+   */
+  GSL_FitData::GSL_FitData(Fit* f):fit(f)
+  {
+    int j = 0;
+    for(int i=0;i<f->getFunction()->nParams();++i)
+    {
+      if (f->getFunction()->isActive(i))
+      {
+        J.m_index.push_back(j);
+        j++;
+      }
+      else
+        J.m_index.push_back(-1);
+    }
+  }
 
 } // namespace CurveFitting
 } // namespace Mantid

Code/Mantid/CurveFitting/src/LevenbergMarquardtMinimizer.cpp

@@ -11,6 +11,53 @@ namespace Mantid
 namespace CurveFitting
 {
 
+// Get a reference to the logger
+Kernel::Logger& LevenbergMarquardtMinimizer::g_log = Kernel::Logger::get("LevenbergMarquardtMinimizer");
+
+void LevenbergMarquardtMinimizer::initialize(double* X, const double* Y, 
+                                             double *sqrtWeight, const int& nData, 
+                                             const int& nParam, gsl_vector* startGuess, 
+                                             Fit* fit, const std::string& costFunction)
+{
+  // set-up GSL container to be used with GSL simplex algorithm
+  m_data = new GSL_FitData(fit);  //,X, Y, sqrtWeight, nData, nParam);
+  m_data->p = nParam;
+  m_data->n = nData; 
+  m_data->X = X;
+  m_data->Y = Y;
+  m_data->sqrtWeightData = sqrtWeight;
+  m_data->holdCalculatedData = new double[nData];
+  m_data->holdCalculatedJacobian =  gsl_matrix_alloc (nData, nParam);
+
+  if ( costFunction.compare("Least squares") == 0 )
+    m_data->costFunc = new CostFuncLeastSquares();
+  else if ( costFunction.compare("Ignore positive peaks") == 0 )
+    m_data->costFunc = new CostFuncIgnorePosPeaks();
+  else
+  {
+    g_log.error("Unrecognised cost function. Default to Least squares\n");
+    m_data->costFunc = new CostFuncLeastSquares();
+  }
+
+  // specify the type of GSL solver to use
+  const gsl_multifit_fdfsolver_type *T = gsl_multifit_fdfsolver_lmsder;
+
+    gslContainer.f = &gsl_f;
+    gslContainer.df = &gsl_df;
+    gslContainer.fdf = &gsl_fdf;
+    gslContainer.n = nData;
+    gslContainer.p = nParam;
+    gslContainer.params = m_data;
+
+
+  // setup GSL solver
+  m_gslSolver = gsl_multifit_fdfsolver_alloc(T, nData, nParam);
+  gsl_multifit_fdfsolver_set(m_gslSolver, &gslContainer, startGuess);
+
+  m_function = fit->getFunction();
+
+}
+
 LevenbergMarquardtMinimizer::LevenbergMarquardtMinimizer(
   gsl_multifit_function_fdf& gslContainer, 
   gsl_vector* startGuess, API::IFunction* func) : m_name("Levenberg Marquardt") 
@@ -27,6 +74,11 @@ LevenbergMarquardtMinimizer::LevenbergMarquardtMinimizer(
 
 LevenbergMarquardtMinimizer::~LevenbergMarquardtMinimizer()
 {
+  //delete [] m_data->holdCalculatedData;
+  //delete m_data->costFunc;
+  //gsl_matrix_free (m_data->holdCalculatedJacobian);
+  //delete m_data;
+
   gsl_multifit_fdfsolver_free(m_gslSolver);
 }