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