Mantid
Loading...
Searching...
No Matches
GSLFunctions.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 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//----------------------------------------------------------------------
13
14namespace Mantid::CurveFitting {
15
22int gsl_f(const gsl_vector *x, void *params, gsl_vector *f) {
23 assert(x->data);
24 auto *p = reinterpret_cast<struct GSL_FitData *>(params);
25
26 // update function parameters
27 size_t ia = 0;
28 for (size_t i = 0; i < p->function->nParams(); ++i) {
29 if (p->function->isActive(i)) {
30 if (ia < x->size) {
31 p->function->setActiveParameter(i, x->data[ia]);
32 ++ia;
33 } else {
34 // The number of active parameters now exceeds the space
35 // originally allocated
37 }
38 }
39 }
40 p->function->applyTies();
41
42 auto values = std::dynamic_pointer_cast<API::FunctionValues>(p->costFunction->getValues());
43 if (!values) {
44 throw std::invalid_argument("FunctionValues expected");
45 }
46 p->function->function(*p->costFunction->getDomain(), *values);
47
48 // Add penalty
49 double penalty = 0.;
50 for (size_t i = 0; i < p->function->nParams(); ++i) {
51 API::IConstraint *c = p->function->getConstraint(i);
52 if (c) {
53 penalty += c->checkDeriv();
54 }
55 }
56
57 size_t n = values->size() - 1;
58 // add penalty to first and last point and every 10th point in between
59 if (penalty != 0.0) {
60 values->addToCalculated(0, penalty);
61 values->addToCalculated(n, penalty);
62
63 for (size_t i = 9; i < n; i += 10) {
64 values->addToCalculated(i, penalty);
65 }
66 }
67
68 // function() return calculated data values. Need to convert this values into
69 // calculated-observed devided by error values used by GSL
70
71 for (size_t i = 0; i < p->n; i++) {
72 f->data[i] = (values->getCalculated(i) - values->getFitData(i)) * values->getFitWeight(i);
73 }
74 return GSL_SUCCESS;
75}
76
84int gsl_df(const gsl_vector *x, void *params, gsl_matrix *J) {
85
86 auto *p = reinterpret_cast<struct GSL_FitData *>(params);
87 gsl_matrix *J_tr = gsl_matrix_calloc(J->size2, J->size1);
88 gsl_matrix_transpose_memcpy(J_tr, J);
89 EigenMatrix m(J_tr->size2, J_tr->size1);
90 p->J.setJ(&m);
91
92 // update function parameters
93 if (x->data) {
94 size_t ia = 0;
95 for (size_t i = 0; i < p->function->nParams(); ++i) {
96 if (p->function->isActive(i)) {
97 p->function->setActiveParameter(i, x->data[ia]);
98 ++ia;
99 }
100 }
101 }
102 p->function->applyTies();
103
104 // calculate the Jacobian
105 p->function->functionDeriv(*p->costFunction->getDomain(), p->J);
106
107 // p->function->addPenaltyDeriv(&p->J);
108 // add penalty
109 size_t n = p->costFunction->getValues()->size() - 1;
110 size_t ia = 0;
111 for (size_t i = 0; i < p->function->nParams(); ++i) {
112 if (!p->function->isActive(i))
113 continue;
114 API::IConstraint *c = p->function->getConstraint(i);
115 if (c) {
116 double penalty = c->checkDeriv2();
117 if (penalty != 0.0) {
118 double deriv = p->J.get(0, ia);
119 p->J.set(0, ia, deriv + penalty);
120 deriv = p->J.get(n, ia);
121 p->J.set(n, ia, deriv + penalty);
122
123 for (size_t j = 9; j < n; j += 10) {
124 deriv = p->J.get(j, ia);
125 p->J.set(j, ia, deriv + penalty);
126 }
127 }
128 } // if (c)
129 ++ia;
130 }
131
132 // functionDeriv() return derivatives of calculated data values. Need to
133 // convert this values into
134 // derivatives of calculated-observed divided by error values used by GSL
135 auto values = std::dynamic_pointer_cast<API::FunctionValues>(p->costFunction->getValues());
136 if (!values) {
137 throw std::invalid_argument("FunctionValues expected");
138 }
139
140 EigenMatrix m_tr = m.tr();
141 std::copy(&m_tr.mutator().data()[0], &m_tr.mutator().data()[J_tr->size1 * J_tr->size2], &J->data[0]);
142
143 for (size_t iY = 0; iY < p->n; iY++)
144 for (size_t iP = 0; iP < p->p; iP++) {
145 J->data[iY * p->p + iP] *= values->getFitWeight(iY);
146 }
147
148 return GSL_SUCCESS;
149}
150
158int gsl_fdf(const gsl_vector *x, void *params, gsl_vector *f, gsl_matrix *J) {
159 gsl_f(x, params, f);
160 gsl_df(x, params, J);
161 return GSL_SUCCESS;
162}
163
168GSL_FitData::GSL_FitData(const std::shared_ptr<CostFunctions::CostFuncLeastSquares> &cf)
169 : function(cf->getFittingFunction()), costFunction(cf) {
170 gsl_set_error_handler_off();
171 // number of active parameters
172 p = 0;
173 for (size_t i = 0; i < function->nParams(); ++i) {
174 if (function->isActive(i))
175 ++p;
176 }
177
178 // number of fitting data
179 n = cf->getValues()->size();
180
181 bool functionFixed = false;
182 if (p == 0) {
183 p = 1;
184 functionFixed = true;
185 }
186
187 // holdCalculatedJacobian = gsl_matrix_alloc (n, p);
188
189 initFuncParams = gsl_vector_alloc(p);
190
191 if (functionFixed) {
192 gsl_vector_set(initFuncParams, 0, 0.0);
193 } else {
194 size_t ia = 0;
195 for (size_t i = 0; i < function->nParams(); ++i) {
196 if (function->isActive(i)) {
197 gsl_vector_set(initFuncParams, ia, function->activeParameter(i));
198 ++ia;
199 }
200 }
201 }
202
203 int j = 0;
204 for (size_t i = 0; i < function->nParams(); ++i) {
205 if (function->isActive(i)) {
206 J.m_index.emplace_back(j);
207 j++;
208 } else
209 J.m_index.emplace_back(-1);
210 }
211}
212
214 // gsl_matrix_free(holdCalculatedJacobian);
215 gsl_vector_free(initFuncParams);
216}
217
218} // namespace Mantid::CurveFitting
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.
A wrapper around Eigen::Matrix.
Definition: EigenMatrix.h:33
map_type & mutator()
Get the map to Eigen matrix.
Definition: EigenMatrix.h:56
Exception thrown when a fitting function changes number of parameters during fit.
Definition: Exception.h:336
int gsl_fdf(const gsl_vector *x, void *params, gsl_vector *f, gsl_matrix *J)
Fit derivatives and function GSL wrapper.
int gsl_f(const gsl_vector *x, void *params, gsl_vector *f)
Fit GSL function wrapper.
int gsl_df(const gsl_vector *x, void *params, gsl_matrix *J)
Fit GSL derivative function wrapper.
Various GSL specific functions used GSL specific minimizers.
Definition: GSLFunctions.h:28
size_t n
number of points to be fitted (size of X, Y and sqrtWeightData arrays)
Definition: GSLFunctions.h:34
API::IFunction_sptr function
Pointer to the function.
Definition: GSLFunctions.h:38
size_t p
number of (active) fit parameters
Definition: GSLFunctions.h:36
JacobianImpl1< EigenMatrix > J
Jacobi matrix interface.
Definition: GSLFunctions.h:43
GSL_FitData(const std::shared_ptr< CostFunctions::CostFuncLeastSquares > &cf)
Constructor.
gsl_vector * initFuncParams
Initial function parameters.
Definition: GSLFunctions.h:41