Skip to content

Commit

Permalink
Re #4158. Returned the old Levenberg-Marquardt minimizer.
Browse files Browse the repository at this point in the history
  • Loading branch information
@mantid-roman
mantid-roman committed Apr 2, 2012
1 parent 4f44d76 commit aa9b247
Show file tree
Hide file tree
Showing 45 changed files with 1,423 additions and 901 deletions.
5 changes: 5 additions & 0 deletions Code/Mantid/Framework/API/inc/MantidAPI/FunctionValues.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,12 +57,17 @@ class MANTID_API_DLL FunctionValues
size_t size() const {return m_calculated.size();}
/// store i-th calculated value. 0 <= i < size()
void setCalculated(size_t i,double value) {m_calculated[i] = value;}
/// set all calculated values to same number
void setCalculated(double value);
/// get i-th calculated value. 0 <= i < size()
double getCalculated(size_t i) const {return m_calculated[i];}
void addToCalculated(size_t i, double value) {m_calculated[i] += value;}
/// Get a pointer to calculated data at index i
double* getPointerToCalculated(size_t i);
/// Add other calculated values
FunctionValues& operator+=(const FunctionValues& values);
/// Multiply by other calculated values
FunctionValues& operator*=(const FunctionValues& values);
/// Set all calculated values to zero
void zeroCalculated();

Expand Down
13 changes: 11 additions & 2 deletions Code/Mantid/Framework/API/inc/MantidAPI/IFunctionMW.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,9 +4,12 @@
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAPI/IFunction1D.h"
#include "MantidAPI/IFunction.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidKernel/Unit.h"

#include <boost/weak_ptr.hpp>

namespace Mantid
{

Expand Down Expand Up @@ -40,12 +43,14 @@ namespace API
File change history is stored at: <https://svn.mantidproject.org/mantid/trunk/Code/Mantid>.
Code Documentation is available at: <http://doxygen.mantidproject.org>
*/
class MANTID_API_DLL IFunctionMW: public virtual IFunction1D
class MANTID_API_DLL IFunctionMW: public virtual IFunction
{
public:

void setMatrixWorkspace(boost::shared_ptr<const API::MatrixWorkspace> workspace,size_t wi,double startX, double endX);

boost::shared_ptr<const API::MatrixWorkspace> getMatrixWorkspace() const;
size_t getWorkspaceIndex() const {return m_workspaceIndex;}
protected:

/// Convert a value from one unit (inUnit) to unit defined in workspace (ws)
Expand All @@ -56,6 +61,10 @@ class MANTID_API_DLL IFunctionMW: public virtual IFunction1D
void convertValue(std::vector<double>& values, Kernel::Unit_sptr& outUnit,
boost::shared_ptr<const MatrixWorkspace> ws,
size_t wsIndex) const;

boost::weak_ptr<const API::MatrixWorkspace> m_workspace;

size_t m_workspaceIndex;

/// Static reference to the logger class
static Kernel::Logger& g_log;
Expand Down
23 changes: 22 additions & 1 deletion Code/Mantid/Framework/API/src/FunctionValues.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,13 @@ namespace API
m_calculated.resize(domain.size());
}


/// set all calculated values to same number
void FunctionValues::setCalculated(double value)
{
std::fill(m_calculated.begin(),m_calculated.end(),value);
}

/**
* Get a pointer to calculated data at index i
* @param i :: Index.
Expand Down Expand Up @@ -66,10 +73,24 @@ namespace API
return *this;
}

/** Multiply this calculated values by others.
* @param values :: Some other values to be added to this calculated values.
*/
FunctionValues& FunctionValues::operator*=(const FunctionValues& values)
{
if (values.size() != size())
{
throw std::runtime_error("Cannot multiply function values: different sizes.");
}
std::transform(m_calculated.begin(),m_calculated.end(),values.m_calculated.begin(),m_calculated.begin(),
std::multiplies<double>());
return *this;
}

/// Set all calculated values to zero
void FunctionValues::zeroCalculated()
{
std::fill(m_calculated.begin(),m_calculated.end(),0.0);
setCalculated(0.0);
}


Expand Down
5 changes: 4 additions & 1 deletion Code/Mantid/Framework/API/src/IFunction.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -472,6 +472,7 @@ void IFunction::calNumericalDeriv(const FunctionDomain& domain, Jacobian& jacobi
m_plusStep.reset(domain);
}

applyTies(); // just in case
function(domain,m_minusStep);

for (size_t iP = 0; iP < nParam; iP++)
Expand All @@ -490,12 +491,14 @@ void IFunction::calNumericalDeriv(const FunctionDomain& domain, Jacobian& jacobi

double paramPstep = val + step;
setActiveParameter(iP, paramPstep);
applyTies();
function(domain,m_plusStep);

step = paramPstep - val;
setActiveParameter(iP, val);

for (size_t i = 0; i < nData; i++) {
for (size_t i = 0; i < nData; i++)
{
jacobian.set(i,iP,
(m_plusStep.getCalculated(i) - m_minusStep.getCalculated(i))/step);
}
Expand Down
217 changes: 11 additions & 206 deletions Code/Mantid/Framework/API/src/IFunctionMW.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,97 +37,10 @@ namespace API
*/
void IFunctionMW::setMatrixWorkspace(boost::shared_ptr<const API::MatrixWorkspace> workspace,size_t wi,double startX, double endX)
{
if (!workspace) return; // unset the workspace

//size_t n = workspace->blocksize(); // length of each Y vector
//size_t xMin = 0;
//size_t xMax = 0;

//if (wi >= workspace->getNumberHistograms())
//{
// throw std::range_error("Workspace index out of range");
//}

//const MantidVec& x = workspace->readX(wi);
//const MantidVec& y = workspace->readY(wi);
//const MantidVec& e = workspace->readE(wi);
//bool isHist = x.size() > y.size();

//if (startX >= endX)
//{
// xMin = 0;
// xMax = n - 1;
//}
//else
//{
// size_t m = isHist? n - 1 : n;
// for(; xMax < m; ++xMax)
// {
// if (x[xMax] > endX)
// {
// if (xMax > 0) xMax--;
// break;
// }
// if (x[xMax] <= startX)
// {
// xMin = xMax;
// }
// }
//}

//if (xMin > xMax)
//{
// std::swap(xMin,xMax);
//}

//if (xMax >= n) xMax = n - 1;

//m_xMinIndex = xMin;
//m_xMaxIndex = xMax;
m_workspace = workspace;
m_workspaceIndex = wi;

//m_dataSize = xMax - xMin + 1;
//m_data = &y[xMin];
//m_xValues.reset(new double[m_dataSize]);
//m_weights.reset(new double[m_dataSize]);

//bool negativeError = false;

//for (size_t 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 if(e[xMin + i] < 0.0)
// {
// negativeError = true;
// m_weights[i] = 1./fabs(e[xMin + i]); //1.0;
// }
// else
// m_weights[i] = 1./e[xMin + i];
//}

//if ( negativeError )
// g_log.warning() << "Negative error values found! These are set to absolute value\n";

//if (workspace->hasMaskedBins(wi))
//{
// const MatrixWorkspace::MaskList& mlist = workspace->maskedBins(wi);
// MatrixWorkspace::MaskList::const_iterator it = mlist.begin();
// for(;it!=mlist.end();++it)
// {
// m_weights[it->first - xMin] = 0.;
// }
//}
if (!workspace) return; // unset the workspace

try
{
Expand Down Expand Up @@ -285,6 +198,14 @@ void IFunctionMW::setMatrixWorkspace(boost::shared_ptr<const API::MatrixWorkspac
}
}

/**
* Get a shared pointer to the saved matrix workspace.
*/
boost::shared_ptr<const API::MatrixWorkspace> IFunctionMW::getMatrixWorkspace() const
{
return m_workspace.lock();
}

/** Convert a value from unit defined in workspace (ws) to outUnit
*
* @param value :: assumed to be in unit of workspace
Expand Down Expand Up @@ -405,121 +326,5 @@ void IFunctionMW::convertValue(std::vector<double>& values, Kernel::Unit_sptr& o
}
}

namespace
{
/**
* 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];
}
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
};
}

/**
* Creates a workspace containing values calculated with this function. It takes a workspace and ws index
* of a spectrum which this function may have been fitted to. The output contains the original spectrum
* (wi = 0), the calculated values (ws = 1), and the difference between them (ws = 2).
* @param inWS :: input workspace
* @param wi :: workspace index
* @param sd :: optional standard deviations of the parameters for calculating the error bars
* @return created workspase
*/
//boost::shared_ptr<API::MatrixWorkspace> IFunctionMW::createCalculatedWorkspace(
// boost::shared_ptr<const API::MatrixWorkspace> inWS,
// size_t wi,
// const std::vector<double>& sd
// )
//{
// const MantidVec& inputX = inWS->readX(wi);
// const MantidVec& inputY = inWS->readY(wi);
// const MantidVec& inputE = inWS->readE(wi);
// size_t nData = dataSize();
//
// size_t 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);
//
// assert(m_xMaxIndex-m_xMinIndex+1 == nData);
//
// for(size_t i=0;i<3;i++)
// {
// ws->dataX(i).assign(inputX.begin()+m_xMinIndex,inputX.begin()+m_xMaxIndex+1+histN);
// }
//
// ws->dataY(0).assign(inputY.begin()+m_xMinIndex,inputY.begin()+m_xMaxIndex+1);
// ws->dataE(0).assign(inputE.begin()+m_xMinIndex,inputE.begin()+m_xMaxIndex+1);
//
// MantidVec& Ycal = ws->dataY(1);
// MantidVec& Ecal = ws->dataE(1);
// MantidVec& E = ws->dataY(2);
//
// double* lOut = new double[nData]; // to capture output from call to function()
// function( lOut );
//
// for(size_t i=0; i<nData; i++)
// {
// Ycal[i] = lOut[i];
// E[i] = m_data[i] - Ycal[i];
// }
//
// delete [] lOut;
//
// if (sd.size() == static_cast<size_t>(this->nParams()))
// {
// SimpleJacobian J(nData,this->nParams());
// try
// {
// this->functionDeriv(&J);
// }
// catch(...)
// {
// this->calNumericalDeriv(&J,&m_xValues[0],nData);
// }
// for(size_t i=0; i<nData; i++)
// {
// double err = 0.0;
// for(size_t j=0;j< static_cast<size_t>(nParams());++j)
// {
// double d = J.get(i,j) * sd[j];
// err += d*d;
// }
// Ecal[i] = sqrt(err);
// }
// }
//
// return ws;
//}


} // namespace API
} // namespace Mantid

0 comments on commit aa9b247

Please sign in to comment.