Skip to content

Commit

Permalink
Re #9782. Refined first draft, added tests
Browse files Browse the repository at this point in the history
  • Loading branch information
Michael Wedel committed Jul 1, 2014
1 parent f3703ba commit 9ca66ea
Show file tree
Hide file tree
Showing 3 changed files with 313 additions and 66 deletions.
59 changes: 39 additions & 20 deletions Code/Mantid/Framework/API/inc/MantidAPI/IFunctionAutoDiff.h
View file
Expand Up @@ -55,11 +55,11 @@ class AutoDiffJacobian : public Jacobian
}

virtual void set(size_t iY, size_t iP, double value) {
m_jacobian[iY + iP * m_nParams] = value;
m_jacobian[iY + iP * m_nValues] = value;
}

virtual double get(size_t iY, size_t iP) {
return m_jacobian[iY + iP * m_nParams];
return m_jacobian[iY + iP * m_nValues];
}

void setValid() {
Expand All @@ -82,35 +82,54 @@ class AutoDiffJacobian : public Jacobian
bool m_isValid;
};


class DLLExport IFunctionAutoDiff : virtual public IFunction, virtual public ParamFunction
class AutoDiffParameters
{
public:
IFunctionAutoDiff();
virtual ~IFunctionAutoDiff();
AutoDiffParameters(const IFunction *fromFunction)
: m_names(), m_parameters()
{
size_t np = fromFunction->nParams();

void function(const FunctionDomain &domain, FunctionValues &values) const;
void functionDeriv(const FunctionDomain &domain, Jacobian &jacobian);
m_parameters.resize(np);

virtual void functionAutoDiff(std::vector<adept::adouble> &x, std::vector<adept::adouble> &y, std::vector<adept::adouble> &parameters) const = 0;
for(size_t i = 0; i < np; ++i) {
m_names[fromFunction->parameterName(i)] = i;
m_parameters[i].set_value(fromFunction->getParameter(i));
}
}

virtual void setParameter(size_t i, const double &value, bool explicitlySet) {
ParamFunction::setParameter(i, value, explicitlySet);
size_t size() const { return m_parameters.size(); }

m_jacobian->setInvalid();
}
adept::adouble &operator[](size_t i) { return m_parameters[i]; }

virtual void setParameter(const std::string& name, const double &value, bool explicitlySet) {
ParamFunction::setParameter(name, value, explicitlySet);
adept::adouble getParameter(size_t i) const
{
return m_parameters[i];
}

m_jacobian->setInvalid();
adept::adouble getParameter(std::string &name) const
{
return m_parameters[m_names.at(name)];
}

protected:
std::vector<adept::adouble> getAdeptParameters() const;
private:
std::map<std::string, size_t> m_names;
std::vector<adept::adouble> m_parameters;
};


class DLLExport IFunctionAutoDiff : virtual public IFunction
{
public:
IFunctionAutoDiff();
virtual ~IFunctionAutoDiff();

void function(const FunctionDomain &domain, FunctionValues &values) const;
void functionDeriv(const FunctionDomain &domain, Jacobian &jacobian);

void functionWrapper(const FunctionDomain &domain, FunctionValues &values, AutoDiffJacobian *jacobian) const;

adept::Stack *m_stack;
AutoDiffJacobian *m_jacobian;
virtual void functionAutoDiff(std::vector<adept::adouble> &x, std::vector<adept::adouble> &y, const AutoDiffParameters &parameters) const = 0;
};


Expand Down
79 changes: 42 additions & 37 deletions Code/Mantid/Framework/API/src/IFunctionAutoDiff.cpp
View file
Expand Up @@ -6,15 +6,11 @@ namespace Mantid
namespace API
{


//----------------------------------------------------------------------------------------------
/** Constructor
*/
IFunctionAutoDiff::IFunctionAutoDiff() :
IFunction(),
ParamFunction(),
m_stack(new adept::Stack),
m_jacobian(new AutoDiffJacobian)
IFunction()
{
}

Expand All @@ -25,59 +21,68 @@ IFunctionAutoDiff::~IFunctionAutoDiff()
{
}

/**
*
*
* @param domain
* @param values
*/
void IFunctionAutoDiff::function(const FunctionDomain &domain, FunctionValues &values) const
{
const FunctionDomain1D &d1d = dynamic_cast<const FunctionDomain1D &>(domain);

const double *domainStart = d1d.getPointerAt(0);
size_t domainSize = d1d.size();

std::vector<adept::adouble> x(domainStart, domainStart + domainSize);
std::vector<adept::adouble> y(d1d.size());
std::vector<adept::adouble> parameters = getAdeptParameters();

m_stack->new_recording();

functionAutoDiff(x, y, parameters);

for(size_t i = 0; i < values.size(); ++i) {
*(values.getPointerToCalculated(i)) = y[i].value();
}

m_stack->independent(&parameters[0], parameters.size());
m_stack->dependent(&y[0], domainSize);
m_stack->jacobian(m_jacobian->rawValues());
functionWrapper(domain, values, 0);
}

void IFunctionAutoDiff::functionDeriv(const FunctionDomain &domain, Jacobian &jacobian)
{
if(!m_jacobian->isValid()) {
FunctionValues values(domain);
function(domain, values);
}
AutoDiffJacobian localJacobian;

FunctionValues values(domain);
functionWrapper(domain, values, &localJacobian);

size_t np = nParams();
size_t ny = domain.size();

for(size_t p = 0; p < np; ++p) {
for(size_t y = 0; y < ny; ++y) {
jacobian.set(y, p, m_jacobian->get(y, p));
jacobian.set(y, p, localJacobian.get(y, p));
}
}
}

std::vector<adept::adouble> IFunctionAutoDiff::getAdeptParameters() const
void IFunctionAutoDiff::functionWrapper(const FunctionDomain &domain, FunctionValues &values, AutoDiffJacobian *jacobian) const
{
size_t np = nParams();
const FunctionDomain1D &d1d = dynamic_cast<const FunctionDomain1D &>(domain);

const double *domainStart = d1d.getPointerAt(0);
size_t domainSize = d1d.size();

adept::Stack stack;

std::vector<adept::adouble> adeptParameters(np);
for(size_t i = 0; i < np; ++np) {
adeptParameters[i].set_value(getParameter(i));
std::vector<adept::adouble> x(domainStart, domainStart + domainSize);
std::vector<adept::adouble> y(d1d.size());

if(jacobian != 0) {
stack.new_recording();
} else {
stack.pause_recording();
}

return adeptParameters;
}
AutoDiffParameters params(this);

functionAutoDiff(x, y, params);

for(size_t i = 0; i < values.size(); ++i) {
*(values.getPointerToCalculated(i)) = y[i].value();
}

if(jacobian != 0) {
stack.independent(&params[0], params.size());
stack.dependent(&y[0], domainSize);

jacobian->setSize(params.size(), domainSize);
stack.jacobian(jacobian->rawValues());
}
}

} // namespace API
} // namespace Mantid

0 comments on commit 9ca66ea

Please sign in to comment.