Commit bd3ab796 authored by Roman Tolchenov's avatar Roman Tolchenov
Browse files

Updated Fit to use GenericFit internally. re #2132

parent 4dc9b59b
......@@ -49,6 +49,7 @@ public:
/* Overriden methods */
void setWorkspace(boost::shared_ptr<Workspace> ws,const std::string& slicing);
/// Set the workspace
void setMatrixWorkspace(boost::shared_ptr<const API::MatrixWorkspace> workspace,int spec,int xMin,int xMax);
/// Returns the function's name
......@@ -153,6 +154,8 @@ public:
/// Check the function.
void checkFunction();
void setUpNewStuff(boost::shared_array<double> xs = boost::shared_array<double>(),boost::shared_array<double> weights = boost::shared_array<double>());
protected:
/// Function initialization. Declare function parameters in this method.
void init();
......
......@@ -224,8 +224,10 @@ public:
/// The string operator
virtual operator std::string()const{return asString();}
/// Set the workspace
/// @param wsIDString A string identifying the data to be fitted, e.g. workspace name and spectrum index separated by a comma
virtual void setWorkspace(const std::string& wsIDString) = 0;
/// @param ws Shared pointer to a workspace
/// @param slicing A string identifying the data in the worspace to be fitted, e.g. spectrum index, starting and ending x values, etc
/// Concrete form is defined by the derived functions.
virtual void setWorkspace(boost::shared_ptr<Workspace> ws,const std::string& slicing) = 0;
/// Get the workspace
virtual boost::shared_ptr<const API::Workspace> getWorkspace()const = 0;
/// Iinialize the function
......
......@@ -10,12 +10,19 @@
#include "MantidAPI/IFitFunction.h"
#include "MantidAPI/FunctionFactory.h"
#include "boost/shared_ptr.hpp"
#include "boost/shared_array.hpp"
#include "boost/variant.hpp"
#include <string>
#include <vector>
namespace Mantid
{
namespace CurveFitting
{
class Fit;
}
namespace API
{
#ifndef IGNORE_IFUNCTION_ARGUMENT
......@@ -119,10 +126,7 @@ public:
/* Overidden methods */
/// Set the workspace
/// @param wsIDString A string identifying the data to be fitted, e.g. workspace name and spectrum index separated by a comma
virtual void setWorkspace(const std::string& wsIDString);
/// Get the workspace
virtual void setWorkspace(boost::shared_ptr<Workspace> ws,const std::string& slicing);
virtual boost::shared_ptr<const API::Workspace> getWorkspace()const;
/// Returns the size of the fitted data (number of double values returned by the function)
......@@ -154,7 +158,8 @@ public:
virtual void calJacobianForCovariance(Jacobian* out, const double* xValues, const int& nData);
/// To be used temporarily with the old Fit algorithm.
void setUpNewStuff();
virtual void setUpNewStuff(boost::shared_array<double> xs = boost::shared_array<double>(),
boost::shared_array<double> weights = boost::shared_array<double>());
protected:
......@@ -163,6 +168,8 @@ protected:
boost::shared_ptr<const MatrixWorkspace> ws,
int wsIndex);
boost::shared_ptr<API::MatrixWorkspace> createCalculatedWorkspace(boost::shared_ptr<const API::MatrixWorkspace> inWS,int wi)const;
/// Shared pointer to the workspace
boost::shared_ptr<const API::MatrixWorkspace> m_workspace;
/// Spectrum index
......@@ -176,13 +183,16 @@ protected:
/// Pointer to the fitted data
const double* m_data;
/// Pointer to the fitting weights
std::vector<double> m_weights;
boost::shared_array<double> m_weights;
/// Pointer to the x (function argument as in f(x)) data
const double* m_xValues;
boost::shared_array<double> m_xValues;
/// Static reference to the logger class
static Kernel::Logger& g_log;
/// Making a friend
friend class CurveFitting::Fit;
};
} // namespace API
......
......@@ -711,7 +711,9 @@ void CompositeFunction::setMatrixWorkspace(boost::shared_ptr<const API::MatrixWo
{
IFunction::setMatrixWorkspace(workspace,spec,xMin,xMax);
for(int i=0;i<nFunctions();i++)
{
getFunction(i)->setMatrixWorkspace(workspace,spec,xMin,xMax);
}
}
/**
......@@ -899,5 +901,23 @@ IFitFunction* CompositeFunction::getContainingFunction(const IFitFunction* fun)
return NULL;
}
void CompositeFunction::setWorkspace(boost::shared_ptr<Workspace> ws,const std::string& slicing)
{
IFunction::setWorkspace(ws,slicing);
for(int iFun=0;iFun<nFunctions();iFun++)
{
getFunction(iFun)->setUpNewStuff(m_xValues,m_weights);
}
}
void CompositeFunction::setUpNewStuff(boost::shared_array<double> xs,boost::shared_array<double> weights)
{
IFunction::setUpNewStuff(xs,weights);
for(int iFun=0;iFun<nFunctions();iFun++)
{
getFunction(iFun)->setUpNewStuff(xs,weights);
}
}
} // namespace API
} // namespace Mantid
......@@ -205,13 +205,70 @@ void Expression::tokenize()
bool inString = false;
int skip = 0;
bool canBeBinary = false;
bool isNumber = false; // if parser is inside a number (important case is 123.45e+67)
bool canDotBeAdded = false;
bool canEBeAdded = false;
bool canPlusBeAdded = false;
Tokens tokens;
for(size_t i=0;i<m_expr.size();i++)
{
char c = m_expr[i];
if (!inString && skip == 0)
{
if (lvl == 0 && is_op_symbol(c))// insert new token
if (isNumber)
{
if (c == '.')
{
if (canDotBeAdded)
{
canDotBeAdded = false;
}
else
{
isNumber = false;
}
}
else if (c == 'e' || c == 'E')
{
if (canEBeAdded)
{
canEBeAdded = false;
canDotBeAdded = false;
canPlusBeAdded = true;
}
else
{
isNumber = false;
}
}
else if (c == '+' || c == '-')
{
if (canPlusBeAdded)
{
canPlusBeAdded = false;
canEBeAdded = false;
canDotBeAdded = false;
}
else
{
isNumber = false;
}
}
else if (!isdigit(c))
{
isNumber = false;
}
}
else if (isdigit(c))
{
isNumber = true;
canDotBeAdded = true;
canEBeAdded = true;
canPlusBeAdded = false;
}
if (lvl == 0 && !isNumber && is_op_symbol(c))// insert new token
{
if (i == last)
{
......
......@@ -31,18 +31,16 @@ namespace API
}
/** Base class implementation of derivative IFitFunction throws error. This is to check if such a function is provided
by derivative class. In the derived classes this method must return the derivatives of the resuduals function
(defined in void Fit1D::function(const double*, double*, const double*, const double*, const double*, const int&))
by derivative class. In the derived classes this method must return the derivatives of the function
with respect to the fit parameters. If this method is not reimplemented the derivative free simplex minimization
algorithm is used.
* @param out Derivatives
* @param xValues X values for data points
* @param nData Number of data points
algorithm is used or the derivatives are computed numerically.
* @param out Pointer to a Jacobian matrix. If it is NULL the method is called in order to check whether it's implemented or not.
* If the derivatives are implemented the method must simply return, otherwise it must throw Kernel::Exception::NotImplementedError.
*/
void IFitFunction::functionDeriv(Jacobian* out)
{
(void) out; //Avoid compiler warning
throw Kernel::Exception::NotImplementedError("No derivative IFitFunction provided");
throw ("No derivative IFitFunction provided");
}
......
......@@ -11,6 +11,8 @@
#include "MantidAPI/ConstraintFactory.h"
#include "MantidAPI/AnalysisDataService.h"
#include "MantidAPI/Expression.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidAPI/TextAxis.h"
#include "MantidGeometry/Instrument/ParameterMap.h"
#include "MantidGeometry/Instrument/Component.h"
#include "MantidGeometry/Instrument/DetectorGroup.h"
......@@ -32,64 +34,70 @@ namespace API
Kernel::Logger& IFunction::g_log = Kernel::Logger::get("IFunction");
/// Set the workspace
/// @param wsIDString A string identifying the data to be fitted. Format for IFunction:
/// "WorkspaceName,WorkspaceIndex=int,StartX=double,EndX=double"
void IFunction::setWorkspace(const std::string& wsIDString)
/** Set the workspace
* @param ws A shared pointer to a workspace. Must be a MatrixWorkspace.
* @param slicing A string identifying the data to be fitted. Format for IFunction:
* "WorkspaceIndex=int,StartX=double,EndX=double". StartX and EndX are optional.
*/
void IFunction::setWorkspace(boost::shared_ptr<Workspace> ws,const std::string& slicing)
{
try
{
Expression expr;
expr.parse(wsIDString);
if (expr.name() != "," || expr.size() < 2 ) // TODO: Expression needs pattern matching
{
throw std::invalid_argument("Syntax error in argument");
}
std::string wsName = expr[0].name();
MatrixWorkspace_sptr ws = boost::dynamic_pointer_cast<MatrixWorkspace>(AnalysisDataService::Instance().retrieve(wsName));
if (!ws)
MatrixWorkspace_sptr mws = boost::dynamic_pointer_cast<MatrixWorkspace>(ws);
if (!mws)
{
throw std::invalid_argument("Workspace has wrong type (not a MatrixWorkspace)");
throw std::invalid_argument("Workspace has a wrong type (not a MatrixWorkspace)");
}
int index = -1;
int xMin = -1;
int xMax = -1;
double startX,endX;
for(int i = 1; i < expr.size(); ++i)
int n = int(mws->blocksize()); // length of each Y vector
Expression expr;
expr.parse(slicing);
if (expr.name() == ",")
{
const Expression& e = expr[i];
if (e.size() == 2 && e.name() == "=")
for(int i = 0; i < expr.size(); ++i)
{
if (e[0].name() == "WorkspaceIndex")
{
index = boost::lexical_cast<int>(e[1].name());
}
else if (e[0].name() == "StartX")
const Expression& e = expr[i];
if (e.size() == 2 && e.name() == "=")
{
startX = boost::lexical_cast<double>(e[1].name());
xMin = 0;
}
else if (e[0].name() == "EndX")
{
endX = boost::lexical_cast<double>(e[1].name());
xMax = 0;
if (e[0].name() == "WorkspaceIndex")
{
index = boost::lexical_cast<int>(e[1].name());
}
else if (e[0].name() == "StartX")
{
startX = boost::lexical_cast<double>(e[1].name());
xMin = 0;
}
else if (e[0].name() == "EndX")
{
endX = boost::lexical_cast<double>(e[1].name());
xMax = 0;
}
}
}
}
else if (expr.size() == 2 && expr.name() == "=" && expr[0].name() == "WorkspaceIndex")
{
index = boost::lexical_cast<int>(expr[1].name());
}
if (index < 0)
{
g_log.warning("WorkspaceIndex not set, defaulting to 0");
}
else if (index >= ws->getNumberHistograms())
else if (index >= mws->getNumberHistograms())
{
throw std::range_error("WorkspaceIndex outside range");
}
const MantidVec& x = ws->readX(index);
const MantidVec& y = ws->readY(index);
const MantidVec& e = ws->readE(index);
int n = int(y.size());
const MantidVec& x = mws->readX(index);
const MantidVec& y = mws->readY(index);
const MantidVec& e = mws->readE(index);
if (xMin < 0)
{
xMin = 0;
......@@ -132,20 +140,29 @@ namespace API
}
m_dataSize = xMax - xMin;
m_data = &y[xMin];
m_xValues = &x[xMin];
m_weights.resize(m_dataSize);
m_xValues.reset(new double[m_dataSize]);
m_weights.reset(new double[m_dataSize]);
bool isHist = x.size() > y.size();
for (int i = 0; i < m_dataSize; ++i)
{
if (isHist)
{
m_xValues[i] = 0.5*(x[xMin + i] + x[xMin + i + 1]);
}
else
{
m_xValues[i] = x[xMin + i];
}
if (e[xMin + i] <= 0.0)
m_weights[i] = 1.0;
else
m_weights[i] = 1./e[xMin + i];
}
if (ws->hasMaskedBins(index))
if (mws->hasMaskedBins(index))
{
const MatrixWorkspace::MaskList& mlist = ws->maskedBins(index);
const MatrixWorkspace::MaskList& mlist = mws->maskedBins(index);
MatrixWorkspace::MaskList::const_iterator it = mlist.begin();
for(;it!=mlist.end();it++)
{
......@@ -153,7 +170,7 @@ namespace API
}
}
setMatrixWorkspace(ws,index,xMin,xMax);
setMatrixWorkspace(mws,index,xMin,xMax);
}
catch(std::exception& e)
......@@ -190,8 +207,8 @@ namespace API
/// @param out The buffer for writing the calculated values. Must be big enough to accept dataSize() values
void IFunction::function(double* out)const
{
if (m_weights.empty()) return;
function(out,m_xValues,m_dataSize);
if (m_dataSize == 0) return;
function(out,m_xValues.get(),m_dataSize);
// Add penalty factor to function if any constraint is violated
double penalty = 0.;
......@@ -221,8 +238,13 @@ namespace API
/// Derivatives of function with respect to active parameters
void IFunction::functionDeriv(Jacobian* out)
{
if (m_weights.empty()) return;
functionDeriv(out,m_xValues,m_dataSize);
if (out == NULL)
{
functionDeriv(out,m_xValues.get(),0);
return;
}
if (m_dataSize == 0) return;
functionDeriv(out,m_xValues.get(),m_dataSize);
if (m_dataSize <= 0) return;
......@@ -480,33 +502,100 @@ void IFunction::calJacobianForCovariance(Jacobian* out, const double* xValues, c
}
/// Called after setMatrixWorkspace if setWorkspace hadn't been called before
void IFunction::setUpNewStuff()
void IFunction::setUpNewStuff(boost::shared_array<double> xs,boost::shared_array<double> weights)
{
m_dataSize = m_xMaxIndex - m_xMinIndex;
m_data = &m_workspace->readY(m_workspaceIndex)[m_xMinIndex];
m_xValues = &m_workspace->readX(m_workspaceIndex)[m_xMinIndex];
m_weights.resize(m_dataSize);
const MantidVec& e = m_workspace->readE(m_workspaceIndex);
for (int i = 0; i < m_dataSize; ++i)
//m_xValues = &m_workspace->readX(m_workspaceIndex)[m_xMinIndex];
if (weights && xs)
{
if (e[m_xMinIndex + i] <= 0.0)
m_weights[i] = 1.0;
else
m_weights[i] = 1./e[m_xMinIndex + i];
m_xValues = xs;
m_weights = weights;
}
if (m_workspace->hasMaskedBins(m_workspaceIndex))
else
{
const MatrixWorkspace::MaskList& mlist = m_workspace->maskedBins(m_workspaceIndex);
MatrixWorkspace::MaskList::const_iterator it = mlist.begin();
for(;it!=mlist.end();it++)
m_xValues.reset(new double[m_dataSize]);
m_weights.reset(new double[m_dataSize]);
const MantidVec& x = m_workspace->readX(m_workspaceIndex);
const MantidVec& e = m_workspace->readE(m_workspaceIndex);
bool isHist = m_workspace->isHistogramData();
for (int i = 0; i < m_dataSize; ++i)
{
if (isHist)
{
m_xValues[i] = 0.5*(x[m_xMinIndex + i] + x[m_xMinIndex + i + 1]);
}
else
{
m_xValues[i] = x[m_xMinIndex + i];
}
if (e[m_xMinIndex + i] <= 0.0)
m_weights[i] = 1.0;
else
m_weights[i] = 1./e[m_xMinIndex + i];
}
if (m_workspace->hasMaskedBins(m_workspaceIndex))
{
m_weights[it->first - m_xMinIndex] = 0.;
const MatrixWorkspace::MaskList& mlist = m_workspace->maskedBins(m_workspaceIndex);
MatrixWorkspace::MaskList::const_iterator it = mlist.begin();
for(;it!=mlist.end();it++)
{
m_weights[it->first - m_xMinIndex] = 0.;
}
}
}
}
boost::shared_ptr<API::MatrixWorkspace> IFunction::createCalculatedWorkspace(boost::shared_ptr<const API::MatrixWorkspace> inWS,int wi)const
{
const MantidVec& inputX = inWS->readX(wi);
const MantidVec& inputY = inWS->readY(wi);
int nData = dataSize();
int histN = inWS->isHistogramData() ? 1 : 0;
API::MatrixWorkspace_sptr ws =
Mantid::API::WorkspaceFactory::Instance().create(
"Workspace2D",
3,
nData + histN,
nData);
ws->setTitle("");
ws->setYUnitLabel(inWS->YUnitLabel());
ws->setYUnit(inWS->YUnit());
ws->getAxis(0)->unit() = inWS->getAxis(0)->unit();
API::TextAxis* tAxis = new API::TextAxis(3);
tAxis->setLabel(0,"Data");
tAxis->setLabel(1,"Calc");
tAxis->setLabel(2,"Diff");
ws->replaceAxis(1,tAxis);
for(int i=0;i<3;i++)
{
ws->dataX(i).assign(inputX.begin()+m_xMinIndex,inputX.begin()+m_xMaxIndex+histN);
}
ws->dataY(0).assign(inputY.begin()+m_xMinIndex,inputY.begin()+m_xMaxIndex);
MantidVec& Ycal = ws->dataY(1);
MantidVec& E = ws->dataY(2);
double* lOut = new double[nData]; // to capture output from call to function()
function( lOut );
for(int i=0; i<nData; i++)
{
Ycal[i] = lOut[i];
E[i] = m_data[i] - Ycal[i];
}
delete [] lOut;
return ws;
}
} // namespace API
} // namespace Mantid
......@@ -84,7 +84,7 @@ set ( TEST_FILES #test/BackToBackExponentialTest.h # TODO tests are commented ou
test/FunctionFactoryTest.h
test/FunctionTest.h
test/Gaussian1DTest.h
#test/GaussianTest.h # TODO tests are commented out
test/GaussianTest.h # TODO tests are commented out
#test/IkedaCarpenterPVTest.h # TODO tests are commented out
test/LinearBackgroundTest.h
test/LinearTest.h
......@@ -96,7 +96,8 @@ set ( TEST_FILES #test/BackToBackExponentialTest.h # TODO tests are commented ou
test/SpecialFunctionSupportTest.h
test/SplineBackgroundTest.h
test/UserFunction1DTest.h
test/UserFunctionTest.h )
test/UserFunctionTest.h
)
# Add the target for this directory
add_library ( CurveFitting ${SRC_FILES} ${INC_FILES})
......
......@@ -54,6 +54,7 @@ public:
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, API::IFitFunction* function, const std::string& costFunction);
void initialize(API::IFitFunction* function, const std::string& costFunction);
private:
/// name of this minimizer
......
......@@ -55,6 +55,8 @@ public:
void initialize(double* X, const double* Y, double *sqrtWeight, const int& nData, const int& nParam,
gsl_vector* startGuess, API::IFitFunction* function, const std::string& costFunction);
void initialize(API::IFitFunction* function, const std::string& costFunction);
private:
/// name of this minimizer
const std::string m_name;
......
......@@ -73,8 +73,6 @@ namespace Mantid
/// Derivatives of function with respect to parameters you are trying to fit
//virtual void functionDeriv(const double* in, API::Jacobian* out, const double* xValues, const int& nData);
/// Set a function for fitting
void setFunction(API::IFunction* fun);
/// Get the function for fitting
API::IFunction* getFunction()const{return m_function;}
......@@ -82,6 +80,10 @@ namespace Mantid
// Overridden Algorithm methods
void init();
void exec();
void exec1();
/// Set a function for fitting
void setFunction(API::IFunction* fun);
/// Option for providing intelligent range starting value based e.g. on the user input parameter values
virtual void modifyStartOfRange(double& startX)
......@@ -120,6 +122,9 @@ namespace Mantid
/// Pointer to the fitting function
API::IFunction* m_function;
/// Function initialization string
std::string m_function_input;