Mantid
Loading...
Searching...
No Matches
BivariateNormal.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#include "MantidAPI/FunctionFactory.h"
8#include "MantidAPI/MatrixWorkspace.h"
9#include "MantidAPI/ParameterTie.h"
10
11#include "MantidCurveFitting/Constraints/BoundaryConstraint.h"
12#include "MantidCurveFitting/Functions/BivariateNormal.h"
13
14#include "MantidHistogramData/HistogramY.h"
15
16#include "MantidKernel/PhysicalConstants.h"
17#include "MantidKernel/System.h"
18
19#include <algorithm>
20#include <cmath>
21#include <cstdio>
22#include <fstream>
23#include <memory>
24#include <sstream>
25#include <string>
26
27using namespace Mantid::API;
28
29namespace Mantid::CurveFitting::Functions {
30
31using namespace CurveFitting;
32using namespace Constraints;
33using namespace HistogramData;
34
35namespace {
37Kernel::Logger g_log("BivariateNormal");
38} // namespace
39
40DECLARE_FUNCTION(BivariateNormal)
41
42// Indicies into Attrib array( local variable in initCommon
43#define S_int 0
44#define S_xint 1
45#define S_yint 2
46#define S_x2int 3
47#define S_y2int 4
48#define S_xyint 5
49#define S_y 6
50#define S_x 7
51#define S_x2 8
52#define S_y2 9
53#define S_xy 10
54#define S_1 11
55
56// Indicies into the LastParams array
57#define IBACK 0
58#define ITINTENS 1
59#define IXMEAN 2
60#define IYMEAN 3
61#define IVXX 4
62#define IVYY 5
63#define IVXY 6
64
65BivariateNormal::BivariateNormal()
66 : API::ParamFunction(), CalcVxx(false), CalcVyy(false), CalcVxy(false), NCells(0), CalcVariances(false), mIx(0.0),
67 mx(0.0), mIy(0.0), my(0.0), SIxx(0.0), SIyy(0.0), SIxy(0.0), Sxx(0.0), Syy(0.0), Sxy(0.0), TotI(0.0), TotN(0.0),
68 Varx0(-1.0), Vary0(-1.0), expVals(nullptr), uu(0.0), coefNorm(0.0), expCoeffx2(0.0), expCoeffy2(0.0),
69 expCoeffxy(0.0) {
70 LastParams[IVXX] = -1;
71 Varx0 = -1;
72}
73
74BivariateNormal::~BivariateNormal() { delete[] expVals; }
75
76// overwrite IFunction base class methods
77
78void BivariateNormal::function1D(double *out, const double *xValues, const size_t nData) const {
79
80 UNUSED_ARG(xValues);
81 UNUSED_ARG(nData);
82 if (nData == 0)
83 return;
84
85 if (Varx0 < 0) {
86 for (size_t i = 0; i < nData; i++)
87 out[i] = 0.0;
88 return;
89 }
90 double coefNorm, expCoeffx2, expCoeffy2, expCoeffxy, Varxx, Varxy, Varyy;
91 int NCells;
92 bool isNaNs;
93
94 API::MatrixWorkspace_const_sptr ws = getMatrixWorkspace();
95
96 const auto &D = ws->y(0);
97 const auto &X = ws->y(1);
98 const auto &Y = ws->y(2);
99 int K = 1;
100
101 if (nParams() > 4)
102 K = 3;
103
104 getConstraint(IBACK)->setPenaltyFactor(K * 3000);
105
106 double badParams = initCoeff(D, X, Y, coefNorm, expCoeffx2, expCoeffy2, expCoeffxy, NCells, Varxx, Varxy, Varyy);
107
108 std::ostringstream inf;
109 inf << "F Parameters=";
110 for (size_t k = 0; k < nParams(); k++)
111 inf << "," << getParameter(k);
112 if (nParams() < 6)
113 inf << "," << Varxx << "," << Varyy << "," << Varxy;
114 inf << '\n';
115
116 NCells = std::min<int>(static_cast<int>(nData), NCells);
117
118 double Background = getParameter(IBACK);
119 double Intensity = getParameter(ITINTENS);
120 double Xmean = getParameter(IXMEAN);
121 double Ymean = getParameter(IYMEAN);
122
123 double DDD = std::min<double>(10, 10 * std::max<double>(0, -Background));
124 int x = 0;
125 isNaNs = false;
126 double chiSq = 0;
127
128 // double penalty =0;
129
130 for (int i = 0; i < NCells; i++) {
131 // double pen =0;
132 if (badParams > 0)
133 out[x] = badParams;
134 else if (isNaNs)
135 out[x] = 10000;
136 else {
137 double dx = X[i] - Xmean;
138 double dy = Y[i] - Ymean;
139 out[x] =
140 Background + coefNorm * Intensity * exp(expCoeffx2 * dx * dx + expCoeffxy * dx * dy + expCoeffy2 * dy * dy);
141 out[x] = out[x] + DDD;
142
143 if (out[x] != out[x]) {
144 out[x] = 100000;
145 isNaNs = true;
146 }
147 }
148 double diff = out[x] - D[x];
149 chiSq += diff * diff;
150
151 x++;
152 }
153 inf << "Constr:";
154 for (size_t i = 0; i < nParams(); i++) {
155 IConstraint *constr = getConstraint(i);
156 if (constr)
157 inf << i << "=" << constr->check() << ";";
158 }
159 inf << "\n\n chiSq =" << chiSq << " nData " << nData << '\n';
160 g_log.debug(inf.str());
161}
162
163void BivariateNormal::functionDeriv1D(API::Jacobian *out, const double *xValues, const size_t nData) {
164 UNUSED_ARG(xValues);
165 UNUSED_ARG(nData);
166 if (nData == 0)
167 return;
168 double penDeriv = initCommon();
169
170 std::ostringstream inf;
171 inf << "***penalty(" << penDeriv << "),Parameters=";
172 for (size_t k = 0; k < 7; k++)
173 inf << "," << LastParams[k];
174 inf << '\n';
175 g_log.debug(inf.str());
176
177 std::vector<double> outf(nData, 0.0);
178 function1D(outf.data(), xValues, nData);
179
180 double uu = LastParams[IVXX] * LastParams[IVYY] - LastParams[IVXY] * LastParams[IVXY];
181 /* if( uu <=0)
182 {
183 penDeriv = -2*uu;
184 uu=1;
185 }
186 */
187 API::MatrixWorkspace_const_sptr ws = getMatrixWorkspace();
188 const auto &X = ws->y(1);
189 const auto &Y = ws->y(2);
190
191 for (int x = 0; x < NCells; x++) {
192
193 double penaltyDeriv = penDeriv;
194
195 double r = Y[x];
196 double c = X[x];
197
198 out->set(x, IBACK, +1.0);
199
200 if (penaltyDeriv <= 0)
201 out->set(x, ITINTENS, expVals[x] * coefNorm);
202 else if (LastParams[ITINTENS] < 0)
203 out->set(x, ITINTENS, -.01);
204 else
205 out->set(x, ITINTENS, 0.01);
206
207 double coefExp = coefNorm * LastParams[ITINTENS];
208
209 double coefxy = LastParams[IVXY] / uu;
210 double coefx2 = -LastParams[IVYY] / 2 / uu;
211
212 if (penaltyDeriv <= 0)
213 out->set(x, IXMEAN,
214 penaltyDeriv +
215 coefExp * expVals[x] * (-2 * coefx2 * (c - LastParams[IXMEAN]) - coefxy * (r - LastParams[IYMEAN])));
216 else // if(LastParams[IXMEAN] < 0)
217 out->set(x, IXMEAN, 0);
218 // else
219 // out->set(x,IXMEAN,0);
220
221 coefExp = coefNorm * LastParams[ITINTENS];
222
223 double coefy2 = -LastParams[IVXX] / 2 / uu;
224
225 if (penaltyDeriv <= 0)
226 out->set(x, IYMEAN,
227 penaltyDeriv +
228 coefExp * expVals[x] * (-coefxy * (c - LastParams[IXMEAN]) - 2 * coefy2 * (r - LastParams[IYMEAN])));
229 else // if(LastParams[IYMEAN] < 0)
230 out->set(x, IYMEAN, 0);
231 // else
232 // out->set(x,IYMEAN,0);
233
234 double M = 1;
235 if (nParams() < 5)
236 M = 10;
237 coefExp = coefNorm * LastParams[ITINTENS];
238
239 double C = -LastParams[IVYY] / 2 / uu;
240
241 double SIVXX;
242 if (penaltyDeriv <= 0)
243 SIVXX = coefExp * expVals[x] *
244 (C +
245 -LastParams[IVYY] / uu *
246 (coefx2 * (c - LastParams[IXMEAN]) * (c - LastParams[IXMEAN]) +
247 coefxy * (r - LastParams[IYMEAN]) * (c - LastParams[IXMEAN]) +
248 coefy2 * (r - LastParams[IYMEAN]) * (r - LastParams[IYMEAN])) -
249 (r - LastParams[IYMEAN]) * (r - LastParams[IYMEAN]) / 2 / uu);
250 else if (LastParams[IVXX] < .01)
251 SIVXX = -M;
252 else
253 SIVXX = 0;
254
255 if (LastParams[IVXX] > 1.2 * Varx0 && !CalcVariances)
256 SIVXX = 0;
257 if (LastParams[IVXX] < .8 * Varx0)
258 SIVXX = 0;
259 coefExp = coefNorm * LastParams[ITINTENS];
260
261 C = -LastParams[IVXX] / 2 / uu;
262
263 double SIVYY;
264 if (penaltyDeriv <= 0)
265 SIVYY = coefExp * expVals[x] *
266 (C +
267 -LastParams[IVXX] / uu *
268 (coefx2 * (c - LastParams[IXMEAN]) * (c - LastParams[IXMEAN]) +
269 coefxy * (r - LastParams[IYMEAN]) * (c - LastParams[IXMEAN]) +
270 coefy2 * (r - LastParams[IYMEAN]) * (r - LastParams[IYMEAN])) -
271 (c - LastParams[IXMEAN]) * (c - LastParams[IXMEAN]) / 2 / uu);
272
273 else if (LastParams[IVYY] < .01)
274 SIVYY = -M;
275 else
276 SIVYY = 0;
277
278 if (LastParams[IVYY] > 1.2 * Vary0 && !CalcVariances)
279 SIVYY = 0;
280
281 if (LastParams[IVYY] < .8 * Vary0)
282 SIVYY = 0;
283 coefExp = coefNorm * LastParams[ITINTENS];
284
285 C = LastParams[IVXY] / uu;
286
287 double SIVXY;
288 if (penaltyDeriv <= 0)
289 SIVXY = coefExp * expVals[x] *
290 (C +
291 2 * LastParams[IVXY] / uu *
292 (coefx2 * (c - LastParams[IXMEAN]) * (c - LastParams[IXMEAN]) +
293 coefxy * (r - LastParams[IYMEAN]) * (c - LastParams[IXMEAN]) +
294 coefy2 * (r - LastParams[IYMEAN]) * (r - LastParams[IYMEAN])) +
295 (r - LastParams[IYMEAN]) * (c - LastParams[IXMEAN]) / uu);
296
297 else if (uu < 0)
298 SIVXY = 2 * LastParams[IVXY];
299 else
300 SIVXY = 0;
301
302 if (!CalcVxx && nParams() > 6)
303 out->set(x, IVXX, penaltyDeriv + SIVXX);
304 else {
305 // out->set(x,IVXX,0.0);
306 double bdderiv = out->get(x, IBACK);
307
308 bdderiv += SIVXX *
309 (-Sxx - (LastParams[IXMEAN] - mx) * (LastParams[IXMEAN] - mx) * TotN + LastParams[IVXX] * TotN) /
310 (TotI - LastParams[IBACK] * TotN);
311
312 out->set(x, IBACK, bdderiv);
313
314 double mxderiv = out->get(x, IXMEAN);
315 mxderiv += SIVXX *
316 (2 * (LastParams[IXMEAN] - mIx) *
317
318 TotI -
319 2 * LastParams[IBACK] * (LastParams[IXMEAN] - mx) * TotN) /
320 (TotI - LastParams[IBACK] * TotN);
321 out->set(x, IXMEAN, mxderiv);
322 }
323 if (!CalcVyy && nParams() > 6)
324 out->set(x, IVYY, penaltyDeriv + SIVYY);
325 else {
326 // out->set(x,IVYY, 0.0);
327 double bdderiv = out->get(x, IBACK);
328
329 bdderiv += SIVYY *
330 (-Syy - (LastParams[IYMEAN] - my) * (LastParams[IYMEAN] - my) * TotN + LastParams[IVYY] * TotN) /
331 (TotI - LastParams[IBACK] * TotN);
332 out->set(x, IBACK, bdderiv);
333 double myderiv = out->get(x, IYMEAN);
334
335 myderiv += SIVYY *
336 (2 * (LastParams[IYMEAN] - mIy) *
337
338 TotI -
339 2 * LastParams[IBACK] * (LastParams[IYMEAN] - my) * TotN) /
340 (TotI - LastParams[IBACK] * TotN);
341 out->set(x, IYMEAN, myderiv);
342 }
343 if (!CalcVxy && nParams() > 6) {
344 out->set(x, IVXY, penaltyDeriv + SIVXY);
345
346 } else {
347 // out->set(x,IVXY, 0.0);
348 double bdderiv = out->get(x, IBACK);
349 bdderiv += SIVXY *
350 (-Sxy - (LastParams[IYMEAN] - my) * (LastParams[IXMEAN] - mx) * TotN + LastParams[IVXY] * TotN) /
351 (TotI - LastParams[IBACK] * TotN);
352 out->set(x, IBACK, bdderiv);
353 double myderiv = out->get(x, IYMEAN);
354 myderiv += SIVXY *
355 ((LastParams[IXMEAN] - mIx) *
356
357 TotI -
358 LastParams[IBACK] * (LastParams[IXMEAN] - mx) * TotN) /
359 (TotI - LastParams[IBACK] * TotN);
360 out->set(x, IYMEAN, myderiv);
361 double mxderiv = out->get(x, IXMEAN);
362 mxderiv += SIVXY * ((LastParams[IYMEAN] - mIy) * TotI - LastParams[IBACK] * (LastParams[IYMEAN] - my) * TotN) /
363 (TotI - LastParams[IBACK] * TotN);
364 out->set(x, IXMEAN, mxderiv);
365 }
366 }
367}
368
369void BivariateNormal::init() {
370 declareParameter("Background", 0.00);
371 declareParameter("Intensity", 0.00);
372 declareParameter("Mcol", 0.00, "Mean column(x) value");
373 declareParameter("Mrow", 0.00, "Mean row(y) value");
374
375 CalcVariances = false;
376
377 NCells = -1;
378 LastParams[IVXX] = -1;
379}
380
381double BivariateNormal::initCommon() {
382
383 double penalty = 0;
384 bool ParamsOK = true;
385 bool CommonsOK = true;
386 if (!expVals)
387 CommonsOK = false;
388
389 API::MatrixWorkspace_const_sptr ws = getMatrixWorkspace();
390 const auto &D = ws->y(0);
391 const auto &X = ws->y(1);
392 const auto &Y = ws->y(2);
393
394 if (NCells < 0) {
395 NCells = static_cast<int>(std::min<size_t>(D.size(), std::min<size_t>(X.size(), Y.size())));
396 CommonsOK = false;
397 }
398
399 double Attrib[12] = {0.0};
400
401 double MinX, MinY, MaxX, MaxY, MaxD, MinD;
402 MinX = MaxX = X[0];
403 MinY = MaxY = Y[0];
404 MaxD = MinD = D[0];
405
406 if (!CommonsOK) {
407
408 for (int i = 0; i < NCells; i++) {
409 Attrib[S_int] += D[i];
410 Attrib[S_xint] += D[i] * X[i];
411 Attrib[S_yint] += D[i] * Y[i];
412 Attrib[S_x2int] += D[i] * X[i] * X[i];
413 Attrib[S_y2int] += D[i] * Y[i] * Y[i];
414 Attrib[S_xyint] += D[i] * X[i] * Y[i];
415
416 Attrib[S_y] += Y[i];
417 Attrib[S_x] += X[i];
418 Attrib[S_x2] += X[i] * X[i];
419 Attrib[S_y2] += Y[i] * Y[i];
420 Attrib[S_xy] += X[i] * Y[i];
421 Attrib[S_1] += 1.0;
422
423 if (X[i] < MinX)
424 MinX = X[i];
425 if (X[i] > MaxX)
426 MaxX = X[i];
427 if (Y[i] < MinY)
428 MinY = Y[i];
429 if (Y[i] > MaxY)
430 MaxY = Y[i];
431 if (D[i] < MinD)
432 MinD = D[i];
433 if (D[i] > MaxD)
434 MaxD = D[i];
435 }
436
437 mIx = Attrib[S_xint] / Attrib[S_int];
438 mIy = Attrib[S_yint] / Attrib[S_int];
439 mx = Attrib[S_x] / Attrib[S_1];
440 my = Attrib[S_y] / Attrib[S_1];
441
442 SIxx = Attrib[S_x2int] - (Attrib[S_xint] * Attrib[S_xint]) / Attrib[S_int];
443 SIyy = Attrib[S_y2int] - (Attrib[S_yint]) * (Attrib[S_yint]) / Attrib[S_int];
444 SIxy = Attrib[S_xyint] - (Attrib[S_xint]) * (Attrib[S_yint]) / Attrib[S_int];
445
446 Sxx = Attrib[S_x2] - (Attrib[S_x]) * (Attrib[S_x]) / Attrib[S_1];
447 Syy = Attrib[S_y2] - (Attrib[S_y]) * (Attrib[S_y]) / Attrib[S_1];
448 Sxy = Attrib[S_xy] - (Attrib[S_x]) * (Attrib[S_y]) / Attrib[S_1];
449
450 // CommonsOK = false;
451
452 TotI = Attrib[S_int];
453 TotN = Attrib[S_1];
454
455 // CommonsOK = false;
456
457 if (getConstraint(0) == nullptr) {
458
459 addConstraint((std::make_unique<BoundaryConstraint>(this, "Background", 0, Attrib[S_int] / Attrib[S_1])));
460 }
461
462 double maxIntensity = Attrib[S_int] + 3 * sqrt(Attrib[S_int]);
463
464 if (maxIntensity < 100)
465 maxIntensity = 100;
466
467 if (getConstraint(1) == nullptr) {
468 addConstraint(std::make_unique<BoundaryConstraint>(this, "Intensity", 0, maxIntensity));
469 }
470
471 double minMeany = MinY * .9 + .1 * MaxY;
472 double maxMeany = MinY * .1 + .9 * MaxY;
473
474 if (getConstraint(3) == nullptr) {
475 addConstraint(std::make_unique<BoundaryConstraint>(this, "Mrow", minMeany, maxMeany));
476 }
477
478 double minMeanx = MinX * .9 + .1 * MaxX;
479 double maxMeanx = MinX * .1 + .9 * MaxX;
480 if (getConstraint(2) == nullptr) {
481 addConstraint(std::make_unique<BoundaryConstraint>(this, "Mcol", minMeanx, maxMeanx));
482 }
483
484 if (CalcVariances && nParams() > 6) {
485 std::ostringstream ssxx, ssyy, ssxy;
486
487 ssyy << std::string("(") << (SIyy) << "+(Mrow-" << (mIy) << ")*(Mrow-" << (mIy) << ")*" << Attrib[S_int]
488 << "-Background*" << (Syy) << "-Background*(Mrow-" << (my) << ")*(Mrow-" << (my) << ")*" << Attrib[S_1]
489 << ")/(" << (Attrib[S_int]) << "-Background*" << (Attrib[S_1]) << ")";
490
491 if (getTie(IVYY) == nullptr) {
492 tie("SSrow", ssyy.str());
493 CalcVxx = true;
494 }
495
496 ssxx << std::string("(") << (SIxx) << "+(Mcol-" << (mIx) << ")*(Mcol-" << (mIx) << ")*" << Attrib[S_int]
497 << "-Background*" << (Sxx) << "-Background*(Mcol-" << (mx) << ")*(Mcol-" << (mx) << ")*" << Attrib[S_1]
498 << ")/(" << (Attrib[S_int]) << "-Background*" << (Attrib[S_1]) << ")";
499
500 if (getTie(IVXX) == nullptr) {
501 tie("SScol", ssxx.str());
502 CalcVyy = true;
503 }
504
505 ssxy << std::string("(") << (SIxy) << "+(Mcol-" << (mIx) << ")*(Mrow-" << (mIy) << ")*" << Attrib[S_int]
506 << "-Background*" << (Sxy) << "-Background*(Mcol-" << (mx) << ")*(Mrow-" << (my) << ")*" << Attrib[S_1]
507 << ")/(" << (Attrib[S_int]) << "-Background*" << (Attrib[S_1]) << ")";
508
509 if (getTie(IVXY) == nullptr) {
510 tie("SSrc", ssxy.str());
511 CalcVxy = true;
512 }
513 }
514 CommonsOK = true;
515 }
516
517 if (LastParams[IVXX] < 0) {
518 ParamsOK = false;
519 CommonsOK = false;
520
521 } else
522 for (size_t i = 0; i < nParams() && ParamsOK; i++)
523 if (getParameter(i) != LastParams[i])
524 ParamsOK = false;
525
526 if (!ParamsOK) {
527
528 for (size_t i = 0; i < nParams(); i++)
529 LastParams[i] = getParameter(i);
530 }
531
532 if (!CommonsOK || !ParamsOK) {
533
534 int NCells1;
535 double Varxx, Varxy, Varyy;
536
537 Varxx = Varxy = Varyy = -1;
538 penalty = initCoeff(D, X, Y, coefNorm, expCoeffx2, expCoeffy2, expCoeffxy, NCells1, Varxx, Varxy, Varyy);
539
540 if (Varx0 < 0 && penalty <= 0) {
541 Varx0 = Varxx;
542 Vary0 = Varyy;
543 }
544
545 LastParams[IVXX] = Varxx;
546 LastParams[IVXY] = Varxy;
547 LastParams[IVYY] = Varyy;
548
549 delete[] expVals;
550 expVals = new double[NCells];
551
552 for (int i = 0; i < NCells; i++) {
553
554 double dx = X[i] - LastParams[IXMEAN];
555 double dy = Y[i] - LastParams[IYMEAN];
556 expVals[i] = exp(expCoeffx2 * dx * dx + expCoeffxy * dx * dy + expCoeffy2 * dy * dy);
557 }
558 }
559 return penalty;
560}
561
562double BivariateNormal::initCoeff(const HistogramY &D, const HistogramY &X, const HistogramY &Y, double &coefNorm,
563 double &expCoeffx2, double &expCoeffy2, double &expCoeffxy, int &NCells,
564 double &Varxx, double &Varxy, double &Varyy) const {
565
566 double Background = getParameter("Background");
567 bool zeroDenom = false;
568 if (TotI == 0 && TotN == 0)
569 zeroDenom = true;
570 else if (TotI - Background * TotN <= 0)
571 zeroDenom = true;
572 if (CalcVxx || nParams() < 6) {
573 Varxx =
574 (SIxx + (getParameter("Mcol") - mIx) * (getParameter("Mcol") - mIx) * TotI - getParameter("Background") * Sxx -
575 getParameter("Background") * (getParameter("Mcol") - (mx)) * (getParameter("Mcol") - (mx)) * TotN) /
576 (TotI - getParameter("Background") * TotN);
577
578 if (Varx0 > 0) {
579 Varxx = std::min<double>(Varxx, 1.21 * Varx0);
580 Varxx = std::max<double>(Varxx, .79 * Varx0);
581 }
582
583 } else {
584 Varxx = getParameter(IVXX);
585 }
586
587 double Mrow = getParameter("Mrow");
588
589 if (CalcVyy || nParams() < 6) {
590
591 Varyy = (SIyy + (Mrow - (mIy)) * (Mrow - (mIy)) * TotI -
592 getParameter("Background") * (Syy)-Background * (Mrow - (my)) * (Mrow - (my)) * TotN) /
593 (TotI - Background * TotN);
594 if (Vary0 > 0) {
595 Varyy = std::min<double>(Varyy, 1.21 * Vary0);
596 Varyy = std::max<double>(Varyy, .79 * Vary0);
597 }
598 } else {
599 Varyy = getParameter(IVYY);
600 }
601
602 if (CalcVxy || nParams() < 6) {
603 Varxy = ((SIxy) + (getParameter("Mcol") - (mIx)) * (getParameter("Mrow") - (mIy)) * TotI -
604 getParameter("Background") * (Sxy)-getParameter("Background") * (getParameter("Mcol") - (mx)) *
605 (getParameter("Mrow") - (my)) * TotN) /
606 (TotI - getParameter("Background") * TotN);
607
608 } else {
609 Varxy = getParameter(IVXY);
610 }
611
612 double uu = Varxx * Varyy - Varxy * Varxy;
613 double penalty;
614 if (zeroDenom)
615 penalty = 200;
616 else
617 penalty = std::max<double>(0, -Varxx + .01) + std::max<double>(0, -Varyy + .01) + std::max<double>(0, -uu + .01);
618
619 if (fabs(uu) < .01) {
620 if (uu < 0)
621 uu = -.01;
622 else
623 uu = .01;
624 }
625
626 NCells = static_cast<int>(std::min<size_t>(D.size(), std::min<size_t>(X.size(), Y.size())));
627 if (zeroDenom) {
628 coefNorm = expCoeffx2 = expCoeffy2 = 1;
629 expCoeffxy = 0;
630 Varxx = Varyy = 5;
631 Varxy = 0;
632 return penalty;
633 }
634 coefNorm = .5 / M_PI / sqrt(fabs(uu));
635
636 expCoeffx2 = -fabs(Varyy) / 2 / fabs(uu);
637 expCoeffxy = Varxy / uu;
638 expCoeffy2 = -fabs(Varxx) / 2 / fabs(uu);
639
640 if (nParams() < 5)
641 penalty *= 10;
642
643 return penalty;
644}
645
646} // namespace Mantid::CurveFitting::Functions
S_y2int
#define S_y2int
Definition: BivariateNormal.cpp:47
S_yint
#define S_yint
Definition: BivariateNormal.cpp:45
S_x2
#define S_x2
Definition: BivariateNormal.cpp:51
S_y
#define S_y
Definition: BivariateNormal.cpp:49
S_xint
#define S_xint
Definition: BivariateNormal.cpp:44
S_1
#define S_1
Definition: BivariateNormal.cpp:54
S_xy
#define S_xy
Definition: BivariateNormal.cpp:53
S_x2int
#define S_x2int
Definition: BivariateNormal.cpp:46
S_xyint
#define S_xyint
Definition: BivariateNormal.cpp:48
S_int
#define S_int
Definition: BivariateNormal.cpp:43
S_x
#define S_x
Definition: BivariateNormal.cpp:50
S_y2
#define S_y2
Definition: BivariateNormal.cpp:52
DECLARE_FUNCTION
#define DECLARE_FUNCTION(classname)
Macro for declaring a new type of function to be used with the FunctionFactory.
Definition: FunctionFactory.h:151
IVYY
#define IVYY
Definition: IntegratePeakTimeSlices.cpp:70
IBACK
#define IBACK
Definition: IntegratePeakTimeSlices.cpp:65
IYMEAN
#define IYMEAN
Definition: IntegratePeakTimeSlices.cpp:68
ITINTENS
#define ITINTENS
Definition: IntegratePeakTimeSlices.cpp:66
IXMEAN
#define IXMEAN
Definition: IntegratePeakTimeSlices.cpp:67
IVXX
#define IVXX
Definition: IntegratePeakTimeSlices.cpp:69
IVXY
#define IVXY
Definition: IntegratePeakTimeSlices.cpp:71
fabs
#define fabs(x)
Definition: Matrix.cpp:22
UNUSED_ARG
#define UNUSED_ARG(x)
Function arguments are sometimes unused in certain implmentations but are required for documentation ...
Definition: System.h:64
Mantid::API::IConstraint
An interface to a constraint.
Definition: IConstraint.h:26
Mantid::API::IConstraint::check
virtual double check()=0
Returns a penalty number which is bigger than or equal to zero If zero it means that the constraint i...
Mantid::API::IConstraint::setPenaltyFactor
virtual void setPenaltyFactor(const double &c)=0
set the penalty factor for the constraint Set panelty factor.
Mantid::API::IFunction1D::g_log
static Kernel::Logger g_log
Logger instance.
Definition: IFunction1D.h:74
Mantid::API::IFunctionMW::getMatrixWorkspace
std::shared_ptr< const API::MatrixWorkspace > getMatrixWorkspace() const
Get shared pointer to the workspace.
Definition: IFunctionMW.cpp:33
Mantid::API::IFunction::tie
virtual void tie(const std::string &parName, const std::string &expr, bool isDefault=false)
Tie a parameter to other parameters (or a constant)
Definition: IFunction.cpp:214
Mantid::API::IFunction::getTie
virtual ParameterTie * getTie(size_t i) const
Get the tie of i-th parameter.
Definition: IFunction.cpp:356
Mantid::API::IFunction::addConstraint
virtual void addConstraint(std::unique_ptr< IConstraint > ic)
Add a constraint to function.
Definition: IFunction.cpp:376
Mantid::API::IFunction::getConstraint
virtual IConstraint * getConstraint(size_t i) const
Get constraint of i-th parameter.
Definition: IFunction.cpp:392
Mantid::API::Jacobian
Represents the Jacobian in IFitFunction::functionDeriv.
Definition: Jacobian.h:22
Mantid::API::Jacobian::get
virtual double get(size_t iY, size_t iP)=0
Get the value to a Jacobian matrix element.
Mantid::API::Jacobian::set
virtual void set(size_t iY, size_t iP, double value)=0
Set a value to a Jacobian matrix element.
Mantid::API::ParamFunction
Implements the part of IFunction interface dealing with parameters.
Definition: ParamFunction.h:33
Mantid::API::ParamFunction::declareParameter
void declareParameter(const std::string &name, double initValue=0, const std::string &description="") override
Declare a new parameter.
Definition: ParamFunction.cpp:239
Mantid::API::ParamFunction::nParams
size_t nParams() const override
Total number of parameters.
Definition: ParamFunction.h:53
Mantid::API::ParamFunction::getParameter
double getParameter(size_t i) const override
Get i-th parameter.
Definition: ParamFunction.cpp:62
Mantid::CurveFitting::Functions::BivariateNormal::CalcVariances
bool CalcVariances
from experimental data versus fit the (Co)Variances
Definition: BivariateNormal.h:132
Mantid::CurveFitting::Functions::BivariateNormal::LastParams
double LastParams[9]
Saves previous/this set of parameters.
Definition: BivariateNormal.h:149
Mantid::CurveFitting::Functions::BivariateNormal::expVals
double * expVals
Save common exponential values for each cell.
Definition: BivariateNormal.h:151
Mantid::CurveFitting::Functions::BivariateNormal::function1D
void function1D(double *out, const double *xValues, const size_t nData) const override
Function you want to fit to.
Definition: BivariateNormal.cpp:78
Mantid::CurveFitting::Functions::BivariateNormal::functionDeriv1D
void functionDeriv1D(API::Jacobian *out, const double *xValues, const size_t nData) override
Derivatives of function with respect to active parameters.
Definition: BivariateNormal.cpp:163
Mantid::CurveFitting::Functions::BivariateNormal::initCommon
double initCommon()
Check for changes in parameters, etc.
Definition: BivariateNormal.cpp:381
Mantid::CurveFitting::Functions::BivariateNormal::init
void init() override
Function initialization. Declare function parameters in this method.
Definition: BivariateNormal.cpp:369
Mantid::CurveFitting::Functions::BivariateNormal::initCoeff
double initCoeff(const HistogramData::HistogramY &D, const HistogramData::HistogramY &X, const HistogramData::HistogramY &Y, double &coefNorm, double &expCoeffx2, double &expCoeffy2, double &expCoeffxy, int &NCells, double &Varxx, double &Varxy, double &Varyy) const
common values
Definition: BivariateNormal.cpp:562
Mantid::Kernel::Logger
The Logger class is in charge of the publishing messages from the framework through various channels.
Definition: Logger.h:52
Mantid::Kernel::Logger::debug
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
Mantid::API::g_log
Kernel::Logger g_log("ExperimentInfo")
static logger object
Mantid::API::MatrixWorkspace_const_sptr
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
Definition: MatrixWorkspace.h:522
Generated on Tue Jan 3 2023 11:57:33 for Mantid by doxygen 1.9.5