Commit 3353f0d4 authored by Roman Tolchenov's avatar Roman Tolchenov
Browse files

Added MultiBG function for doing a multispectral fit. Updated...

Added MultiBG function for doing a multispectral fit. Updated FitPropertyBrowser to be able to use it. re #966
parent 927827bb
......@@ -14,6 +14,7 @@ set ( SRC_FILES
src/Column.cpp
src/ColumnFactory.cpp
src/CompositeFunction.cpp
src/CompositeFunctionMD.cpp
src/CompositeFunctionMW.cpp
src/ConstraintFactory.cpp
src/CostFunctionFactory.cpp
......
......@@ -107,6 +107,13 @@ namespace Mantid
* @param newName :: The new name
*/
void rename(const std::string& newName);
/**
* Make sure the expression is a list of expression separated by sep, eg "term1,term2,..."
* If it's not a list turn it into one, eg "expr,"
* @param sep :: Separator
*/
void toList(const std::string& sep = ",");
private:
/// copy contructor
......
......@@ -71,6 +71,10 @@ void CompositeFunction::init()
std::string CompositeFunction::asString()const
{
std::ostringstream ostr;
if (name() != "CompositeFunctionMW")
{
ostr << "composite=" <<name() << ";";
}
for(size_t i=0;i<nFunctions();i++)
{
IFitFunction* fun = getFunction(i);
......
......@@ -145,7 +145,8 @@ void CompositeFunctionMW::setMatrixWorkspace(boost::shared_ptr<const API::Matrix
for(size_t i=0;i<nFunctions();i++)
{
IFunctionMW* fun = dynamic_cast<IFunctionMW*>(getFunction(i));
fun->setMatrixWorkspace(workspace,spec,xMin,xMax);
if (fun)
fun->setMatrixWorkspace(workspace,spec,xMin,xMax);
}
}
......@@ -157,7 +158,8 @@ void CompositeFunctionMW::setWorkspace(boost::shared_ptr<const Workspace> ws,con
{
IFunctionMW* fun = dynamic_cast<IFunctionMW*>(getFunction(iFun));
//fun->setWorkspace(ws, slicing, copyData); // TODO: This was added by JZ May 13, 2011, to fix tests. Does this make sense to someone who knows?
fun->setUpNewStuff(m_xValues,m_weights);
if (fun)
fun->setUpNewStuff(m_xValues,m_weights);
}
}
......
......@@ -649,6 +649,15 @@ void Expression::renameAll(const std::string& oldName,const std::string& newName
}
}
void Expression::toList(const std::string& sep)
{
if (name() == sep) return;
Expression term(*this);
m_terms.resize(1);
m_terms[0] = term;
setFunct(sep);
}
}//API
}//Mantid
......@@ -758,6 +758,5 @@ void IFunctionMW::calNumericalDeriv(Jacobian* out, const double* xValues, const
}
}
} // namespace API
} // namespace Mantid
......@@ -379,6 +379,14 @@ public:
TS_ASSERT(vars.find("x")!=vars.end());
}
void testToList()
{
Expression e;
e.parse("x");
e.toList();
TS_ASSERT_EQUALS(e.name(),",");
}
};
#endif /*EXPRESSIONTEST_H_*/
......@@ -25,6 +25,7 @@ set ( SRC_FILES
src/LinearBackground.cpp
src/Lorentzian.cpp
src/Lorentzian1D.cpp
src/MultiBG.cpp
src/PRConjugateGradientMinimizer.cpp
src/PlotPeakByLogValue.cpp
src/ProductFunctionMW.cpp
......@@ -68,6 +69,7 @@ set ( INC_FILES
inc/MantidCurveFitting/LinearBackground.h
inc/MantidCurveFitting/Lorentzian.h
inc/MantidCurveFitting/Lorentzian1D.h
inc/MantidCurveFitting/MultiBG.h
inc/MantidCurveFitting/PRConjugateGradientMinimizer.h
inc/MantidCurveFitting/PlotPeakByLogValue.h
inc/MantidCurveFitting/ProductFunctionMW.h
......@@ -99,6 +101,7 @@ set ( TEST_FILES
test/LinearTest.h
test/Lorentzian1DTest.h
test/LorentzianTest.h
test/MultiBGTest.h
test/PlotPeakByLogValueTest.h
test/ProductFunctionMWTest.h
test/QuadraticTest.h
......
#ifndef MANTID_CURVEFITTING_MULTIBG_H_
#define MANTID_CURVEFITTING_MULTIBG_H_
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAPI/CompositeFunction.h"
#include "MantidAPI/IFunctionMW.h"
#include "MantidAPI/MatrixWorkspace.h"
#ifdef _WIN32
#pragma warning( disable: 4250 )
#endif
namespace Mantid
{
namespace CurveFitting
{
/** A composite function.
@author Roman Tolchenov, Tessella Support Services plc
@date 19/08/2011
Copyright &copy; 2009 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National 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 MultiBG : public API::CompositeFunction
{
public:
/// Default constructor
MultiBG():API::CompositeFunction(){}
///Destructor
virtual ~MultiBG();
/* Overriden methods */
/// Set the workspace
void setWorkspace(boost::shared_ptr<const API::Workspace> ws,const std::string& slicing, bool copyData = true);
virtual boost::shared_ptr<const API::Workspace> getWorkspace()const{return m_spectra[0].first;}
/// Returns the function's name
std::string name()const{return "MultiBG";}
/// Returns the size of the fitted data (number of double values returned by the function)
virtual int dataSize()const{return static_cast<int>(m_data.size());}
/// Returns a reference to the fitted data. These data are taken from the workspace set by setWorkspace() method.
virtual const double* getData()const{return &m_data[0];}
virtual const double* getWeights()const{return &m_weights[0];}
/// Function you want to fit to.
void function(double* out)const;
/// Derivatives of function with respect to active parameters
void functionDeriv(API::Jacobian* out);
protected:
/// to collect different workspaces found in child functions
std::vector< std::pair< boost::shared_ptr<const API::MatrixWorkspace>, size_t> > m_spectra;
/// corresponding MatrixWorkspace indexes
//std::vector< size_t > m_wsIndex;
/// to store function indices to workspaces: m_funIndex[i] gives vector of indexes of m_spectra for function i
std::vector< std::vector<size_t> > m_funIndex;
/// the data vector which is a composition of all fitted spectra
std::vector<double> m_data;
/// the vector of fitting weights, one for each value in m_data
std::vector<double> m_weights;
/// offsets of particular workspaces in the m_data and m_weights arrays
std::vector<size_t> m_offset;
};
} // namespace CurveFitting
} // namespace Mantid
#endif /*MANTID_CURVEFITTING_MULTIBG_H_*/
......@@ -71,7 +71,7 @@ void Gaussian::setActiveParameter(int i,double value)
int j = indexOfActive(i);
if (parameterName(j) == "Sigma")
setParameter(j,sqrt(1./value),false);
setParameter(j,sqrt(fabs(1./value)),false);
else
setParameter(j,value,false);
}
......
......@@ -169,7 +169,9 @@ namespace CurveFitting
while (status == GSL_CONTINUE && iter < maxInterations)
{
iter++;
status = minimizer->iterate();
if (status != GSL_SUCCESS && minimizer->hasConverged() != GSL_SUCCESS)
{
// From experience it is found that gsl_multifit_fdfsolver_iterate occasionally get
......
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidCurveFitting/MultiBG.h"
#include "MantidAPI/IFunctionMW.h"
#include "MantidAPI/Expression.h"
#include "MantidAPI/AnalysisDataService.h"
#include <boost/lambda/lambda.hpp>
#include <boost/lexical_cast.hpp>
#include <sstream>
#include <iostream>
#include <algorithm>
#include <iterator>
#include <float.h>
using namespace boost::lambda;
namespace Mantid
{
namespace CurveFitting
{
DECLARE_FUNCTION(MultiBG)
///Destructor
MultiBG::~MultiBG()
{
}
/**
* Function you want to fit to.
*/
void MultiBG::function(double* out)const
{
std::vector<double> tmpOut(dataSize());
std::fill_n(out,dataSize(),0);
for(size_t i = 0; i < nFunctions(); i++)
{
IFitFunction* fun = getFunction(i);
size_t nWS = m_funIndex[i].size();
for(size_t k = 0;k < nWS; ++k)
{
size_t j = m_funIndex[i][k];
fun->setWorkspace(m_spectra[k].first,"WorkspaceIndex="+boost::lexical_cast<std::string>(m_spectra[j].second),false);
//std::cerr << i << ' ' << k << " Function " << fun->name() << " ws " << fun->getWorkspace()->getName() << " wi "
// << dynamic_cast<Mantid::API::IFunctionMW*>(fun)->getWorkspaceIndex() << std::endl;
double *out1 = out + m_offset[j];
double *tmp1 = &tmpOut[0] + m_offset[j];
size_t nData = 0;
if (j < m_offset.size() - 1 )
nData = m_offset[j + 1] - m_offset[j];
else
nData = dataSize() - m_offset[j];
if (i == 0)
{
fun->function(out1);
}
else
{
fun->function(tmp1);
std::transform(out1, out1 + nData, tmp1, out1, std::plus<double>());
}
}
}
//std::cerr << "Function:\n";
//for(size_t i = 0; i<nParams();++i)
//{
// std::cerr << getParameter(i) << ' ' ;
//}
//std::cerr << std::endl;
//std::for_each(out,out+m_dataSize,std::cerr << _1 << '\n');
//std::cerr << std::endl;
}
/// Derivatives of function with respect to active parameters
void MultiBG::functionDeriv(API::Jacobian* out)
{
// it is possible that out is NULL
if (!out) return;
// claculate numerically
double stepPercentage = DBL_EPSILON*1000; // step percentage
double step; // real step
const int nParam = nParams();
const int nData = dataSize();
//for(size_t i=0; i < nParams(); ++i)
// std::cerr << i << ' ' << getParameter(i) << std::endl;
//std::cerr << std::endl;
std::vector<double> tmpFunctionOutputMinusStep(nData);
std::vector<double> tmpFunctionOutputPlusStep(nData);
function(&tmpFunctionOutputMinusStep[0]);
for (int iP = 0; iP < nParam; iP++)
{
if ( isActive(iP) )
{
const double& val = getParameter(iP);
if (fabs(val) < stepPercentage)
{
step = stepPercentage;
}
else
{
step = val*stepPercentage;
}
double paramPstep = val + step;
setParameter(iP, paramPstep);
function(&tmpFunctionOutputPlusStep[0]);
step = paramPstep - val;
setParameter(iP, val);
for (int i = 0; i < nData; i++)
{
double value = (tmpFunctionOutputPlusStep[i]-tmpFunctionOutputMinusStep[i])/step;
out->set(i,iP,value);
}
}
}
}
/**
* Sets workspaces to member functions. Constructs the data set for fitting.
* @param ws :: Pointer to a workspace, not used. Workspaces are taken either from member functions or slicing.
* @param slicing :: A map between member functions and workspaces or empty string. Format:
* "f0,Workspace0,i0;f1,Workspace1,i1;f2,Workspace2,i2;..."
*/
void MultiBG::setWorkspace(boost::shared_ptr<const API::Workspace> ws,const std::string& slicing, bool)
{
boost::shared_ptr<const API::MatrixWorkspace> mws = boost::dynamic_pointer_cast<const API::MatrixWorkspace>(ws);
if (ws && !mws)
{
throw std::invalid_argument("Workspace in MultiBG has a wrong type (not a MatrixWorkspace)");
}
m_funIndex.resize(nFunctions());
if (! slicing.empty())
{
Mantid::API::Expression expr;
expr.parse(slicing);
// expr can be treated as a list even if it has only 1 term
expr.toList(";");
for(size_t i=0;i<expr.size();++i)
{
const Mantid::API::Expression& e = expr[i];
if (e.name() != "," || e.size() != 3)
{
// slicing has a wrong format - ignore it
break;
}
try
{
std::string wsName = e[1].name();
Mantid::API::MatrixWorkspace_sptr ws = boost::dynamic_pointer_cast<Mantid::API::MatrixWorkspace>(
Mantid::API::AnalysisDataService::Instance().retrieve(wsName));
size_t iFun = boost::lexical_cast<size_t>(e[0].name().substr(1));
size_t wi = boost::lexical_cast<size_t>(e[2].name());
if (iFun >= nFunctions())
{
throw std::invalid_argument("MultiBG::setWorkspace: function "+e[0].name()+" not found");
}
std::pair< boost::shared_ptr<const API::MatrixWorkspace>, size_t> spectrum = std::make_pair(ws,wi);
m_funIndex[iFun].push_back(m_spectra.size());
m_spectra.push_back(spectrum);
getFunction(iFun)->setWorkspace(ws,"WorkspaceIndex="+e[2].name());
}
catch(...)
{
break;
}
}
}
// examine the member functions and fill in the m_funIndex array
for(size_t iFun=0;iFun<nFunctions();iFun++)
{
API::IFunctionMW* fun = dynamic_cast<API::IFunctionMW*>(getFunction(iFun));
if (! fun )
{
throw std::runtime_error("MultiBG works with IFunctionMW only");
}
if (fun->getWorkspace())
{
boost::shared_ptr<const API::MatrixWorkspace> iws = fun->getMatrixWorkspace();
std::pair< boost::shared_ptr<const API::MatrixWorkspace>, size_t> spectrum = std::make_pair(iws,fun->getWorkspaceIndex());
std::vector< std::pair< boost::shared_ptr<const API::MatrixWorkspace>, size_t> >::iterator it = std::find(m_spectra.begin(),m_spectra.end(),spectrum);
size_t i;
if (it == m_spectra.end())
{
i = m_spectra.size();
m_spectra.push_back(spectrum);
}
else
{
i = size_t(std::distance(it,m_spectra.begin()));
}
m_funIndex[iFun].push_back(i);
//fun->setWorkspace(boost::static_pointer_cast<const API::Workspace>(iws),slicing,false);
}
}
// make functions without set workspace fit to all workspaces
for(size_t iFun=0;iFun<nFunctions();iFun++)
{
std::vector<size_t>& index = m_funIndex[iFun];
if (index.empty())
{
index.resize(m_spectra.size());
int i = 0;
std::for_each(index.begin(),index.end(),_1 = var(i)++);
getFunction(iFun)->setWorkspace(m_spectra[0].first,"WorkspaceIndex="+boost::lexical_cast<std::string>(m_spectra[0].second));
}
}
// set dimensions and calculate ws's contribution to m_dataSize
//IFunctionMW::setWorkspace(ws,slicing,false);
// add other workspaces
m_offset.resize(m_spectra.size(),0);
size_t nData = 0;
for(size_t i = 0; i < m_spectra.size(); ++i)
{
mws = m_spectra[i].first;
size_t n = mws->blocksize();
m_offset[i] = nData;
nData += static_cast<int>(n);
}
m_data.resize(nData);
m_weights.resize(nData);
//... fill in the data and the weights ...
for(size_t i = 0; i < m_spectra.size(); ++i)
{
mws = m_spectra[i].first;
size_t wi = m_spectra[i].second;
const Mantid::MantidVec& Y = mws->readY(wi);
const Mantid::MantidVec& E = mws->readE(wi);
size_t j0 = m_offset[i];
for(size_t j = 0; j < Y.size(); ++j)
{
m_data[j0 + j] = Y[j];
double err = E[j];
m_weights[j0 + j] = err != 0.0 ? 1./err : 1.0;
}
}
//std::cerr << "Workspace:\n";
//std::for_each(&m_data[0],&m_data[0]+m_dataSize,std::cerr << _1 << '\n');
}
} // namespace API
} // namespace Mantid
......@@ -577,7 +577,7 @@ public:
{
// Should be all of them
// TODO: Should this be 35 or 36?
doFunctionNameTest<IFitFunction>(35, "", "");
doFunctionNameTest<IFitFunction>(37, "", "");
}
void test_PeakFunction_Name_Retrieval()
......
#ifndef MULTIBGTEST_H_
#define MULTIBGTEST_H_
#include <cxxtest/TestSuite.h>
#include "MantidCurveFitting/MultiBG.h"
#include "MantidCurveFitting/Fit.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidAPI/AnalysisDataService.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/AlgorithmFactory.h"
#include <fstream>
using namespace Mantid::API;
using namespace Mantid::CurveFitting;
/**
* A simple implementation of Jacobian.
*/
class SimpleJacobian: public Jacobian
{
public:
/// Constructor
SimpleJacobian(size_t nData,size_t nParams):m_nData(nData),m_nParams(nParams),m_data(nData*nParams){}
/// Setter
virtual void set(size_t iY, size_t iP, double value)
{
m_data[iY * m_nParams + iP] = value;
}
/// Getter
virtual double get(size_t iY, size_t iP)
{
return m_data[iY * m_nParams + iP];
}
double *array(){return &m_data[0];}
size_t n1()const{return m_nData;}
size_t n2()const{return m_nParams;}
private:
size_t m_nData; ///< size of the data / first dimension
size_t m_nParams; ///< number of parameters / second dimension
std::vector<double> m_data; ///< data storage
};
class MultiBGTest : public CxxTest::TestSuite
{
public:
MultiBGTest()
{
a.resize(2); b.resize(2);
a[0] = 2.; b[0] = -.10;
a[1] = -1.0; b[1] = .10;
h = 20.0;
s = 1.0;
c = 0.1;
numBins = 31;
x0 = -10.0;
dx = 20.0 / numBins;
Mantid::DataObjects::Workspace2D_sptr WS(new Mantid::DataObjects::Workspace2D);
WS->initialize(2,numBins+1,numBins);
for (int i = 0; i < numBins + 1; ++i)
{
WS->dataX(0)[i] = x0 + dx * i;
}
for(size_t spec = 0; spec < WS->getNumberHistograms(); ++spec)
{
WS->dataX(spec) = WS->readX(0);
for (int i = 0; i < numBins; ++i)
{
const double x = (WS->readX(0)[i] + WS->readX(0)[i+1])/2;
WS->dataY(spec)[i] = a[spec] + b[spec] * x + h * exp(-0.5 * (x - c) * (x - c) /s/s );// + noise();
WS->dataE(spec)[i] = 1.0;//sqrt(fabs(WS->dataY(spec)[i]));
}
}
std::ofstream fil("MultiBGTestWS.txt");
for(int i = 0; i < numBins; ++i)
{
fil << (WS->readX(0)[i] + WS->readX(0)[i+1])/2 << ',' << WS->readY(0)[i] << ',' << WS->readE(0)[i] << ',' << WS->readY(1)[i] << ',' << WS->readE(1)[i] << std::endl;
}
AnalysisDataService::Instance().add("MultiBGTestWS",WS);
Mantid::API::IPeakFunction::setPeakRadius(1000);
}
void testCorrectDataVectorConstrucion()
{
std::string funIni = "composite=MultiBG;name=Gaussian,Height=100.,Sigma=0.00100,PeakCentre=0;"