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

Re #4158. Added MultiDomainFunction

parent 12892a02
......@@ -53,6 +53,7 @@ set ( SRC_FILES
src/ImplictFunctionFactory.cpp
src/InstrumentDataService.cpp
src/LiveListenerFactory.cpp
src/JointDomain.cpp
src/LoadAlgorithmFactory.cpp
src/LocatedDataRef.cpp
src/LocatedDataValue.cpp
......@@ -63,6 +64,8 @@ set ( SRC_FILES
src/MemoryManager.cpp
src/MultipleExperimentInfos.cpp
src/MultipleFileProperty.cpp
src/MultiDomainFunction.cpp
src/NumericAxis.cpp
src/NullCoordTransform.cpp
src/NumericAxis.cpp
src/ParamFunction.cpp
......@@ -170,6 +173,7 @@ set ( INC_FILES
inc/MantidAPI/InstrumentDataService.h
inc/MantidAPI/Jacobian.h
inc/MantidAPI/LiveListenerFactory.h
inc/MantidAPI/JointDomain.h
inc/MantidAPI/LoadAlgorithmFactory.h
inc/MantidAPI/LocatedDataRef.h
inc/MantidAPI/LocatedDataValue.h
......@@ -180,6 +184,8 @@ set ( INC_FILES
inc/MantidAPI/MemoryManager.h
inc/MantidAPI/MultipleExperimentInfos.h
inc/MantidAPI/MultipleFileProperty.h
inc/MantidAPI/MultiDomainFunction.h
inc/MantidAPI/NumericAxis.h
inc/MantidAPI/NullCoordTransform.h
inc/MantidAPI/NumericAxis.h
inc/MantidAPI/ParamFunction.h
......@@ -247,6 +253,7 @@ set ( TEST_FILES
test/MDGeometryTest.h
test/MatrixWorkspaceMDIteratorTest.h
test/MemoryManagerTest.h
test/MultiDomainFunctionTest.h
test/MultipleExperimentInfosTest.h
test/MultipleFilePropertyTest.h
test/NumericAxisTest.h
......
......@@ -197,15 +197,21 @@ typedef boost::shared_ptr<const CompositeFunction> CompositeFunction_const_sptr;
class PartialJacobian: public Jacobian
{
Jacobian* m_J; ///< pointer to the overall Jacobian
size_t m_iP0; ///< offset in the overall Jacobian for a particular function
//size_t m_iaP0; ///< offset in the active Jacobian for a particular function
size_t m_iY0; ///< fitting data index offset in the overall Jacobian for a particular function
size_t m_iP0; ///< parameter index offset in the overall Jacobian for a particular function
public:
/** Constructor
* @param J :: A pointer to the overall Jacobian
* @param iP0 :: The parameter index (declared) offset for a particular function
* @param iap0 :: The active parameter index (declared) offset for a particular function
*/
PartialJacobian(Jacobian* J,size_t iP0):m_J(J),m_iP0(iP0)//,m_iaP0(iap0)
PartialJacobian(Jacobian* J,size_t iP0):m_J(J),m_iY0(0),m_iP0(iP0)//,m_iaP0(iap0)
{}
/** Constructor
* @param J :: A pointer to the overall Jacobian
* @param iY0 :: The data index offset for a particular function
* @param iP0 :: The parameter index offset for a particular function
*/
PartialJacobian(Jacobian* J,size_t iY0,size_t iP0):m_J(J),m_iY0(iY0),m_iP0(iP0)
{}
/**
* Overridden Jacobian::set(...).
......@@ -215,7 +221,7 @@ public:
*/
void set(size_t iY, size_t iP, double value)
{
m_J->set(iY,m_iP0 + iP,value);
m_J->set(m_iY0 + iY,m_iP0 + iP,value);
}
/**
* Overridden Jacobian::get(...).
......@@ -224,11 +230,11 @@ public:
*/
double get(size_t iY, size_t iP)
{
return m_J->get(iY,m_iP0 + iP);
return m_J->get(m_iY0 + iY,m_iP0 + iP);
}
/** Add number to all iY (data) Jacobian elements for a given iP (parameter)
* @param value :: Value to add
* @param iActiveP :: The index of an active parameter.
* @param iP :: The index of an active parameter.
*/
virtual void addNumberToColumn(const double& value, const size_t& iP)
{
......@@ -236,6 +242,7 @@ public:
}
};
} // namespace API
} // namespace Mantid
......
......@@ -53,7 +53,7 @@ public:
FunctionValues(const FunctionValues& values);
/// Reset the values to match a new domain.
void reset(const FunctionDomain& domain);
/// Return the number of points, values, etc in the domain
/// Return the number of values
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;}
......@@ -63,6 +63,7 @@ public:
double getCalculated(size_t i) const {return m_calculated[i];}
double operator[](size_t i) const {return m_calculated[i];}
void addToCalculated(size_t i, double value) {m_calculated[i] += value;}
void addToCalculated(size_t i, const FunctionValues& values);
/// Get a pointer to calculated data at index i
double* getPointerToCalculated(size_t i);
/// Add other calculated values
......
#ifndef MANTID_API_JOINTDOMAIN_H_
#define MANTID_API_JOINTDOMAIN_H_
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAPI/DllConfig.h"
#include "MantidAPI/CompositeDomain.h"
#include <stdexcept>
#include <vector>
#include <algorithm>
namespace Mantid
{
namespace API
{
/** Base class that represents the domain of a function.
A domain is a generalisation of x (argument) and y (value) arrays.
A domain consists at least of a list of function arguments for which a function should
be evaluated and a buffer for the calculated values. If used in fitting also contains
the fit data and weights.
@author Roman Tolchenov, Tessella plc
@date 15/11/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 MANTID_API_DLL JointDomain: public CompositeDomain
{
public:
/// Return the number of points in the domain
virtual size_t size() const;
/// Return the number of parts in the domain
virtual size_t getNParts() const;
/// Return i-th domain
virtual const FunctionDomain& getDomain(size_t i) const;
void addDomain(FunctionDomain_sptr domain);
protected:
std::vector< FunctionDomain_sptr > m_domains;
};
} // namespace API
} // namespace Mantid
#endif /*MANTID_API_JOINTDOMAIN_H_*/
#ifndef MANTID_API_MULTIDOMAINFUNCTION_H_
#define MANTID_API_MULTIDOMAINFUNCTION_H_
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAPI/CompositeFunction.h"
#include <map>
namespace Mantid
{
namespace API
{
class CompositeDomain;
/** A composite function defined on a CompositeDomain. Member functions can be applied to
one or more member domains of the CompositeDomain. If two functions applied to the same domain
the results are added (+).
@author Roman Tolchenov, Tessella plc
@date 13/03/2012
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 MANTID_API_DLL MultiDomainFunction : public CompositeFunction
{
public:
/// Constructor
MultiDomainFunction():m_nDomains(0),m_maxIndex(0){}
/// Function you want to fit to.
/// @param domain :: The buffer for writing the calculated values. Must be big enough to accept dataSize() values
virtual void function(const FunctionDomain& domain, FunctionValues& values)const;
/// Derivatives of function with respect to active parameters
virtual void functionDeriv(const FunctionDomain& domain, Jacobian& jacobian);
/// Associate a function and a domain
void setDomainIndex(size_t funIndex, size_t domainIndex);
/// Associate a function and a list of domains
void setDomainIndices(size_t funIndex, const std::vector<size_t>& domainIndices);
/// Clear all domain indices
void clearDomainIndices();
protected:
/// Counts number of the domains
void countNumberOfDomains();
void countValueOffsets(const CompositeDomain& domain)const;
/// Domain index map: finction -> domain
std::map<size_t, std::vector<size_t> > m_domains;
/// Number of domains this MultiDomainFunction operates on. == number of different values in m_domains
size_t m_nDomains;
/// Maximum domain index
size_t m_maxIndex;
mutable std::vector<size_t> m_valueOffsets;
};
} // namespace API
} // namespace Mantid
#endif /*MANTID_API_MULTIDOMAINFUNCTION_H_*/
......@@ -126,7 +126,6 @@ void CompositeFunction::function(const FunctionDomain& domain, FunctionValues& v
values.zeroCalculated();
for(size_t iFun = 0; iFun < nFunctions(); ++iFun)
{
domain.reset();
m_functions[ iFun ]->function(domain,tmp);
values += tmp;
}
......@@ -141,7 +140,6 @@ void CompositeFunction::functionDeriv(const FunctionDomain& domain, Jacobian& ja
{
for(size_t iFun = 0; iFun < nFunctions(); ++iFun)
{
domain.reset();
PartialJacobian J(&jacobian,paramOffset(iFun));
getFunction(iFun)->functionDeriv(domain,J);
}
......
......@@ -93,6 +93,23 @@ namespace API
setCalculated(0.0);
}
/**
* Add values starting at index i.
*/
void FunctionValues::addToCalculated(size_t i, const FunctionValues& values)
{
if (i + values.size() > size())
{
throw std::runtime_error("Cannot add function values: different sizes.");
}
std::transform(
m_calculated.begin() + i,
m_calculated.begin() + i + values.size(),
values.m_calculated.begin(),
m_calculated.begin() + i,
std::plus<double>()
);
}
/**
* Set a fitting data value.
......
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAPI/JointDomain.h"
namespace Mantid
{
namespace API
{
/// Return the number of points in the domain
size_t JointDomain::size() const
{
size_t n = 0;
std::for_each(m_domains.begin(),m_domains.end(),[&n](const FunctionDomain_sptr d){
n += d->size();
});
return n;
}
/// Return the number of parts in the domain
size_t JointDomain::getNParts() const
{
return m_domains.size();
}
/// Return i-th domain
const FunctionDomain& JointDomain::getDomain(size_t i) const
{
return *m_domains.at(i);
}
void JointDomain::addDomain(FunctionDomain_sptr domain)
{
m_domains.push_back(domain);
}
} // namespace API
} // namespace Mantid
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAPI/MultiDomainFunction.h"
#include "MantidAPI/CompositeDomain.h"
#include <boost/lexical_cast.hpp>
#include <set>
namespace Mantid
{
namespace API
{
/**
* Associate a member function and a domain. The function will only be applied
* to this domain.
* @param funIndex :: Index of a member function.
* @param domainIndex :: Index of a domain to be associated with the function.
*/
void MultiDomainFunction::setDomainIndex(size_t funIndex, size_t domainIndex)
{
m_domains[funIndex] = std::vector<size_t>(1,domainIndex);
countNumberOfDomains();
}
/**
* Associate a member function and a list of domains. The function will only be applied
* to the listed domains.
* @param funIndex :: Index of a member function.
* @param domainIndex :: A vector with indices of domains to be associated with the function.
*/
void MultiDomainFunction::setDomainIndices(size_t funIndex, const std::vector<size_t>& domainIndices)
{
m_domains[funIndex] = domainIndices;
countNumberOfDomains();
}
/**
* Counts number of the domains.
*/
void MultiDomainFunction::countNumberOfDomains()
{
std::set<size_t> dSet;
for(auto it = m_domains.begin(); it != m_domains.end(); ++it)
{
if (it->second.size())
{
dSet.insert(it->second.begin(),it->second.end());
}
}
m_nDomains = dSet.size();
m_maxIndex = dSet.empty() ? 0 : *dSet.rbegin();
}
/**
* Count value offsets for each member domain in a CompositeDomain.
*/
void MultiDomainFunction::countValueOffsets(const CompositeDomain& domain)const
{
m_valueOffsets.clear();
m_valueOffsets.push_back(0);
for(size_t i = 0; i < domain.getNParts(); ++i)
{
const FunctionDomain& d = domain.getDomain(i);
m_valueOffsets.push_back(m_valueOffsets.back() + d.size());
}
}
/// Clear all domain indices
void MultiDomainFunction::clearDomainIndices()
{
m_domains.clear();
countNumberOfDomains();
}
/// Function you want to fit to.
/// @param domain :: The buffer for writing the calculated values. Must be big enough to accept dataSize() values
void MultiDomainFunction::function(const FunctionDomain& domain, FunctionValues& values)const
{
// works only on CompositeDomain
if (!dynamic_cast<const CompositeDomain*>(&domain))
{
throw std::invalid_argument("Non-CompositeDomain passed to MultiDomainFunction.");
}
const CompositeDomain& cd = dynamic_cast<const CompositeDomain&>(domain);
// domain must not have less parts than m_maxIndex
if (cd.getNParts() < m_maxIndex)
{
throw std::invalid_argument("CompositeDomain has too few parts ("
+ boost::lexical_cast<std::string>(cd.getNParts()) +
") for MultiDomainFunction (max index " + boost::lexical_cast<std::string>(m_maxIndex) + ").");
}
// domain and values must be consistent
if (cd.size() != values.size())
{
throw std::invalid_argument("MultiDomainFunction: domain and values have different sizes.");
}
// evaluate member functions
values.zeroCalculated();
for(size_t iFun = 0; iFun < nFunctions(); ++iFun)
{
// find the domains member function must be applied to
std::vector<size_t> domains;
auto it = m_domains.find(iFun);
if (it == m_domains.end())
{// apply to all domains
domains.resize(cd.getNParts());
size_t i = 0;
std::generate(domains.begin(),domains.end(),[&i]()->size_t{return i++;});
}
else
{// apply to selected domains
domains.assign(it->second.begin(),it->second.end());
}
countValueOffsets(cd);
for(auto i = domains.begin(); i != domains.end(); ++i)
{
const FunctionDomain& d = cd.getDomain(*i);
FunctionValues tmp(d);
getFunction(iFun)->function( d, tmp );
values.addToCalculated(m_valueOffsets[*i],tmp);
}
}
}
/// Derivatives of function with respect to active parameters
void MultiDomainFunction::functionDeriv(const FunctionDomain& domain, Jacobian& jacobian)
{
// works only on CompositeDomain
if (!dynamic_cast<const CompositeDomain*>(&domain))
{
throw std::invalid_argument("Non-CompositeDomain passed to MultiDomainFunction.");
}
const CompositeDomain& cd = dynamic_cast<const CompositeDomain&>(domain);
// domain must not have less parts than m_maxIndex
if (cd.getNParts() < m_maxIndex)
{
throw std::invalid_argument("CompositeDomain has too few parts ("
+ boost::lexical_cast<std::string>(cd.getNParts()) +
") for MultiDomainFunction (max index " + boost::lexical_cast<std::string>(m_maxIndex) + ").");
}
// evaluate member functions derivatives
for(size_t iFun = 0; iFun < nFunctions(); ++iFun)
{
// find the domains member function must be applied to
std::vector<size_t> domains;
auto it = m_domains.find(iFun);
if (it == m_domains.end())
{// apply to all domains
domains.resize(cd.getNParts());
size_t i = 0;
std::generate(domains.begin(),domains.end(),[&i]()->size_t{return i++;});
}
else
{// apply to selected domains
domains.assign(it->second.begin(),it->second.end());
}
countValueOffsets(cd);
for(auto i = domains.begin(); i != domains.end(); ++i)
{
const FunctionDomain& d = cd.getDomain(*i);
PartialJacobian J(&jacobian,m_valueOffsets[*i],paramOffset(iFun));
getFunction(iFun)->functionDeriv(d,J);
}
}
}
} // namespace API
} // namespace Mantid
#ifndef MULTIDOMAINFUNCTIONTEST_H_
#define MULTIDOMAINFUNCTIONTEST_H_
#include "MantidAPI/FunctionDomain1D.h"
#include "MantidAPI/FunctionValues.h"
#include "MantidAPI/MultiDomainFunction.h"
#include "MantidAPI/JointDomain.h"
#include "MantidAPI/IFunction1D.h"
#include "MantidAPI/ParamFunction.h"
#include <cxxtest/TestSuite.h>
#include <boost/make_shared.hpp>
#include <algorithm>
#include <iostream>
using namespace Mantid;
using namespace Mantid::API;
class MultiDomainFunctionTest_Function: public virtual IFunction1D, public virtual ParamFunction
{
public:
MultiDomainFunctionTest_Function():IFunction1D(),ParamFunction()
{
this->declareParameter("A",0);
this->declareParameter("B",0);
}
virtual std::string name() const {return "MultiDomainFunctionTest_Function";}
protected:
virtual void function1D(double* out, const double* xValues, const size_t nData)const
{
const double A = getParameter(0);
const double B = getParameter(1);
for(size_t i = 0; i < nData; ++i)
{
double x = xValues[i];
out[i] = A + B * x;
}
}
virtual void functionDeriv1D(Jacobian* out, const double* xValues, const size_t nData)
{
for(size_t i = 0; i < nData; ++i)
{
double x = xValues[i];
out->set(i,0,x);
out->set(i,1,1.0);
}
}
};
class MultiDomainFunctionTest : public CxxTest::TestSuite
{
public:
MultiDomainFunctionTest()
{
multi.addFunction(boost::make_shared<MultiDomainFunctionTest_Function>());
multi.addFunction(boost::make_shared<MultiDomainFunctionTest_Function>());
multi.addFunction(boost::make_shared<MultiDomainFunctionTest_Function>());
multi.getFunction(0)->setParameter("A",0);
multi.getFunction(0)->setParameter("B",1);
multi.getFunction(1)->setParameter("A",1);
multi.getFunction(1)->setParameter("B",2);
multi.getFunction(2)->setParameter("A",2);
multi.getFunction(2)->setParameter("B",3);
domain.addDomain(boost::make_shared<FunctionDomain1D>(0,1,9));
domain.addDomain(boost::make_shared<FunctionDomain1D>(1,2,10));
domain.addDomain(boost::make_shared<FunctionDomain1D>(2,3,11));
}
void test_calc_domain0_only()
{
multi.setDomainIndex(0,0);
multi.setDomainIndices(1,std::vector<size_t>());
multi.setDomainIndices(2,std::vector<size_t>());
//multi.setDomainIndex(1,1);
//multi.setDomainIndex(2,2);
//std::vector<size_t> ii;
//ii.push_back(0);
//ii.push_back(1);
//multi.setDomainIndices(1,ii);