Mantid
Loading...
Searching...
No Matches
CostFuncFitting.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2019 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
7//----------------------------------------------------------------------
8// Includes
9//----------------------------------------------------------------------
16
17#include <gsl/gsl_multifit_nlin.h>
18#include <limits>
19#include <utility>
20
22
27 : m_dirtyVal(true), m_dirtyDeriv(true), m_dirtyHessian(true), m_includePenalty(true), m_value(0), m_pushed(false),
28 m_pushedValue(false) {}
29
34 m_dirtyVal = true;
35 m_dirtyDeriv = true;
36 m_dirtyHessian = true;
37}
38
42double CostFuncFitting::getParameter(size_t i) const {
44 return m_function->activeParameter(m_indexMap[i]);
45}
46
50void CostFuncFitting::setParameter(size_t i, const double &value) {
52 m_function->setActiveParameter(m_indexMap[i], value);
53 setDirty();
54}
55
61std::string CostFuncFitting::parameterName(size_t i) const {
63 return m_function->parameterName(m_indexMap[i]);
64}
65
69 return m_indexMap.size();
70}
71
82 m_function = std::move(function);
83 m_domain = std::move(domain);
84 m_values = std::move(values);
85 reset();
86}
87
93 return false;
94 }
95 if (m_function->nParams() != m_numberFunParams) {
96 reset();
97 }
98 auto nActive = m_indexMap.size();
99 for (size_t i = 0, j = 0; i < m_numberFunParams; ++i) {
100 if (m_function->isActive(i)) {
101 if (j >= nActive || m_indexMap[j] != i) {
102 reset();
103 break;
104 }
105 ++j;
106 }
107 }
108 return true;
109}
110
115 if (!isValid()) {
116 throw std::runtime_error("Fitting cost function isn't set");
117 }
118}
119
124 // construct the jacobian
125 EigenJacobian J(*m_function, m_values->size());
126 size_t na = this->nParams(); // number of active parameters
127 const auto &j = J.getJ();
128 assert(static_cast<size_t>(j.cols()) == na);
129 covar.resize(na, na);
130
131 // calculate the derivatives
132 m_function->functionDeriv(*m_domain, J);
133
134 // let the GSL to compute the covariance matrix
135 EigenMatrix covarTr(covar.tr());
136 EigenMatrix tempJTr;
137 tempJTr = j.transpose();
138 gsl_matrix_const_view tempJTrGSL = getGSLMatrixView_const(tempJTr.inspector());
139 gsl_matrix_view covarTrGSL = getGSLMatrixView(covarTr.mutator());
140
141 gsl_multifit_covar(&tempJTrGSL.matrix, epsrel, &covarTrGSL.matrix);
142 covar = covarTr.tr();
143}
144
152 EigenMatrix c;
153 calActiveCovarianceMatrix(c, epsrel);
154
155 size_t np = m_function->nParams();
156
157 bool isTransformationIdentity = true;
158 for (size_t i = 0; i < np; ++i) {
159 if (!m_function->isActive(i))
160 continue;
161 isTransformationIdentity =
162 isTransformationIdentity && (m_function->activeParameter(i) == m_function->getParameter(i));
163 }
164
165 if (isTransformationIdentity) {
166 // if the transformation is identity simply copy the matrix
167 covar = c;
168 } else {
169 // else do the transformation
170 EigenMatrix tm;
172 covar = tm.tr() * c * tm;
173 }
174}
175
182void CostFuncFitting::calFittingErrors(const EigenMatrix &covar, double chi2) {
184 size_t np = m_function->nParams();
185 auto covarMatrix = std::shared_ptr<Kernel::Matrix<double>>(new Kernel::Matrix<double>(np, np));
186 m_function->setCovarianceMatrix(covarMatrix);
187 size_t ia = 0;
188 for (size_t i = 0; i < np; ++i) {
189 if (!m_function->isActive(i)) {
190 m_function->setError(i, 0);
191 } else {
192 size_t ja = 0;
193 for (size_t j = 0; j < np; ++j) {
194 if (m_function->isActive(j)) {
195 (*covarMatrix)[i][j] = covar.get(ia, ja);
196 ++ja;
197
198 if (ja >= covar.size2())
199 break;
200 }
201 }
202 double err = sqrt(covar.get(ia, ia));
203 m_function->setError(i, err);
204 ++ia;
205
206 if (ia >= covar.size1())
207 break;
208 }
209 }
210 m_function->setCovarianceMatrix(covarMatrix);
211 m_function->setReducedChiSquared(chi2 / static_cast<double>((m_values->size() - np)));
212}
213
219 const double epsilon = std::numeric_limits<double>::epsilon() * 100;
220 size_t np = m_function->nParams();
221 size_t na = nParams();
222 tm.resize(na, na);
223 size_t ia = 0;
224 for (size_t i = 0; i < np; ++i) {
225 if (!m_function->isActive(i))
226 continue;
227 double p0 = m_function->getParameter(i);
228 for (size_t j = 0; j < na; ++j) {
229 double ap = getParameter(j);
230 double step = ap == 0.0 ? epsilon : ap * epsilon;
231 setParameter(j, ap + step);
232 tm.set(ia, j, (m_function->getParameter(i) - p0) / step);
233 setParameter(j, ap);
234 }
235 ++ia;
236 }
237}
238
241 if (m_function) {
242 m_function->applyTies();
243 }
244}
245
248 m_numberFunParams = m_function->nParams();
249 m_indexMap.clear();
250 for (size_t i = 0; i < m_numberFunParams; ++i) {
251 if (m_function->isActive(i)) {
252 m_indexMap.emplace_back(i);
253 }
254 API::IConstraint *c = m_function->getConstraint(i);
255 if (c) {
257 }
258 }
259 m_dirtyDeriv = true;
260 m_dirtyHessian = true;
261}
262
268 auto np = nParams();
269 if (np != params.size()) {
270 throw Kernel::Exception::FitSizeWarning(params.size(), np);
271 }
272 for (size_t i = 0; i < np; ++i) {
273 setParameter(i, params.get(i));
274 }
275 m_function->applyTies();
276}
277
283 auto np = nParams();
284 if (params.size() != np) {
285 params.resize(np);
286 }
287 for (size_t i = 0; i < np; ++i) {
288 params.set(i, getParameter(i));
289 }
290}
291
295double CostFuncFitting::val() const {
296 if (!m_dirtyVal) {
297 return m_value;
298 }
299
301
302 m_value = 0.0;
303
304 auto seqDomain = std::dynamic_pointer_cast<SeqDomain>(m_domain);
305
306 if (seqDomain) {
307 seqDomain->additiveCostFunctionVal(*this);
308 } else {
309 if (!m_values) {
310 throw std::runtime_error("CostFunction: undefined FunctionValues.");
311 }
313 }
314
315 // add penalty
316 if (m_includePenalty) {
317 for (size_t i = 0; i < m_function->nParams(); ++i) {
318 if (!m_function->isActive(i))
319 continue;
320 API::IConstraint *c = m_function->getConstraint(i);
321 if (c) {
322 m_value += c->check();
323 }
324 }
325 }
326
327 m_dirtyVal = false;
328 return m_value;
329}
330
334void CostFuncFitting::deriv(std::vector<double> &der) const {
335 valDerivHessian(true, false);
336
337 if (der.size() != nParams()) {
338 der.resize(nParams());
339 }
340 for (size_t i = 0; i < nParams(); ++i) {
341 der[i] = m_der.get(i);
342 }
343}
344
349double CostFuncFitting::valAndDeriv(std::vector<double> &der) const {
350 valDerivHessian(true, false);
351
352 if (der.size() != nParams()) {
353 der.resize(nParams());
354 }
355 for (size_t i = 0; i < nParams(); ++i) {
356 der[i] = m_der.get(i);
357 }
358 return m_value;
359}
360
366double CostFuncFitting::valDerivHessian(bool evalDeriv, bool evalHessian) const {
367 if (m_pushed || !evalDeriv) {
368 return val();
369 }
370
372 return m_value;
373 }
374
375 const size_t numParams = nParams();
376
378
379 m_value = 0.0;
380
381 m_der.resize(numParams);
382 m_der.zero();
383 if (evalHessian) {
384 m_hessian.resize(numParams, numParams);
385 m_hessian.zero();
386 }
387
388 auto seqDomain = std::dynamic_pointer_cast<SeqDomain>(m_domain);
389
390 if (seqDomain) {
391 seqDomain->additiveCostFunctionValDerivHessian(*this, evalDeriv, evalHessian);
392 } else {
393 if (!m_values) {
394 throw std::runtime_error("CostFunction: undefined FunctionValues.");
395 }
396 addValDerivHessian(m_function, m_domain, m_values, evalDeriv, evalHessian);
397 }
398
399 // Add constraints penalty
400 size_t np = m_function->nParams();
401 if (m_includePenalty) {
402 for (size_t i = 0; i < np; ++i) {
403 API::IConstraint *c = m_function->getConstraint(i);
404 if (c) {
405 m_value += c->check();
406 }
407 }
408 }
409 m_dirtyVal = false;
410
411 if (evalDeriv) {
412 if (m_includePenalty) {
413 size_t i = 0;
414 for (size_t ip = 0; ip < np; ++ip) {
415 if (!m_function->isActive(ip))
416 continue;
417 API::IConstraint *c = m_function->getConstraint(ip);
418 if (c) {
419 double d = m_der.get(i) + c->checkDeriv();
420 m_der.set(i, d);
421 }
422 ++i;
423 }
424 }
425 m_dirtyDeriv = false;
426
427 if (m_includePenalty) {
428 size_t i = 0;
429 for (size_t ip = 0; ip < np; ++ip) {
430 if (!m_function->isActive(ip))
431 continue;
432 API::IConstraint *c = m_function->getConstraint(ip);
433 if (c && !m_hessian.isEmpty()) {
434 double d = m_hessian.get(i, i) + c->checkDeriv2();
435 m_hessian.set(i, i, d);
436 }
437 ++i;
438 }
439 }
440 // clear the dirty flag if hessian was actually calculated
442 }
443
444 return m_value;
445}
446
451 if (m_pushed) {
452 return m_der;
453 }
456 }
457 return m_der;
458}
459
464 if (m_pushed) {
465 return m_hessian;
466 }
469 }
470 return m_hessian;
471}
472
477 if (m_pushed) {
478 throw std::runtime_error("Cost Function: double push.");
479 }
480 // make sure we are not dirty
483 m_pushed = true;
484}
485
490 if (!m_pushed) {
491 throw std::runtime_error("Cost Function: empty stack.");
492 }
495 m_pushed = false;
496 m_dirtyVal = false;
497 m_dirtyDeriv = false;
498 m_dirtyHessian = false;
499}
500
505 if (!m_pushed) {
506 throw std::runtime_error("Cost Function: empty stack.");
507 }
508 m_pushed = false;
509 setDirty();
510}
511
512} // namespace Mantid::CurveFitting::CostFunctions
const std::string & m_value
Definition: Algorithm.cpp:71
double value
The value of the point.
Definition: FitMW.cpp:51
An interface to a constraint.
Definition: IConstraint.h:26
virtual double checkDeriv()=0
Returns the derivative of the penalty for each active parameter.
virtual double checkDeriv2()=0
Returns the derivative of the penalty for each active parameter.
virtual double check()=0
Returns a penalty number which is bigger than or equal to zero If zero it means that the constraint i...
virtual void setParamToSatisfyConstraint()=0
Set the parameters of IFitFunction to satisfy constraint.
bool m_includePenalty
Flag to include constraint in cost function value.
API::FunctionValues_sptr m_values
Shared poinetr to the function values.
virtual void deriv(std::vector< double > &der) const override
Calculate the derivatives of the cost function.
void setParameter(size_t i, const double &value) override
Set i-th parameter.
virtual void calFittingErrors(const EigenMatrix &covar, double chi2)
Calculate fitting errors.
virtual double valDerivHessian(bool evalDeriv=true, bool evalHessian=true) const
Calculate the value, the first and the second derivatives of the cost function.
void pop()
Restore saved parameters, derivatives and hessian.
void reset() const
Reset the fitting function (neccessary if parameters get fixed/unfixed)
std::vector< size_t > m_indexMap
maps the cost function's parameters to the ones of the fitting function.
API::FunctionDomain_sptr m_domain
Shared pointer to the function domain.
void setParameters(const EigenVector &params)
Set all parameters.
void push()
Save current parameters, derivatives and hessian.
virtual void setFittingFunction(API::IFunction_sptr function, API::FunctionDomain_sptr domain, API::FunctionValues_sptr values)
Set fitting function.
void getParameters(EigenVector &params) const
Get all parameters into a GSLVector.
const EigenMatrix & getHessian() const
Return cached or calculate the Hessian.
virtual void addVal(API::FunctionDomain_sptr domain, API::FunctionValues_sptr values) const =0
Increment to the cost function by evaluating it on a domain.
size_t nParams() const override
Number of parameters.
void calTransformationMatrixNumerically(EigenMatrix &tm)
Calculate the transformation matrix T by numeric differentiation.
API::IFunction_sptr m_function
Shared pointer to the fitting function.
virtual double val() const override
Calculate value of cost function.
virtual double valAndDeriv(std::vector< double > &der) const override
Calculate the value and the derivatives of the cost function.
std::string parameterName(size_t i) const
Get name of i-th parameter.
bool isValid() const
Is the function set and valid?
virtual void calActiveCovarianceMatrix(EigenMatrix &covar, double epsrel=1e-8)
Calculates covariance matrix for fitting function's active parameters.
double getParameter(size_t i) const override
Get i-th parameter.
virtual void addValDerivHessian(API::IFunction_sptr function, API::FunctionDomain_sptr domain, API::FunctionValues_sptr values, bool evalDeriv=true, bool evalHessian=true) const =0
Increments the cost function and its derivatives by evaluating them on a domain.
size_t m_numberFunParams
Number of all parameters in the fitting function.
void drop()
Discard saved parameters, derivatives and hessian.
void checkValidity() const
Throw a runtime_error if function is invalid.
const EigenVector & getDeriv() const
Return cached or calculate the drivatives.
virtual void calCovarianceMatrix(EigenMatrix &covar, double epsrel=1e-8)
Calculates covariance matrix.
void applyTies()
Apply ties in the fitting function.
Two implementations of Jacobian.
Definition: EigenJacobian.h:23
map_type & getJ()
Get the map to the jacobian.
Definition: EigenJacobian.h:48
A wrapper around Eigen::Matrix.
Definition: EigenMatrix.h:33
bool isEmpty() const
Is matrix empty.
void set(size_t i, size_t j, double value)
Set an element.
double get(size_t i, size_t j) const
Get an element.
size_t size1() const
First size of the matrix.
EigenMatrix tr() const
Calculate the eigensystem of a symmetric matrix.
map_type & mutator()
Get the map to Eigen matrix.
Definition: EigenMatrix.h:56
void zero()
Set all elements to zero.
size_t size2() const
Second size of the matrix.
const map_type inspector() const
Get a const copy of the Eigen matrix.
Definition: EigenMatrix.h:58
void resize(const size_t nx, const size_t ny)
Resize the matrix.
A wrapper around Eigen::Vector.
Definition: EigenVector.h:27
void set(const size_t i, const double value)
Set an element.
double get(const size_t i) const
Get an element.
void resize(const size_t n)
Resize the vector.
Definition: EigenVector.cpp:96
size_t size() const
Size of the vector.
Exception thrown when a fitting function changes number of parameters during fit.
Definition: Exception.h:336
Numerical Matrix class.
Definition: Matrix.h:42
std::shared_ptr< FunctionValues > FunctionValues_sptr
typedef for a shared pointer
std::shared_ptr< IFunction > IFunction_sptr
shared pointer to the function base class
Definition: IFunction.h:732
std::shared_ptr< FunctionDomain > FunctionDomain_sptr
typedef for a shared pointer
gsl_matrix_const_view const getGSLMatrixView_const(const map_type m)
take data from a constEigen Matrix and return a transposed gsl view.
Definition: GSLFunctions.h:63
gsl_matrix_view getGSLMatrixView(map_type &tr)
take data from an Eigen Matrix and return a transposed a gsl view.
Definition: GSLFunctions.h:56