diff --git a/Code/Mantid/CurveFitting/CurveFitting.vcproj b/Code/Mantid/CurveFitting/CurveFitting.vcproj index 7c01024270765e62ceb0529891ee8bec658743c2..00762c61dbedbaeb86f65e2f4df03afd11fda185 100644 --- a/Code/Mantid/CurveFitting/CurveFitting.vcproj +++ b/Code/Mantid/CurveFitting/CurveFitting.vcproj @@ -258,6 +258,10 @@ RelativePath=".\src\Quadratic.cpp" > </File> + <File + RelativePath=".\src\SimplexMinimizer.cpp" + > + </File> <File RelativePath=".\src\SpecialFunctionHelper.cpp" > @@ -400,6 +404,10 @@ RelativePath=".\test\QuadraticTest.h" > </File> + <File + RelativePath=".\inc\MantidCurveFitting\SimplexMinimizer.h" + > + </File> <File RelativePath=".\inc\MantidCurveFitting\SpecialFunctionSupport.h" > diff --git a/Code/Mantid/CurveFitting/inc/MantidCurveFitting/SimplexMinimizer.h b/Code/Mantid/CurveFitting/inc/MantidCurveFitting/SimplexMinimizer.h new file mode 100644 index 0000000000000000000000000000000000000000..a29b39e1aa19a6dfc1c1ada9bc8c4031d55a1e5d --- /dev/null +++ b/Code/Mantid/CurveFitting/inc/MantidCurveFitting/SimplexMinimizer.h @@ -0,0 +1,68 @@ +#ifndef MANTID_CURVEFITTING_SIMPLEXMINIMIZER_H_ +#define MANTID_CURVEFITTING_SIMPLEXMINIMIZER_H_ + +//---------------------------------------------------------------------- +// Includes +//---------------------------------------------------------------------- +#include "MantidCurveFitting/IFuncMinimizer.h" +#include <gsl/gsl_multimin.h> + +namespace Mantid +{ +namespace CurveFitting +{ +/** Implementing Simplex by wrapping the IFuncMinimizer interface + around the GSL implementation of this algorithm. + + @author Anders Markvardsen, ISIS, RAL + @date 8/1/2010 + + Copyright © 2009 STFC Rutherford Appleton Laboratory + + This file is part of Mantid. + + Mantid 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. + + Mantid 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 this program. If not, see <http://www.gnu.org/licenses/>. + + File change history is stored at: <https://svn.mantidproject.org/mantid/trunk/Code/Mantid>. + Code Documentation is available at: <http://doxygen.mantidproject.org> +*/ +class DLLExport SimplexMinimizer : public IFuncMinimizer +{ +public: + /// constructor and destructor + ~SimplexMinimizer(); + SimplexMinimizer(gsl_multimin_function& gslContainer, gsl_vector* startGuess, const double& size); + + /// Overloading base class methods + std::string name()const; + int iterate(); + int hasConverged(); + double costFunctionVal(); + +private: + /// name of this minimizer + const std::string m_name; + + /// pointer to the GSL solver doing the work + gsl_multimin_fminimizer *m_gslSolver; + + /// used by GSL + gsl_vector *m_simplexStepSize; +}; + + +} // namespace CurveFitting +} // namespace Mantid + +#endif /*MANTID_CURVEFITTING_SIMPLEXMINIMIZER_H_*/ diff --git a/Code/Mantid/CurveFitting/src/Fit.cpp b/Code/Mantid/CurveFitting/src/Fit.cpp index 618e816817f212141fdb2e42fb873803fbc298c4..06cb394eb486b09026fc9ac89a1b96195c6a5932 100644 --- a/Code/Mantid/CurveFitting/src/Fit.cpp +++ b/Code/Mantid/CurveFitting/src/Fit.cpp @@ -13,6 +13,7 @@ #include "MantidKernel/UnitFactory.h" #include "MantidCurveFitting/BoundaryConstraint.h" #include "MantidCurveFitting/LevenbergMarquardtMinimizer.h" +#include "MantidCurveFitting/SimplexMinimizer.h" #include <gsl/gsl_statistics.h> #include <gsl/gsl_multifit_nlin.h> @@ -430,26 +431,18 @@ namespace CurveFitting f.p = l_data.p; f.params = &l_data; - // set-up remaining GSL machinery for least squared + // set-up minimizer - IFuncMinimizer* minimizer; + IFuncMinimizer* minimizer = NULL; if (isDerivDefined) { minimizer = new LevenbergMarquardtMinimizer(f, initFuncArg); } - - // set-up remaining GSL machinery to use simplex algorithm - // always set this algorithm up since in case levenberg-marquardt fail - // then simplex is used as fall-back algorithm - - const gsl_multimin_fminimizer_type *simplexType = gsl_multimin_fminimizer_nmsimplex; - gsl_multimin_fminimizer *simplexMinimizer = NULL; - gsl_vector *simplexStepSize = NULL; - simplexMinimizer = gsl_multimin_fminimizer_alloc(simplexType, l_data.p); - simplexStepSize = gsl_vector_alloc(l_data.p); - gsl_vector_set_all (simplexStepSize, 1); // is this always a sensible starting step size? - gsl_multimin_fminimizer_set(simplexMinimizer, &gslSimplexContainer, initFuncArg, simplexStepSize); + else + { + minimizer = new SimplexMinimizer(gslSimplexContainer, initFuncArg, 1.0); + } // finally do the fitting @@ -457,7 +450,6 @@ namespace CurveFitting int iter = 0; int status; bool simplexFallBack = false; // set to true if levenberg-marquardt fails - double size; // for simplex algorithm double finalCostFuncVal; double dof = l_data.n - l_data.p; // dof stands for degrees of freedom @@ -481,6 +473,8 @@ namespace CurveFitting if (iter < 3) { simplexFallBack = true; + delete minimizer; + minimizer = new SimplexMinimizer(gslSimplexContainer, initFuncArg, 1.0); iter = 0; g_log.warning() << "Fit algorithm using Levenberg-Marquardt failed " << "reporting the following: " << gsl_strerror(status) << "\n" @@ -502,7 +496,7 @@ namespace CurveFitting while (status == GSL_CONTINUE && iter < maxInterations) { iter++; - status = gsl_multimin_fminimizer_iterate(simplexMinimizer); + status = minimizer->iterate(); if (status) // break if error { @@ -510,8 +504,8 @@ namespace CurveFitting if (iter == 1) { g_log.information() << "Simplex step size reduced to 0.1\n"; - gsl_vector_set_all (simplexStepSize, 0.1); - gsl_multimin_fminimizer_set(simplexMinimizer, &gslSimplexContainer, initFuncArg, simplexStepSize); + delete minimizer; + minimizer = new SimplexMinimizer(gslSimplexContainer, initFuncArg, 0.1); //iter = 0; status = GSL_CONTINUE; continue; @@ -519,12 +513,11 @@ namespace CurveFitting break; } - size = gsl_multimin_fminimizer_size(simplexMinimizer); - status = gsl_multimin_test_size(size, 1e-2); + status = minimizer->hasConverged(); prog.report(); } - finalCostFuncVal = simplexMinimizer->fval / dof; + finalCostFuncVal = minimizer->costFunctionVal() / dof; } // Output summary to log file @@ -547,12 +540,6 @@ namespace CurveFitting setProperty("Output Status", reportOfFit); setProperty("Output Chi^2/DoF", finalCostFuncVal); - - - // cleanup memory allocated for solvers - - gsl_vector_free(simplexStepSize); - gsl_multimin_fminimizer_free(simplexMinimizer); // if Output property is specified output additional workspaces @@ -623,6 +610,10 @@ namespace CurveFitting } + // minimizer may have dynamically allocated memory hence make sure this memory is freed up + + delete minimizer; + // clean up dynamically allocated gsl stuff delete [] l_data.X; diff --git a/Code/Mantid/CurveFitting/src/LevenbergMarquardtMinimizer.cpp b/Code/Mantid/CurveFitting/src/LevenbergMarquardtMinimizer.cpp index 191a256ac3938534c59c2edccdab9b58239d6afc..d44585a238660877a0957c2dfca7e5929d71c75d 100644 --- a/Code/Mantid/CurveFitting/src/LevenbergMarquardtMinimizer.cpp +++ b/Code/Mantid/CurveFitting/src/LevenbergMarquardtMinimizer.cpp @@ -2,8 +2,6 @@ // Includes //---------------------------------------------------------------------- #include "MantidCurveFitting/LevenbergMarquardtMinimizer.h" -#include <gsl/gsl_errno.h> -#include <gsl/gsl_fit.h> #include <gsl/gsl_blas.h> #include "MantidKernel/Exception.h" diff --git a/Code/Mantid/CurveFitting/src/SimplexMinimizer.cpp b/Code/Mantid/CurveFitting/src/SimplexMinimizer.cpp new file mode 100644 index 0000000000000000000000000000000000000000..533d10823ff221ef8f27590d3868bec211c6f3f1 --- /dev/null +++ b/Code/Mantid/CurveFitting/src/SimplexMinimizer.cpp @@ -0,0 +1,54 @@ +//---------------------------------------------------------------------- +// Includes +//---------------------------------------------------------------------- +#include "MantidCurveFitting/SimplexMinimizer.h" +#include "MantidKernel/Exception.h" + +namespace Mantid +{ +namespace CurveFitting +{ + +SimplexMinimizer::SimplexMinimizer( gsl_multimin_function& gslContainer, + gsl_vector* startGuess, const double& size) : m_name("Simplex") +{ + const gsl_multimin_fminimizer_type *simplexType = gsl_multimin_fminimizer_nmsimplex; + + // step size for simplex + m_simplexStepSize = gsl_vector_alloc(gslContainer.n); + gsl_vector_set_all (m_simplexStepSize, size); // is this always a sensible starting step size? + + // setup minimizer + m_gslSolver = gsl_multimin_fminimizer_alloc(simplexType, gslContainer.n); + gsl_multimin_fminimizer_set(m_gslSolver, &gslContainer, startGuess, m_simplexStepSize); +} + +SimplexMinimizer::~SimplexMinimizer() +{ + gsl_vector_free(m_simplexStepSize); + gsl_multimin_fminimizer_free(m_gslSolver); +} + +std::string SimplexMinimizer::name()const +{ + return m_name; +} + +int SimplexMinimizer::iterate() +{ + return gsl_multimin_fminimizer_iterate(m_gslSolver); +} + +int SimplexMinimizer::hasConverged() +{ + double size = gsl_multimin_fminimizer_size(m_gslSolver); + return gsl_multimin_test_size(size, 1e-2); +} + +double SimplexMinimizer::costFunctionVal() +{ + return m_gslSolver->fval; +} + +} // namespace CurveFitting +} // namespace Mantid