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 "MantidCurveFitting/Functions/BivariateNormal.h"
8
9#include "MantidAPI/FunctionFactory.h"
10
11#include "MantidAPI/MatrixWorkspace.h"
12#include "MantidAPI/ParameterTie.h"
13
14#include "MantidCurveFitting/Constraints/BoundaryConstraint.h"
15#include "MantidHistogramData/HistogramY.h"
16
17#include "MantidKernel/PhysicalConstants.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 coefficientNorm, exponentialCoeffx2, exponentialCoeffy2, exponentialCoeffxy, Varxx, Varxy, Varyy;
91 int numberOfCells;
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, coefficientNorm, exponentialCoeffx2, exponentialCoeffy2, exponentialCoeffxy,
107 numberOfCells, Varxx, Varxy, Varyy);
108
109 std::ostringstream inf;
110 inf << "F Parameters=";
111 for (size_t k = 0; k < nParams(); k++)
112 inf << "," << getParameter(k);
113 if (nParams() < 6)
114 inf << "," << Varxx << "," << Varyy << "," << Varxy;
115 inf << '\n';
116
117 numberOfCells = std::min<int>(static_cast<int>(nData), numberOfCells);
118
119 double Background = getParameter(IBACK);
120 double Intensity = getParameter(ITINTENS);
121 double Xmean = getParameter(IXMEAN);
122 double Ymean = getParameter(IYMEAN);
123
124 double DDD = std::min<double>(10, 10 * std::max<double>(0, -Background));
125 int x = 0;
126 isNaNs = false;
127 double chiSq = 0;
128
129 // double penalty =0;
130
131 for (int i = 0; i < numberOfCells; i++) {
132 // double pen =0;
133 if (badParams > 0)
134 out[x] = badParams;
135 else if (isNaNs)
136 out[x] = 10000;
137 else {
138 double dx = X[i] - Xmean;
139 double dy = Y[i] - Ymean;
140 out[x] = Background +
141 coefficientNorm * Intensity *
142 exp(exponentialCoeffx2 * dx * dx + exponentialCoeffxy * dx * dy + exponentialCoeffy2 * dy * dy);
143 out[x] = out[x] + DDD;
144
145 if (out[x] != out[x]) {
146 out[x] = 100000;
147 isNaNs = true;
148 }
149 }
150 double diff = out[x] - D[x];
151 chiSq += diff * diff;
152
153 x++;
154 }
155 inf << "Constr:";
156 for (size_t i = 0; i < nParams(); i++) {
157 IConstraint *constr = getConstraint(i);
158 if (constr)
159 inf << i << "=" << constr->check() << ";";
160 }
161 inf << "\n\n chiSq =" << chiSq << " nData " << nData << '\n';
162 g_log.debug(inf.str());
163}
164
165void BivariateNormal::functionDeriv1D(API::Jacobian *out, const double *xValues, const size_t nData) {
166 UNUSED_ARG(xValues);
167 UNUSED_ARG(nData);
168 if (nData == 0)
169 return;
170 double penDeriv = initCommon();
171
172 std::ostringstream inf;
173 inf << "***penalty(" << penDeriv << "),Parameters=";
174 for (size_t k = 0; k < 7; k++)
175 inf << "," << LastParams[k];
176 inf << '\n';
177 g_log.debug(inf.str());
178
179 std::vector<double> outf(nData, 0.0);
180 function1D(outf.data(), xValues, nData);
181
182 double paramUU = LastParams[IVXX] * LastParams[IVYY] - LastParams[IVXY] * LastParams[IVXY];
183 /* if( uu <=0)
184 {
185 penDeriv = -2*uu;
186 uu=1;
187 }
188 */
189 API::MatrixWorkspace_const_sptr ws = getMatrixWorkspace();
190 const auto &X = ws->y(1);
191 const auto &Y = ws->y(2);
192
193 for (int x = 0; x < NCells; x++) {
194
195 double penaltyDeriv = penDeriv;
196
197 double r = Y[x];
198 double c = X[x];
199
200 out->set(x, IBACK, +1.0);
201
202 if (penaltyDeriv <= 0)
203 out->set(x, ITINTENS, expVals[x] * coefNorm);
204 else if (LastParams[ITINTENS] < 0)
205 out->set(x, ITINTENS, -.01);
206 else
207 out->set(x, ITINTENS, 0.01);
208
209 double coefExp = coefNorm * LastParams[ITINTENS];
210
211 double coefxy = LastParams[IVXY] / paramUU;
212 double coefx2 = -LastParams[IVYY] / 2 / paramUU;
213
214 if (penaltyDeriv <= 0)
215 out->set(x, IXMEAN,
216 penaltyDeriv +
217 coefExp * expVals[x] * (-2 * coefx2 * (c - LastParams[IXMEAN]) - coefxy * (r - LastParams[IYMEAN])));
218 else // if(LastParams[IXMEAN] < 0)
219 out->set(x, IXMEAN, 0);
220 // else
221 // out->set(x,IXMEAN,0);
222
223 coefExp = coefNorm * LastParams[ITINTENS];
224
225 double coefy2 = -LastParams[IVXX] / 2 / paramUU;
226
227 if (penaltyDeriv <= 0)
228 out->set(x, IYMEAN,
229 penaltyDeriv +
230 coefExp * expVals[x] * (-coefxy * (c - LastParams[IXMEAN]) - 2 * coefy2 * (r - LastParams[IYMEAN])));
231 else // if(LastParams[IYMEAN] < 0)
232 out->set(x, IYMEAN, 0);
233 // else
234 // out->set(x,IYMEAN,0);
235
236 double M = 1;
237 if (nParams() < 5)
238 M = 10;
239 coefExp = coefNorm * LastParams[ITINTENS];
240
241 double C = -LastParams[IVYY] / 2 / paramUU;
242
243 double SIVXX;
244 if (penaltyDeriv <= 0)
245 SIVXX = coefExp * expVals[x] *
246 (C +
247 -LastParams[IVYY] / paramUU *
248 (coefx2 * (c - LastParams[IXMEAN]) * (c - LastParams[IXMEAN]) +
249 coefxy * (r - LastParams[IYMEAN]) * (c - LastParams[IXMEAN]) +
250 coefy2 * (r - LastParams[IYMEAN]) * (r - LastParams[IYMEAN])) -
251 (r - LastParams[IYMEAN]) * (r - LastParams[IYMEAN]) / 2 / paramUU);
252 else if (LastParams[IVXX] < .01)
253 SIVXX = -M;
254 else
255 SIVXX = 0;
256
257 if (LastParams[IVXX] > 1.2 * Varx0 && !CalcVariances)
258 SIVXX = 0;
259 if (LastParams[IVXX] < .8 * Varx0)
260 SIVXX = 0;
261 coefExp = coefNorm * LastParams[ITINTENS];
262
263 C = -LastParams[IVXX] / 2 / paramUU;
264
265 double SIVYY;
266 if (penaltyDeriv <= 0)
267 SIVYY = coefExp * expVals[x] *
268 (C +
269 -LastParams[IVXX] / paramUU *
270 (coefx2 * (c - LastParams[IXMEAN]) * (c - LastParams[IXMEAN]) +
271 coefxy * (r - LastParams[IYMEAN]) * (c - LastParams[IXMEAN]) +
272 coefy2 * (r - LastParams[IYMEAN]) * (r - LastParams[IYMEAN])) -
273 (c - LastParams[IXMEAN]) * (c - LastParams[IXMEAN]) / 2 / paramUU);
274
275 else if (LastParams[IVYY] < .01)
276 SIVYY = -M;
277 else
278 SIVYY = 0;
279
280 if (LastParams[IVYY] > 1.2 * Vary0 && !CalcVariances)
281 SIVYY = 0;
282
283 if (LastParams[IVYY] < .8 * Vary0)
284 SIVYY = 0;
285 coefExp = coefNorm * LastParams[ITINTENS];
286
287 C = LastParams[IVXY] / paramUU;
288
289 double SIVXY;
290 if (penaltyDeriv <= 0)
291 SIVXY = coefExp * expVals[x] *
292 (C +
293 2 * LastParams[IVXY] / paramUU *
294 (coefx2 * (c - LastParams[IXMEAN]) * (c - LastParams[IXMEAN]) +
295 coefxy * (r - LastParams[IYMEAN]) * (c - LastParams[IXMEAN]) +
296 coefy2 * (r - LastParams[IYMEAN]) * (r - LastParams[IYMEAN])) +
297 (r - LastParams[IYMEAN]) * (c - LastParams[IXMEAN]) / paramUU);
298
299 else if (paramUU < 0)
300 SIVXY = 2 * LastParams[IVXY];
301 else
302 SIVXY = 0;
303
304 if (!CalcVxx && nParams() > 6)
305 out->set(x, IVXX, penaltyDeriv + SIVXX);
306 else {
307 // out->set(x,IVXX,0.0);
308 double bdderiv = out->get(x, IBACK);
309
310 bdderiv += SIVXX *
311 (-Sxx - (LastParams[IXMEAN] - mx) * (LastParams[IXMEAN] - mx) * TotN + LastParams[IVXX] * TotN) /
312 (TotI - LastParams[IBACK] * TotN);
313
314 out->set(x, IBACK, bdderiv);
315
316 double mxderiv = out->get(x, IXMEAN);
317 mxderiv += SIVXX *
318 (2 * (LastParams[IXMEAN] - mIx) *
319
320 TotI -
321 2 * LastParams[IBACK] * (LastParams[IXMEAN] - mx) * TotN) /
322 (TotI - LastParams[IBACK] * TotN);
323 out->set(x, IXMEAN, mxderiv);
324 }
325 if (!CalcVyy && nParams() > 6)
326 out->set(x, IVYY, penaltyDeriv + SIVYY);
327 else {
328 // out->set(x,IVYY, 0.0);
329 double bdderiv = out->get(x, IBACK);
330
331 bdderiv += SIVYY *
332 (-Syy - (LastParams[IYMEAN] - my) * (LastParams[IYMEAN] - my) * TotN + LastParams[IVYY] * TotN) /
333 (TotI - LastParams[IBACK] * TotN);
334 out->set(x, IBACK, bdderiv);
335 double myderiv = out->get(x, IYMEAN);
336
337 myderiv += SIVYY *
338 (2 * (LastParams[IYMEAN] - mIy) *
339
340 TotI -
341 2 * LastParams[IBACK] * (LastParams[IYMEAN] - my) * TotN) /
342 (TotI - LastParams[IBACK] * TotN);
343 out->set(x, IYMEAN, myderiv);
344 }
345 if (!CalcVxy && nParams() > 6) {
346 out->set(x, IVXY, penaltyDeriv + SIVXY);
347
348 } else {
349 // out->set(x,IVXY, 0.0);
350 double bdderiv = out->get(x, IBACK);
351 bdderiv += SIVXY *
352 (-Sxy - (LastParams[IYMEAN] - my) * (LastParams[IXMEAN] - mx) * TotN + LastParams[IVXY] * TotN) /
353 (TotI - LastParams[IBACK] * TotN);
354 out->set(x, IBACK, bdderiv);
355 double myderiv = out->get(x, IYMEAN);
356 myderiv += SIVXY *
357 ((LastParams[IXMEAN] - mIx) *
358
359 TotI -
360 LastParams[IBACK] * (LastParams[IXMEAN] - mx) * TotN) /
361 (TotI - LastParams[IBACK] * TotN);
362 out->set(x, IYMEAN, myderiv);
363 double mxderiv = out->get(x, IXMEAN);
364 mxderiv += SIVXY * ((LastParams[IYMEAN] - mIy) * TotI - LastParams[IBACK] * (LastParams[IYMEAN] - my) * TotN) /
365 (TotI - LastParams[IBACK] * TotN);
366 out->set(x, IXMEAN, mxderiv);
367 }
368 }
369}
370
371void BivariateNormal::init() {
372 declareParameter("Background", 0.00);
373 declareParameter("Intensity", 0.00);
374 declareParameter("Mcol", 0.00, "Mean column(x) value");
375 declareParameter("Mrow", 0.00, "Mean row(y) value");
376
377 CalcVariances = false;
378
379 NCells = -1;
380 LastParams[IVXX] = -1;
381}
382
383double BivariateNormal::initCommon() {
384
385 double penalty = 0;
386 bool ParamsOK = true;
387 bool CommonsOK = true;
388 if (!expVals)
389 CommonsOK = false;
390
391 API::MatrixWorkspace_const_sptr ws = getMatrixWorkspace();
392 const auto &D = ws->y(0);
393 const auto &X = ws->y(1);
394 const auto &Y = ws->y(2);
395
396 if (NCells < 0) {
397 NCells = static_cast<int>(std::min<size_t>(D.size(), std::min<size_t>(X.size(), Y.size())));
398 CommonsOK = false;
399 }
400
401 double Attrib[12] = {0.0};
402
403 double MinX, MinY, MaxX, MaxY, MaxD, MinD;
404 MinX = MaxX = X[0];
405 MinY = MaxY = Y[0];
406 MaxD = MinD = D[0];
407
408 if (!CommonsOK) {
409
410 for (int i = 0; i < NCells; i++) {
411 Attrib[S_int] += D[i];
412 Attrib[S_xint] += D[i] * X[i];
413 Attrib[S_yint] += D[i] * Y[i];
414 Attrib[S_x2int] += D[i] * X[i] * X[i];
415 Attrib[S_y2int] += D[i] * Y[i] * Y[i];
416 Attrib[S_xyint] += D[i] * X[i] * Y[i];
417
418 Attrib[S_y] += Y[i];
419 Attrib[S_x] += X[i];
420 Attrib[S_x2] += X[i] * X[i];
421 Attrib[S_y2] += Y[i] * Y[i];
422 Attrib[S_xy] += X[i] * Y[i];
423 Attrib[S_1] += 1.0;
424
425 if (X[i] < MinX)
426 MinX = X[i];
427 if (X[i] > MaxX)
428 MaxX = X[i];
429 if (Y[i] < MinY)
430 MinY = Y[i];
431 if (Y[i] > MaxY)
432 MaxY = Y[i];
433 if (D[i] < MinD)
434 MinD = D[i];
435 if (D[i] > MaxD)
436 MaxD = D[i];
437 }
438
439 mIx = Attrib[S_xint] / Attrib[S_int];
440 mIy = Attrib[S_yint] / Attrib[S_int];
441 mx = Attrib[S_x] / Attrib[S_1];
442 my = Attrib[S_y] / Attrib[S_1];
443
444 SIxx = Attrib[S_x2int] - (Attrib[S_xint] * Attrib[S_xint]) / Attrib[S_int];
445 SIyy = Attrib[S_y2int] - (Attrib[S_yint]) * (Attrib[S_yint]) / Attrib[S_int];
446 SIxy = Attrib[S_xyint] - (Attrib[S_xint]) * (Attrib[S_yint]) / Attrib[S_int];
447
448 Sxx = Attrib[S_x2] - (Attrib[S_x]) * (Attrib[S_x]) / Attrib[S_1];
449 Syy = Attrib[S_y2] - (Attrib[S_y]) * (Attrib[S_y]) / Attrib[S_1];
450 Sxy = Attrib[S_xy] - (Attrib[S_x]) * (Attrib[S_y]) / Attrib[S_1];
451
452 // CommonsOK = false;
453
454 TotI = Attrib[S_int];
455 TotN = Attrib[S_1];
456
457 // CommonsOK = false;
458
459 if (getConstraint(0) == nullptr) {
460
461 addConstraint((std::make_unique<BoundaryConstraint>(this, "Background", 0, Attrib[S_int] / Attrib[S_1])));
462 }
463
464 double maxIntensity = Attrib[S_int] + 3 * sqrt(Attrib[S_int]);
465
466 if (maxIntensity < 100)
467 maxIntensity = 100;
468
469 if (getConstraint(1) == nullptr) {
470 addConstraint(std::make_unique<BoundaryConstraint>(this, "Intensity", 0, maxIntensity));
471 }
472
473 if (getConstraint(3) == nullptr) {
474 double minMeany = MinY * .9 + .1 * MaxY;
475 double maxMeany = MinY * .1 + .9 * MaxY;
476 addConstraint(std::make_unique<BoundaryConstraint>(this, "Mrow", minMeany, maxMeany));
477 }
478
479 if (getConstraint(2) == nullptr) {
480 double minMeanx = MinX * .9 + .1 * MaxX;
481 double maxMeanx = MinX * .1 + .9 * MaxX;
482 addConstraint(std::make_unique<BoundaryConstraint>(this, "Mcol", minMeanx, maxMeanx));
483 }
484
485 if (CalcVariances && nParams() > 6) {
486 std::ostringstream ssxx, ssyy, ssxy;
487
488 ssyy << std::string("(") << (SIyy) << "+(Mrow-" << (mIy) << ")*(Mrow-" << (mIy) << ")*" << Attrib[S_int]
489 << "-Background*" << (Syy) << "-Background*(Mrow-" << (my) << ")*(Mrow-" << (my) << ")*" << Attrib[S_1]
490 << ")/(" << (Attrib[S_int]) << "-Background*" << (Attrib[S_1]) << ")";
491
492 if (getTie(IVYY) == nullptr) {
493 tie("SSrow", ssyy.str());
494 CalcVxx = true;
495 }
496
497 ssxx << std::string("(") << (SIxx) << "+(Mcol-" << (mIx) << ")*(Mcol-" << (mIx) << ")*" << Attrib[S_int]
498 << "-Background*" << (Sxx) << "-Background*(Mcol-" << (mx) << ")*(Mcol-" << (mx) << ")*" << Attrib[S_1]
499 << ")/(" << (Attrib[S_int]) << "-Background*" << (Attrib[S_1]) << ")";
500
501 if (getTie(IVXX) == nullptr) {
502 tie("SScol", ssxx.str());
503 CalcVyy = true;
504 }
505
506 ssxy << std::string("(") << (SIxy) << "+(Mcol-" << (mIx) << ")*(Mrow-" << (mIy) << ")*" << Attrib[S_int]
507 << "-Background*" << (Sxy) << "-Background*(Mcol-" << (mx) << ")*(Mrow-" << (my) << ")*" << Attrib[S_1]
508 << ")/(" << (Attrib[S_int]) << "-Background*" << (Attrib[S_1]) << ")";
509
510 if (getTie(IVXY) == nullptr) {
511 tie("SSrc", ssxy.str());
512 CalcVxy = true;
513 }
514 }
515 CommonsOK = true;
516 }
517
518 if (LastParams[IVXX] < 0) {
519 ParamsOK = false;
520 CommonsOK = false;
521
522 } else
523 for (size_t i = 0; i < nParams() && ParamsOK; i++)
524 if (getParameter(i) != LastParams[i])
525 ParamsOK = false;
526
527 if (!ParamsOK) {
528
529 for (size_t i = 0; i < nParams(); i++)
530 LastParams[i] = getParameter(i);
531 }
532
533 if (!CommonsOK || !ParamsOK) {
534
535 int NCells1;
536 double Varxx, Varxy, Varyy;
537
538 Varxx = Varxy = Varyy = -1;
539 penalty = initCoeff(D, X, Y, coefNorm, expCoeffx2, expCoeffy2, expCoeffxy, NCells1, Varxx, Varxy, Varyy);
540
541 if (Varx0 < 0 && penalty <= 0) {
542 Varx0 = Varxx;
543 Vary0 = Varyy;
544 }
545
546 LastParams[IVXX] = Varxx;
547 LastParams[IVXY] = Varxy;
548 LastParams[IVYY] = Varyy;
549
550 delete[] expVals;
551 expVals = new double[NCells];
552
553 for (int i = 0; i < NCells; i++) {
554
555 double dx = X[i] - LastParams[IXMEAN];
556 double dy = Y[i] - LastParams[IYMEAN];
557 expVals[i] = exp(expCoeffx2 * dx * dx + expCoeffxy * dx * dy + expCoeffy2 * dy * dy);
558 }
559 }
560 return penalty;
561}
562
563double BivariateNormal::initCoeff(const HistogramY &D, const HistogramY &X, const HistogramY &Y, double &coefNorm,
564 double &expCoeffx2, double &expCoeffy2, double &expCoeffxy, int &NCells,
565 double &Varxx, double &Varxy, double &Varyy) const {
566
567 double Background = getParameter("Background");
568 bool zeroDenom = false;
569 if (TotI == 0 && TotN == 0)
570 zeroDenom = true;
571 else if (TotI - Background * TotN <= 0)
572 zeroDenom = true;
573 if (CalcVxx || nParams() < 6) {
574 Varxx =
575 (SIxx + (getParameter("Mcol") - mIx) * (getParameter("Mcol") - mIx) * TotI - getParameter("Background") * Sxx -
576 getParameter("Background") * (getParameter("Mcol") - (mx)) * (getParameter("Mcol") - (mx)) * TotN) /
577 (TotI - getParameter("Background") * TotN);
578
579 if (Varx0 > 0) {
580 Varxx = std::min<double>(Varxx, 1.21 * Varx0);
581 Varxx = std::max<double>(Varxx, .79 * Varx0);
582 }
583
584 } else {
585 Varxx = getParameter(IVXX);
586 }
587
588 double Mrow = getParameter("Mrow");
589
590 if (CalcVyy || nParams() < 6) {
591
592 Varyy = (SIyy + (Mrow - (mIy)) * (Mrow - (mIy)) * TotI -
593 getParameter("Background") * (Syy)-Background * (Mrow - (my)) * (Mrow - (my)) * TotN) /
594 (TotI - Background * TotN);
595 if (Vary0 > 0) {
596 Varyy = std::min<double>(Varyy, 1.21 * Vary0);
597 Varyy = std::max<double>(Varyy, .79 * Vary0);
598 }
599 } else {
600 Varyy = getParameter(IVYY);
601 }
602
603 if (CalcVxy || nParams() < 6) {
604 Varxy = ((SIxy) + (getParameter("Mcol") - (mIx)) * (getParameter("Mrow") - (mIy)) * TotI -
605 getParameter("Background") * (Sxy)-getParameter("Background") * (getParameter("Mcol") - (mx)) *
606 (getParameter("Mrow") - (my)) * TotN) /
607 (TotI - getParameter("Background") * TotN);
608
609 } else {
610 Varxy = getParameter(IVXY);
611 }
612
613 double paramUU = Varxx * Varyy - Varxy * Varxy;
614 double penalty;
615 if (zeroDenom)
616 penalty = 200;
617 else
618 penalty =
619 std::max<double>(0, -Varxx + .01) + std::max<double>(0, -Varyy + .01) + std::max<double>(0, -paramUU + .01);
620
621 if (fabs(paramUU) < .01) {
622 if (paramUU < 0)
623 paramUU = -.01;
624 else
625 paramUU = .01;
626 }
627
628 NCells = static_cast<int>(std::min<size_t>(D.size(), std::min<size_t>(X.size(), Y.size())));
629 if (zeroDenom) {
630 coefNorm = expCoeffx2 = expCoeffy2 = 1;
631 expCoeffxy = 0;
632 Varxx = Varyy = 5;
633 Varxy = 0;
634 return penalty;
635 }
636 coefNorm = .5 / M_PI / sqrt(fabs(paramUU));
637
638 expCoeffx2 = -fabs(Varyy) / 2 / fabs(paramUU);
639 expCoeffxy = Varxy / paramUU;
640 expCoeffy2 = -fabs(Varxx) / 2 / fabs(paramUU);
641
642 if (nParams() < 5)
643 penalty *= 10;
644
645 return penalty;
646}
647
648} // 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:48
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:233
Mantid::API::IFunction::getTie
virtual ParameterTie * getTie(size_t i) const
Get the tie of i-th parameter.
Definition IFunction.cpp:432
Mantid::API::IFunction::addConstraint
virtual void addConstraint(std::unique_ptr< IConstraint > ic)
Add a constraint to function.
Definition IFunction.cpp:452
Mantid::API::IFunction::getConstraint
virtual IConstraint * getConstraint(size_t i) const
Get constraint of i-th parameter.
Definition IFunction.cpp:468
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:165
Mantid::CurveFitting::Functions::BivariateNormal::initCommon
double initCommon()
Check for changes in parameters, etc.
Definition BivariateNormal.cpp:383
Mantid::CurveFitting::Functions::BivariateNormal::init
void init() override
Function initialization. Declare function parameters in this method.
Definition BivariateNormal.cpp:371
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:563
Mantid::Kernel::Logger
The Logger class is in charge of the publishing messages from the framework through various channels.
Definition Logger.h:51
Mantid::Kernel::Logger::debug
void debug(const std::string &msg)
Logs at debug level.
Definition Logger.cpp:145
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:587
Generated by doxygen 1.9.8