Skip to content

Commit

Permalink
Refs #9782. Added Pearson-VII, changed AutoDiffTestAlg
Browse files Browse the repository at this point in the history
  • Loading branch information
Michael Wedel committed Aug 18, 2014
1 parent d96c282 commit 3b6f16b
Show file tree
Hide file tree
Showing 9 changed files with 259 additions and 49 deletions.
3 changes: 3 additions & 0 deletions Code/Mantid/Framework/CurveFitting/CMakeLists.txt
Expand Up @@ -67,6 +67,7 @@ set ( SRC_FILES
src/NormaliseByPeakArea.cpp
src/PRConjugateGradientMinimizer.cpp
src/ParDomain.cpp
src/PearsonVIIFamily.cpp
src/PlotPeakByLogValue.cpp
src/Polynomial.cpp
src/ProcessBackground.cpp
Expand Down Expand Up @@ -177,6 +178,7 @@ set ( INC_FILES
inc/MantidCurveFitting/NormaliseByPeakArea.h
inc/MantidCurveFitting/PRConjugateGradientMinimizer.h
inc/MantidCurveFitting/ParDomain.h
inc/MantidCurveFitting/PearsonVIIFamily.h
inc/MantidCurveFitting/PlotPeakByLogValue.h
inc/MantidCurveFitting/Polynomial.h
inc/MantidCurveFitting/ProcessBackground.h
Expand Down Expand Up @@ -276,6 +278,7 @@ set ( TEST_FILES
NeutronBk2BkExpConvPVoigtTest.h
NormaliseByPeakAreaTest.h
PRConjugateGradientTest.h
PearsonVIIFamilyTest.h
PlotPeakByLogValueTest.h
PolynomialTest.h
ProcessBackgroundTest.h
Expand Down
@@ -0,0 +1,85 @@
#ifndef MANTID_CURVEFITTING_PEARSONVIIFAMILY_H_
#define MANTID_CURVEFITTING_PEARSONVIIFAMILY_H_

#include "MantidKernel/System.h"
#include "MantidAPI/IFunction1D.h"
#include "MantidAPI/ParamFunction.h"
#include "MantidAPI/IFunction1DAutoDiff.h"


namespace Mantid
{
namespace CurveFitting
{

/** PearsonVIIFamily : TODO: DESCRIPTION
Copyright © 2014 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://github.com/mantidproject/mantid>
Code Documentation is available at: <http://doxygen.mantidproject.org>
*/

namespace Pearsons {

class DLLExport PearsonVIIHandCoded : public API::IFunction1D, public API::ParamFunction
{
public:
virtual ~PearsonVIIHandCoded() {}

std::string name() const { return "PearsonVIIHandCoded"; }

void function1D(double *out, const double *xValues, const size_t nData) const;
void functionDeriv1D(API::Jacobian *out, const double *xValues, const size_t nData);

protected:
void init();
};

class DLLExport PearsonVIINumDiff : public API::IFunction1D, public API::ParamFunction
{
public:
virtual ~PearsonVIINumDiff() {}

std::string name() const { return "PearsonVIINumDiff"; }

void function1D(double *out, const double *xValues, const size_t nData) const;

protected:
void init();
};

class DLLExport PearsonVIIAutoDiff : virtual public API::IFunction1DAutoDiff, virtual public API::ParamFunction
{
public:
virtual ~PearsonVIIAutoDiff() {}

std::string name() const { return "PearsonVIIAutoDiff"; }

void function1DAutoDiff(const API::FunctionDomain1D &domain, std::vector<adept::adouble> &y, const API::AutoDiffParameterAdapter &parameters) const;

protected:
void init();

};

} // namespace PearsonVIIFamily
} // namespace CurveFitting
} // namespace Mantid

#endif /* MANTID_CURVEFITTING_PEARSONVIIFAMILY_H_ */
62 changes: 24 additions & 38 deletions Code/Mantid/Framework/CurveFitting/src/AutoDiffTestAlg.cpp
Expand Up @@ -3,6 +3,7 @@
#include "MantidCurveFitting/GaussianAutoDiff.h"
#include "MantidCurveFitting/GaussianHandCoded.h"
#include "MantidCurveFitting/LorentzianFamily.h"
#include "MantidCurveFitting/PearsonVIIFamily.h"
#include "MantidAPI/FunctionDomain1D.h"
#include "MantidAPI/FunctionValues.h"
#include "MantidCurveFitting/Jacobian.h"
Expand Down Expand Up @@ -60,7 +61,7 @@ const std::string AutoDiffTestAlg::summary() const { return "Test"; }
void AutoDiffTestAlg::init()
{
declareProperty(new WorkspaceProperty<MatrixWorkspace>("InputWorkspace","",Direction::Input), "Data with 20 gaussian peaks.");
declareProperty("GaussianParameters", "", "Gaussian parameters (h, x0, s)");
declareProperty("GaussianParameters", "", "Function parameters");
declareProperty(new WorkspaceProperty<MatrixWorkspace>("OutputWorkspace","",Direction::Output), "Data with 20 gaussian peaks.");
declareProperty("DerivativeType","adept","How to calculate derivatives.");

Expand All @@ -77,14 +78,17 @@ namespace
{
if ( type == "adept" )
{
return IFunction_sptr(new Lorentzians::LorentzianAutoDiff);
//return IFunction_sptr(new GaussianAutoDiff);
//return IFunction_sptr(new Lorentzians::LorentzianAutoDiff);
return IFunction_sptr(new Pearsons::PearsonVIIAutoDiff);
return IFunction_sptr(new GaussianAutoDiff);
} else if ( type == "num") {
return IFunction_sptr(new Lorentzians::LorentzianNumDiff);
//return IFunction_sptr(new GaussianNumDiff);
//return IFunction_sptr(new Lorentzians::LorentzianNumDiff);
return IFunction_sptr(new Pearsons::PearsonVIINumDiff);
return IFunction_sptr(new GaussianNumDiff);
}
return IFunction_sptr(new Lorentzians::LorentzianHandCoded);
//return IFunction_sptr(new GaussianHandCoded);
//return IFunction_sptr(new Lorentzians::LorentzianHandCoded);
return IFunction_sptr(new Pearsons::PearsonVIIHandCoded);
return IFunction_sptr(new GaussianHandCoded);
}
}

Expand All @@ -101,10 +105,10 @@ void AutoDiffTestAlg::exec()
std::string parameterString = getProperty("GaussianParameters");
std::vector<double> paramValues = parameterValues(parameterString);

if(fitData->getNumberHistograms() > 4) {
IFunction_sptr f = getFunction(m_derType);
f->initialize();
IFunction_sptr f = getFunction(m_derType);
f->initialize();

if(fitData->getNumberHistograms() > (f->nParams() + 1)) {
for(size_t i = 0; i < paramValues.size(); ++i) {
f->setParameter(i, paramValues[i] - 0.001 * paramValues[i]);
}
Expand All @@ -115,7 +119,7 @@ void AutoDiffTestAlg::exec()
fitAlgorithm->setProperty("CalcErrors", true);
fitAlgorithm->setProperty("Function", f);
fitAlgorithm->setProperty("InputWorkspace", fitData);
fitAlgorithm->setProperty("WorkspaceIndex", 4);
fitAlgorithm->setProperty("WorkspaceIndex", static_cast<int>(f->nParams() + 1));

fitAlgorithm->execute();

Expand All @@ -134,7 +138,7 @@ void AutoDiffTestAlg::exec()

FunctionDomain1DVector x(fitData->readX(0));
FunctionValues y(x);
CurveFitting::Jacobian J(x.size(), 3);
CurveFitting::Jacobian J(x.size(), g->nParams());

Kernel::Timer timerF;

Expand Down Expand Up @@ -166,33 +170,15 @@ void AutoDiffTestAlg::exec()
}
}

MantidVec &dfdx0 = t->dataY(1);
const MantidVec &dfdx0r = fitData->readY(1);
for(size_t i = 0; i < x.size(); ++i) {
dfdx0[i] = 1.0 - J.get(i, 1) / dfdx0r[i];

if(std::isinf(dfdx0[i])) {
dfdx0[i] = 0.0;
}
}

MantidVec &dfdh = t->dataY(2);
const MantidVec &dfdhr = fitData->readY(2);
for(size_t i = 0; i < x.size(); ++i) {
dfdh[i] = 1.0 - J.get(i, 0) / dfdhr[i];

if(std::isinf(dfdh[i])) {
dfdh[i] = 0.0;
}
}

MantidVec &dfds = t->dataY(3);
const MantidVec &dfdsr = fitData->readY(3);
for(size_t i = 0; i < x.size(); ++i) {
dfds[i] = 1.0 - J.get(i, 2) / dfdsr[i];
for(size_t j = 1; (j < f->nParams() + 1); ++j) {
MantidVec &partialDeriv = t->dataY(j);
const MantidVec &reference = fitData->readY(j);
for(size_t i = 0; i < partialDeriv.size(); ++i) {
partialDeriv[i] = 1.0 - J.get(i, j - 1) / reference[i];

if(std::isinf(dfds[i])) {
dfds[i] = 0.0;
if(std::isinf(partialDeriv[i])) {
partialDeriv[i] = 0.0;
}
}
}

Expand Down
Expand Up @@ -30,8 +30,8 @@ void GaussianAutoDiff::function1DAutoDiff(const FunctionDomain1D &domain, std::v
}

void GaussianAutoDiff::init() {
declareParameter("Height");
declareParameter("PeakCentre");
declareParameter("Height");
declareParameter("Sigma");
}

Expand Down
6 changes: 3 additions & 3 deletions Code/Mantid/Framework/CurveFitting/src/GaussianHandCoded.cpp
Expand Up @@ -43,15 +43,15 @@ using namespace API;
double term = (xDiff)/sigma;
double expTerm = exp(-0.5 * term * term);

out->set(i, 0, expTerm);
out->set(i, 1, (height*expTerm*(xDiff))/(sigma * sigma));
out->set(i, 0, (height*expTerm*(xDiff))/(sigma * sigma));
out->set(i, 1, expTerm);
out->set(i, 2, (height*expTerm*term * term)/(sigma));
}
}

void GaussianHandCoded::init() {
declareParameter("Height");
declareParameter("PeakCentre");
declareParameter("Height");
declareParameter("Sigma");
}

Expand Down
2 changes: 1 addition & 1 deletion Code/Mantid/Framework/CurveFitting/src/GaussianNumDiff.cpp
Expand Up @@ -33,8 +33,8 @@ namespace CurveFitting
}

void GaussianNumDiff::init() {
declareParameter("Height");
declareParameter("PeakCentre");
declareParameter("Height");
declareParameter("Sigma");
}

Expand Down
12 changes: 6 additions & 6 deletions Code/Mantid/Framework/CurveFitting/src/LorentzianFamily.cpp
Expand Up @@ -32,16 +32,16 @@ void LorentzianHandCoded::functionDeriv1D(Jacobian *out, const double *xValues,
double num = diff*diff + s * s;
double ssquared = s * s;

out->set(i, 0, ssquared / (num));
out->set(i, 1, 2 * h * ssquared * diff / (num * num));
out->set(i, 0, 2 * h * ssquared * diff / (num * num));
out->set(i, 1, ssquared / (num));
out->set(i, 2, 2 * h * s / num * (1 - ssquared / num));
}
}

void LorentzianHandCoded::init()
{
declareParameter("Height");
declareParameter("Centre");
declareParameter("Height");
declareParameter("Gamma");
}

Expand All @@ -60,8 +60,8 @@ void LorentzianNumDiff::function1D(double *out, const double *xValues, const siz

void LorentzianNumDiff::init()
{
declareParameter("Height");
declareParameter("Centre");
declareParameter("Height");
declareParameter("Gamma");
}

Expand All @@ -71,7 +71,7 @@ void LorentzianAutoDiff::function1DAutoDiff(const FunctionDomain1D &domain, std:
adept::adouble x0 = parameters.getParameter("Centre");
adept::adouble s = parameters.getParameter("Gamma");

for(size_t i = 0; i < domain.size(); ++i) {
for(size_t i = 0; i < y.size(); ++i) {
adept::adouble diff = domain[i] - x0;
y[i] = h * s * s / (diff*diff + s * s);
}
Expand All @@ -80,8 +80,8 @@ void LorentzianAutoDiff::function1DAutoDiff(const FunctionDomain1D &domain, std:

void LorentzianAutoDiff::init()
{
declareParameter("Height");
declareParameter("Centre");
declareParameter("Height");
declareParameter("Gamma");
}

Expand Down

0 comments on commit 3b6f16b

Please sign in to comment.