From c542f000b9c2e54b2861bd7fef3fcd53f514dc6c Mon Sep 17 00:00:00 2001
From: Roman Tolchenov <roman.tolchenov@stfc.ac.uk>
Date: Wed, 14 Mar 2012 15:50:21 +0000
Subject: [PATCH] Re #4158. Added MultiDomainFunction
---
Code/Mantid/Framework/API/CMakeLists.txt | 7 +
.../API/inc/MantidAPI/CompositeFunction.h | 21 +-
.../API/inc/MantidAPI/FunctionValues.h | 3 +-
.../Framework/API/inc/MantidAPI/JointDomain.h | 64 +++++
.../API/inc/MantidAPI/MultiDomainFunction.h | 80 ++++++
.../Framework/API/src/CompositeFunction.cpp | 2 -
.../Framework/API/src/FunctionValues.cpp | 17 ++
Code/Mantid/Framework/API/src/JointDomain.cpp | 35 +++
.../Framework/API/src/MultiDomainFunction.cpp | 174 +++++++++++++
.../API/test/MultiDomainFunctionTest.h | 245 ++++++++++++++++++
.../Framework/CurveFitting/CMakeLists.txt | 1 +
.../MantidCurveFitting/CostFuncLeastSquares.h | 9 +
.../inc/MantidCurveFitting/GSLVector.h | 6 +
.../CurveFitting/src/CostFuncLeastSquares.cpp | 178 +++++++++----
.../src/LevenbergMarquardtMDMinimizer.cpp | 1 +
.../test/MultiDomainFunctionTest.h | 140 ++++++++++
16 files changed, 919 insertions(+), 64 deletions(-)
create mode 100644 Code/Mantid/Framework/API/inc/MantidAPI/JointDomain.h
create mode 100644 Code/Mantid/Framework/API/inc/MantidAPI/MultiDomainFunction.h
create mode 100644 Code/Mantid/Framework/API/src/JointDomain.cpp
create mode 100644 Code/Mantid/Framework/API/src/MultiDomainFunction.cpp
create mode 100644 Code/Mantid/Framework/API/test/MultiDomainFunctionTest.h
create mode 100644 Code/Mantid/Framework/CurveFitting/test/MultiDomainFunctionTest.h
diff --git a/Code/Mantid/Framework/API/CMakeLists.txt b/Code/Mantid/Framework/API/CMakeLists.txt
index d600963e290..a664cac83c9 100644
--- a/Code/Mantid/Framework/API/CMakeLists.txt
+++ b/Code/Mantid/Framework/API/CMakeLists.txt
@@ -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
diff --git a/Code/Mantid/Framework/API/inc/MantidAPI/CompositeFunction.h b/Code/Mantid/Framework/API/inc/MantidAPI/CompositeFunction.h
index ab75819c402..d8b405259b5 100644
--- a/Code/Mantid/Framework/API/inc/MantidAPI/CompositeFunction.h
+++ b/Code/Mantid/Framework/API/inc/MantidAPI/CompositeFunction.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
diff --git a/Code/Mantid/Framework/API/inc/MantidAPI/FunctionValues.h b/Code/Mantid/Framework/API/inc/MantidAPI/FunctionValues.h
index 4e9373ec4da..07a1b40740e 100644
--- a/Code/Mantid/Framework/API/inc/MantidAPI/FunctionValues.h
+++ b/Code/Mantid/Framework/API/inc/MantidAPI/FunctionValues.h
@@ -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
diff --git a/Code/Mantid/Framework/API/inc/MantidAPI/JointDomain.h b/Code/Mantid/Framework/API/inc/MantidAPI/JointDomain.h
new file mode 100644
index 00000000000..cb6d03740e4
--- /dev/null
+++ b/Code/Mantid/Framework/API/inc/MantidAPI/JointDomain.h
@@ -0,0 +1,64 @@
+#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 © 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_*/
diff --git a/Code/Mantid/Framework/API/inc/MantidAPI/MultiDomainFunction.h b/Code/Mantid/Framework/API/inc/MantidAPI/MultiDomainFunction.h
new file mode 100644
index 00000000000..2607e1561f9
--- /dev/null
+++ b/Code/Mantid/Framework/API/inc/MantidAPI/MultiDomainFunction.h
@@ -0,0 +1,80 @@
+#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 © 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_*/
diff --git a/Code/Mantid/Framework/API/src/CompositeFunction.cpp b/Code/Mantid/Framework/API/src/CompositeFunction.cpp
index db71f7b0e22..831aa49c922 100644
--- a/Code/Mantid/Framework/API/src/CompositeFunction.cpp
+++ b/Code/Mantid/Framework/API/src/CompositeFunction.cpp
@@ -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);
}
diff --git a/Code/Mantid/Framework/API/src/FunctionValues.cpp b/Code/Mantid/Framework/API/src/FunctionValues.cpp
index 3ac0d60fc32..6936e349067 100644
--- a/Code/Mantid/Framework/API/src/FunctionValues.cpp
+++ b/Code/Mantid/Framework/API/src/FunctionValues.cpp
@@ -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.
diff --git a/Code/Mantid/Framework/API/src/JointDomain.cpp b/Code/Mantid/Framework/API/src/JointDomain.cpp
new file mode 100644
index 00000000000..d031a3f3678
--- /dev/null
+++ b/Code/Mantid/Framework/API/src/JointDomain.cpp
@@ -0,0 +1,35 @@
+//----------------------------------------------------------------------
+// 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
diff --git a/Code/Mantid/Framework/API/src/MultiDomainFunction.cpp b/Code/Mantid/Framework/API/src/MultiDomainFunction.cpp
new file mode 100644
index 00000000000..a2a4f5632dc
--- /dev/null
+++ b/Code/Mantid/Framework/API/src/MultiDomainFunction.cpp
@@ -0,0 +1,174 @@
+//----------------------------------------------------------------------
+// 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
diff --git a/Code/Mantid/Framework/API/test/MultiDomainFunctionTest.h b/Code/Mantid/Framework/API/test/MultiDomainFunctionTest.h
new file mode 100644
index 00000000000..a255a9911dc
--- /dev/null
+++ b/Code/Mantid/Framework/API/test/MultiDomainFunctionTest.h
@@ -0,0 +1,245 @@
+#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);
+ //ii.clear();
+ //ii.push_back(0);
+ //ii.push_back(2);
+ //multi.setDomainIndices(2,ii);
+
+ FunctionValues values(domain);
+ multi.function(domain,values);
+
+ const double A = multi.getFunction(0)->getParameter("A");
+ const double B = multi.getFunction(0)->getParameter("B");
+ auto d = static_cast<const FunctionDomain1D&>(domain.getDomain(0));
+ for(size_t i = 0; i < 9; ++i)
+ {
+ TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d[i]);
+ }
+ for(size_t i = 9; i < 19; ++i)
+ {
+ TS_ASSERT_EQUALS(values.getCalculated(i), 0);
+ }
+ for(size_t i = 19; i < 30; ++i)
+ {
+ TS_ASSERT_EQUALS(values.getCalculated(i), 0);
+ }
+
+ }
+
+ void test_calc_domain1_only()
+ {
+ multi.setDomainIndex(0,1);
+ multi.setDomainIndices(1,std::vector<size_t>());
+ multi.setDomainIndices(2,std::vector<size_t>());
+
+ FunctionValues values(domain);
+ multi.function(domain,values);
+
+ const double A = multi.getFunction(0)->getParameter("A");
+ const double B = multi.getFunction(0)->getParameter("B");
+ auto d = static_cast<const FunctionDomain1D&>(domain.getDomain(1));
+ for(size_t i = 0; i < 9; ++i)
+ {
+ TS_ASSERT_EQUALS(values.getCalculated(i), 0);
+ }
+ for(size_t i = 9; i < 19; ++i)
+ {
+ TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d[i-9]);
+ }
+ for(size_t i = 19; i < 30; ++i)
+ {
+ TS_ASSERT_EQUALS(values.getCalculated(i), 0);
+ }
+ }
+
+ void test_calc_domain2_only()
+ {
+ multi.setDomainIndex(0,2);
+ multi.setDomainIndices(1,std::vector<size_t>());
+ multi.setDomainIndices(2,std::vector<size_t>());
+
+ FunctionValues values(domain);
+ multi.function(domain,values);
+
+ const double A = multi.getFunction(0)->getParameter("A");
+ const double B = multi.getFunction(0)->getParameter("B");
+ auto d = static_cast<const FunctionDomain1D&>(domain.getDomain(2));
+ for(size_t i = 0; i < 9; ++i)
+ {
+ TS_ASSERT_EQUALS(values.getCalculated(i), 0);
+ }
+ for(size_t i = 9; i < 19; ++i)
+ {
+ TS_ASSERT_EQUALS(values.getCalculated(i), 0);
+ }
+ for(size_t i = 19; i < 30; ++i)
+ {
+ TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d[i-19]);
+ }
+
+ }
+
+ void test_calc_all_domains()
+ {
+ multi.clearDomainIndices();
+ multi.setDomainIndices(1,std::vector<size_t>());
+ multi.setDomainIndices(2,std::vector<size_t>());
+
+ FunctionValues values(domain);
+ multi.function(domain,values);
+
+ const double A = multi.getFunction(0)->getParameter("A");
+ const double B = multi.getFunction(0)->getParameter("B");
+ auto d0 = static_cast<const FunctionDomain1D&>(domain.getDomain(0));
+ for(size_t i = 0; i < 9; ++i)
+ {
+ TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d0[i]);
+ }
+ auto d1 = static_cast<const FunctionDomain1D&>(domain.getDomain(1));
+ for(size_t i = 9; i < 19; ++i)
+ {
+ TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d1[i-9]);
+ }
+ auto d2 = static_cast<const FunctionDomain1D&>(domain.getDomain(2));
+ for(size_t i = 19; i < 30; ++i)
+ {
+ TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d2[i-19]);
+ }
+
+ }
+
+ void test_calc()
+ {
+ multi.setDomainIndex(0,0);
+ std::vector<size_t> ii;
+ ii.push_back(0);
+ ii.push_back(1);
+ multi.setDomainIndices(1,ii);
+ ii.clear();
+ ii.push_back(0);
+ ii.push_back(2);
+ multi.setDomainIndices(2,ii);
+
+ FunctionValues values(domain);
+ multi.function(domain,values);
+
+ double A = multi.getFunction(0)->getParameter("A") +
+ multi.getFunction(1)->getParameter("A") +
+ multi.getFunction(2)->getParameter("A");
+ double B = multi.getFunction(0)->getParameter("B") +
+ multi.getFunction(1)->getParameter("B") +
+ multi.getFunction(2)->getParameter("B");
+ auto d0 = static_cast<const FunctionDomain1D&>(domain.getDomain(0));
+ for(size_t i = 0; i < 9; ++i)
+ {
+ TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d0[i]);
+ }
+
+ A = multi.getFunction(1)->getParameter("A");
+ B = multi.getFunction(1)->getParameter("B");
+ auto d1 = static_cast<const FunctionDomain1D&>(domain.getDomain(1));
+ for(size_t i = 9; i < 19; ++i)
+ {
+ TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d1[i-9]);
+ }
+
+ A = multi.getFunction(2)->getParameter("A");
+ B = multi.getFunction(2)->getParameter("B");
+ auto d2 = static_cast<const FunctionDomain1D&>(domain.getDomain(2));
+ for(size_t i = 19; i < 30; ++i)
+ {
+ TS_ASSERT_EQUALS(values.getCalculated(i), A + B * d2[i-19]);
+ }
+
+ }
+
+private:
+ MultiDomainFunction multi;
+ JointDomain domain;
+};
+
+#endif /*MULTIDOMAINFUNCTIONTEST_H_*/
diff --git a/Code/Mantid/Framework/CurveFitting/CMakeLists.txt b/Code/Mantid/Framework/CurveFitting/CMakeLists.txt
index 2be6cef1ee6..902104b9471 100644
--- a/Code/Mantid/Framework/CurveFitting/CMakeLists.txt
+++ b/Code/Mantid/Framework/CurveFitting/CMakeLists.txt
@@ -168,6 +168,7 @@ set ( TEST_FILES
test/LogNormalTest.h
test/Lorentzian1DTest.h
test/LorentzianTest.h
+ test/MultiDomainFunctionTest.h
# test/MultiBGTest.h
test/PlotPeakByLogValueTest.h
test/PRConjugateGradientTest.h
diff --git a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/CostFuncLeastSquares.h b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/CostFuncLeastSquares.h
index bb070e8d1ef..97495df7b40 100644
--- a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/CostFuncLeastSquares.h
+++ b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/CostFuncLeastSquares.h
@@ -77,6 +77,15 @@ protected:
virtual void calActiveCovarianceMatrix(GSLMatrix& covar, double epsrel = 1e-8);
+ void addVal(
+ API::FunctionDomain_sptr domain,
+ API::FunctionValues_sptr values
+ )const;
+ void addValDerivHessian(
+ API::FunctionDomain_sptr domain,
+ API::FunctionValues_sptr values,
+ bool evalFunction = true, bool evalDeriv = true, bool evalHessian = true) const;
+
private:
mutable double m_value;
diff --git a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/GSLVector.h b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/GSLVector.h
index cdb7d405ce4..c54cf4f6839 100644
--- a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/GSLVector.h
+++ b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/GSLVector.h
@@ -109,6 +109,12 @@ namespace Mantid
throw std::out_of_range("GSLVector index is out of range.");
}
+ // Set all elements to zero
+ void zero()
+ {
+ gsl_vector_set_zero( m_vector );
+ }
+
GSLVector& operator+=(const GSLVector& v)
{
gsl_vector_add(m_vector, v.gsl());
diff --git a/Code/Mantid/Framework/CurveFitting/src/CostFuncLeastSquares.cpp b/Code/Mantid/Framework/CurveFitting/src/CostFuncLeastSquares.cpp
index cd5f3369183..dfd85fc99f6 100644
--- a/Code/Mantid/Framework/CurveFitting/src/CostFuncLeastSquares.cpp
+++ b/Code/Mantid/Framework/CurveFitting/src/CostFuncLeastSquares.cpp
@@ -4,6 +4,7 @@
#include "MantidCurveFitting/CostFuncLeastSquares.h"
#include "MantidCurveFitting/Jacobian.h"
#include "MantidAPI/IConstraint.h"
+#include "MantidAPI/CompositeDomain.h"
namespace
{
@@ -24,34 +25,52 @@ double CostFuncLeastSquares::val() const
if ( !m_dirtyVal ) return m_value;
checkValidity();
- m_domain->reset();
- m_function->function(*m_domain,*m_values);
- size_t ny = m_values->size();
- double retVal = 0.0;
+ m_value = 0.0;
- for (size_t i = 0; i < ny; i++)
+ auto compDomain = boost::dynamic_pointer_cast<API::CompositeDomain>(m_domain);
+
+ //if (compDomain)
+ //{
+ //}
+ //else
{
- double val = ( m_values->getCalculated(i) - m_values->getFitData(i) ) * m_values->getFitWeight(i);
- retVal += val * val;
+ addVal(m_domain,m_values);
}
-
- retVal *= 0.5;
+ // add penalty
for(size_t i=0;i<m_function->nParams();++i)
{
if ( !m_function->isActive(i) ) continue;
API::IConstraint* c = m_function->getConstraint(i);
if (c)
{
- retVal += c->check();
+ m_value += c->check();
}
}
m_dirtyVal = false;
- return retVal;
+ return m_value;
}
+void CostFuncLeastSquares::addVal(API::FunctionDomain_sptr domain, API::FunctionValues_sptr values)const
+{
+ m_function->function(*m_domain,*m_values);
+ size_t ny = m_values->size();
+
+ double retVal = 0.0;
+
+ for (size_t i = 0; i < ny; i++)
+ {
+ double val = ( m_values->getCalculated(i) - m_values->getFitData(i) ) * m_values->getFitWeight(i);
+ retVal += val * val;
+ }
+
+ m_value += 0.5 * retVal;
+
+}
+
+
/// Calculate the derivatives of the cost function
/// @param der :: Container to output the derivatives
void CostFuncLeastSquares::deriv(std::vector<double>& der) const
@@ -104,21 +123,102 @@ double CostFuncLeastSquares::valDerivHessian(bool evalFunction, bool evalDeriv,
if (m_dirtyVal) evalFunction = true;
checkValidity();
+
+ if (evalFunction)
+ {
+ m_value = 0.0;
+ }
+ if (evalDeriv)
+ {
+ m_der.resize(nParams());
+ m_der.zero();
+ }
+ if (evalHessian)
+ {
+ m_hessian.resize(nParams(),nParams());
+ m_hessian.zero();
+ }
+
+ auto compDomain = boost::dynamic_pointer_cast<API::CompositeDomain>(m_domain);
+
+ //if (compDomain)
+ //{
+ //}
+ //else
+ {
+ addValDerivHessian(m_domain,m_values,evalFunction,evalDeriv,evalHessian);
+ }
+
+ // Add constraints penaly
+ size_t np = m_function->nParams();
+ if (evalFunction)
+ {
+ for(size_t i = 0; i < np; ++i)
+ {
+ API::IConstraint* c = m_function->getConstraint(i);
+ if (c)
+ {
+ m_value += c->check();
+ }
+ }
+ m_dirtyVal = false;
+ }
+
+ if (evalDeriv)
+ {
+ size_t i = 0;
+ for(size_t ip = 0; ip < np; ++ip)
+ {
+ API::IConstraint* c = m_function->getConstraint(ip);
+ if (c)
+ {
+ if ( !m_function->isActive(ip) ) continue;
+ double d = m_der.get(i) + c->checkDeriv();
+ m_der.set(i,d);
+ }
+ ++i;
+ }
+ m_dirtyDeriv = false;
+ }
+
+ if (evalDeriv)
+ {
+ size_t i = 0;
+ for(size_t ip = 0; ip < np; ++ip)
+ {
+ API::IConstraint* c = m_function->getConstraint(ip);
+ if (c)
+ {
+ if ( !m_function->isActive(ip) ) continue;
+ double d = m_hessian.get(i,i) + c->checkDeriv2();
+ m_hessian.set(i,i,d);
+ }
+ ++i;
+ }
+ m_dirtyHessian = false;
+ }
+
+ return m_value;
+}
+
+/**
+ * Update the cost function, derivatives and hessian by adding values calculated
+ * on a domain.
+ */
+void CostFuncLeastSquares::addValDerivHessian(
+ API::FunctionDomain_sptr domain,
+ API::FunctionValues_sptr values,
+ bool evalFunction , bool evalDeriv, bool evalHessian) const
+{
size_t np = m_function->nParams(); // number of parameters
size_t ny = m_domain->size(); // number of data points
Jacobian jacobian(ny,np);
if (evalFunction)
{
- m_domain->reset();
- m_function->function(*m_domain,*m_values);
+ m_function->function(*domain,*values);
}
- m_domain->reset();
- m_function->functionDeriv(*m_domain,jacobian);
+ m_function->functionDeriv(*domain,jacobian);
- if (m_der.size() != nParams())
- {
- m_der.resize(nParams());
- }
size_t iActiveP = 0;
double fVal = 0.0;
if (debug)
@@ -137,7 +237,7 @@ double CostFuncLeastSquares::valDerivHessian(bool evalFunction, bool evalDeriv,
for(size_t ip = 0; ip < np; ++ip)
{
if ( !m_function->isActive(ip) ) continue;
- double d = 0.0;
+ double d = m_der.get(iActiveP);
for(size_t i = 0; i < ny; ++i)
{
double calc = m_values->getCalculated(i);
@@ -150,11 +250,6 @@ double CostFuncLeastSquares::valDerivHessian(bool evalFunction, bool evalDeriv,
fVal += y * y;
}
}
- API::IConstraint* c = m_function->getConstraint(ip);
- if (c)
- {
- d += c->checkDeriv();
- }
m_der.set(iActiveP, d);
//std::cerr << "der " << ip << ' ' << der[iActiveP] << std::endl;
++iActiveP;
@@ -162,27 +257,12 @@ double CostFuncLeastSquares::valDerivHessian(bool evalFunction, bool evalDeriv,
if (evalFunction)
{
- fVal *= 0.5;
- for(size_t i = 0; i < np; ++i)
- {
- API::IConstraint* c = m_function->getConstraint(i);
- if (c)
- {
- fVal += c->check();
- }
- }
- m_value = fVal;
- m_dirtyVal = false;
+ m_value += 0.5 * fVal;
}
- m_dirtyDeriv = false;
- if (!evalHessian) return m_value;
+ if (!evalHessian) return;
size_t na = m_der.size(); // number of active parameters
- if (m_hessian.size1() != na || m_hessian.size2() != na)
- {
- m_hessian.resize(na,na);
- }
size_t i1 = 0; // active parameter index
for(size_t i = 0; i < np; ++i) // over parameters
@@ -192,20 +272,12 @@ double CostFuncLeastSquares::valDerivHessian(bool evalFunction, bool evalDeriv,
for(size_t j = 0; j <= i; ++j) // over ~ half of parameters
{
if ( !m_function->isActive(j) ) continue;
- double d = 0.0;
+ double d = m_hessian.get(i1,i2);
for(size_t k = 0; k < ny; ++k) // over fitting data
{
double w = m_values->getFitWeight(k);
d += jacobian.get(k,i) * jacobian.get(k,j) * w * w;
}
- if (i == j)
- {
- API::IConstraint* c = m_function->getConstraint(i);
- if (c)
- {
- d += c->checkDeriv2();
- }
- }
m_hessian.set(i1,i2,d);
//std::cerr << "hess " << i1 << ' ' << i2 << std::endl;
if (i1 != i2)
@@ -216,11 +288,9 @@ double CostFuncLeastSquares::valDerivHessian(bool evalFunction, bool evalDeriv,
}
++i1;
}
-
- m_dirtyHessian = false;
- return m_value;
}
+
const GSLVector& CostFuncLeastSquares::getDeriv() const
{
if (m_pushed)
diff --git a/Code/Mantid/Framework/CurveFitting/src/LevenbergMarquardtMDMinimizer.cpp b/Code/Mantid/Framework/CurveFitting/src/LevenbergMarquardtMDMinimizer.cpp
index 0300d837aa0..56bbbe4834b 100644
--- a/Code/Mantid/Framework/CurveFitting/src/LevenbergMarquardtMDMinimizer.cpp
+++ b/Code/Mantid/Framework/CurveFitting/src/LevenbergMarquardtMDMinimizer.cpp
@@ -33,6 +33,7 @@ m_mu(0),
m_nu(2.0),
m_rho(1.0)
{
+ gsl_set_error_handler_off();
}
/// Initialize minimizer, i.e. pass a function to minimize.
diff --git a/Code/Mantid/Framework/CurveFitting/test/MultiDomainFunctionTest.h b/Code/Mantid/Framework/CurveFitting/test/MultiDomainFunctionTest.h
new file mode 100644
index 00000000000..546bb7091fc
--- /dev/null
+++ b/Code/Mantid/Framework/CurveFitting/test/MultiDomainFunctionTest.h
@@ -0,0 +1,140 @@
+#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 "MantidCurveFitting/CostFuncLeastSquares.h"
+#include "MantidCurveFitting/LevenbergMarquardtMDMinimizer.h"
+
+#include <cxxtest/TestSuite.h>
+#include <boost/make_shared.hpp>
+#include <algorithm>
+#include <iostream>
+
+using namespace Mantid;
+using namespace Mantid::API;
+using namespace Mantid::CurveFitting;
+
+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,1,x);
+ out->set(i,0,1.0);
+ }
+ }
+};
+
+class MultiDomainFunctionTest : public CxxTest::TestSuite
+{
+public:
+ MultiDomainFunctionTest()
+ {
+ multi = boost::make_shared<MultiDomainFunction>();
+ 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",0);
+
+ multi->getFunction(1)->setParameter("A",0);
+ multi->getFunction(1)->setParameter("B",0);
+
+ multi->getFunction(2)->setParameter("A",0);
+ multi->getFunction(2)->setParameter("B",0);
+
+ domain = boost::make_shared<JointDomain>();
+ 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_fit()
+ {
+
+ //system("pause");
+ auto values = boost::make_shared<FunctionValues>(*domain);
+ const double A0 = 0, A1 = 1, A2 = 2;
+ const double B0 = 1, B1 = 2, B2 = 3;
+
+ auto d0 = static_cast<const FunctionDomain1D&>(domain->getDomain(0));
+ for(size_t i = 0; i < d0.size(); ++i)
+ {
+ values->setFitData(i, A0 + A1 + A2 + (B0 + B1 + B2) * d0[i]);
+ }
+
+ auto d1 = static_cast<const FunctionDomain1D&>(domain->getDomain(1));
+ for(size_t i = 0; i < d1.size(); ++i)
+ {
+ values->setFitData(9 + i, A0 + A1 + (B0 + B1) * d1[i]);
+ }
+
+ auto d2 = static_cast<const FunctionDomain1D&>(domain->getDomain(2));
+ for(size_t i = 0; i < d2.size(); ++i)
+ {
+ values->setFitData(19 + i, A0 + A2 + (B0 + B2) * d2[i]);
+ }
+ values->setFitWeights(1);
+
+ multi->clearDomainIndices();
+ std::vector<size_t> ii(2);
+ ii[0] = 0;
+ ii[1] = 1;
+ multi->setDomainIndices(1,ii);
+ ii[0] = 0;
+ ii[1] = 2;
+ multi->setDomainIndices(2,ii);
+
+ boost::shared_ptr<CostFuncLeastSquares> costFun(new CostFuncLeastSquares);
+ costFun->setFittingFunction(multi,domain,values);
+ TS_ASSERT_EQUALS(costFun->nParams(),6);
+
+ LevenbergMarquardtMDMinimizer s;
+ s.initialize(costFun);
+ TS_ASSERT(s.minimize());
+
+ TS_ASSERT_EQUALS(s.getError(),"success");
+ TS_ASSERT_DELTA(s.costFunctionVal(),0,1e-4);
+
+ TS_ASSERT_DELTA(multi->getFunction(0)->getParameter("A"),0,1e-8);
+ TS_ASSERT_DELTA(multi->getFunction(0)->getParameter("B"),1,1e-8);
+ TS_ASSERT_DELTA(multi->getFunction(1)->getParameter("A"),1,1e-8);
+ TS_ASSERT_DELTA(multi->getFunction(1)->getParameter("B"),2,1e-8);
+ TS_ASSERT_DELTA(multi->getFunction(2)->getParameter("A"),2,1e-8);
+ TS_ASSERT_DELTA(multi->getFunction(2)->getParameter("B"),3,1e-8);
+ }
+
+private:
+ boost::shared_ptr<MultiDomainFunction> multi;
+ boost::shared_ptr<JointDomain> domain;
+};
+
+#endif /*MULTIDOMAINFUNCTIONTEST_H_*/
--
GitLab