diff --git a/Code/Mantid/CurveFitting/inc/MantidCurveFitting/Fit.h b/Code/Mantid/CurveFitting/inc/MantidCurveFitting/Fit.h
index 4aabe93425168847b68a46619c278794295fd373..8c65a61c550c418961e81303fded68483252536b 100644
--- a/Code/Mantid/CurveFitting/inc/MantidCurveFitting/Fit.h
+++ b/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
diff --git a/Code/Mantid/CurveFitting/inc/MantidCurveFitting/GSLFunctions.h b/Code/Mantid/CurveFitting/inc/MantidCurveFitting/GSLFunctions.h
index 5fe4277e5ab0484bcebdd081320832693573b38e..d6ca1f30f38836942c2226c1a014d7de5a6c2aa4 100644
--- a/Code/Mantid/CurveFitting/inc/MantidCurveFitting/GSLFunctions.h
+++ b/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
diff --git a/Code/Mantid/CurveFitting/inc/MantidCurveFitting/LevenbergMarquardtMinimizer.h b/Code/Mantid/CurveFitting/inc/MantidCurveFitting/LevenbergMarquardtMinimizer.h
index e35ac554824fb970c47dd7201b87d8ba68124bbd..d78cd42a45b2404796d44a81eaee721fdda023f2 100644
--- a/Code/Mantid/CurveFitting/inc/MantidCurveFitting/LevenbergMarquardtMinimizer.h
+++ b/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;
};
diff --git a/Code/Mantid/CurveFitting/inc/MantidCurveFitting/SimplexMinimizer.h b/Code/Mantid/CurveFitting/inc/MantidCurveFitting/SimplexMinimizer.h
index 747285c43c31b11751403d1efccf4e83bba74fc5..49c4c5f3f6905505dc5134d38fda4a3f0bb29d76 100644
--- a/Code/Mantid/CurveFitting/inc/MantidCurveFitting/SimplexMinimizer.h
+++ b/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;
};
diff --git a/Code/Mantid/CurveFitting/src/Fit.cpp b/Code/Mantid/CurveFitting/src/Fit.cpp
index 6fbc54dd62f54280e07a7cb5db0e8502735c4425..6973c36d7beaa00242c771f3533bfae02725caca 100644
--- a/Code/Mantid/CurveFitting/src/Fit.cpp
+++ b/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
diff --git a/Code/Mantid/CurveFitting/src/GSLFunctions.cpp b/Code/Mantid/CurveFitting/src/GSLFunctions.cpp
index 9f152ada5c6772e3144a7e31782c66813793dbab..d0911b0e4067cd11e845b8bec91bee26e376d121 100644
--- a/Code/Mantid/CurveFitting/src/GSLFunctions.cpp
+++ b/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
diff --git a/Code/Mantid/CurveFitting/src/LevenbergMarquardtMinimizer.cpp b/Code/Mantid/CurveFitting/src/LevenbergMarquardtMinimizer.cpp
index 15971d5999f4d75706469ee5a022c1d867c04df4..cf0c798db10efd8de7048ac5e791e36e41c765d0 100644
--- a/Code/Mantid/CurveFitting/src/LevenbergMarquardtMinimizer.cpp
+++ b/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);
}