File indexing completed on 2023-10-25 10:05:23
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037 #include "TopQuarkAnalysis/TopHitFit/interface/Chisq_Constrainer.h"
0038 #include "TopQuarkAnalysis/TopHitFit/interface/Defaults.h"
0039 #include <cmath>
0040 #include <cassert>
0041 #include <iostream>
0042 #include <iomanip>
0043
0044 using std::abs;
0045 using std::cout;
0046 using std::fixed;
0047 using std::ios_base;
0048 using std::ostream;
0049 using std::resetiosflags;
0050 using std::setiosflags;
0051 using std::sqrt;
0052
0053
0054
0055 namespace hitfit {
0056
0057 Chisq_Constrainer_Args::Chisq_Constrainer_Args(const Defaults& defs)
0058
0059
0060
0061
0062
0063
0064 : _base_constrainer_args(defs) {
0065 _use_G = defs.get_bool("use_G");
0066 _printfit = defs.get_bool("printfit");
0067 _constraint_sum_eps = defs.get_float("constraint_sum_eps");
0068 _chisq_diff_eps = defs.get_float("chisq_diff_eps");
0069 _maxit = defs.get_int("maxit");
0070 _max_cut = defs.get_int("maxcut");
0071 _cutsize = defs.get_float("cutsize");
0072 _min_tot_cutsize = defs.get_float("min_tot_cutsize");
0073 _chisq_test_eps = defs.get_float("chisq_test_eps");
0074 }
0075
0076 bool Chisq_Constrainer_Args::printfit() const
0077
0078
0079
0080
0081 {
0082 return _printfit;
0083 }
0084
0085 bool Chisq_Constrainer_Args::use_G() const
0086
0087
0088
0089
0090 {
0091 return _use_G;
0092 }
0093
0094 double Chisq_Constrainer_Args::constraint_sum_eps() const
0095
0096
0097
0098
0099 {
0100 return _constraint_sum_eps;
0101 }
0102
0103 double Chisq_Constrainer_Args::chisq_diff_eps() const
0104
0105
0106
0107
0108 {
0109 return _chisq_diff_eps;
0110 }
0111
0112 unsigned Chisq_Constrainer_Args::maxit() const
0113
0114
0115
0116
0117 {
0118 return _maxit;
0119 }
0120
0121 unsigned Chisq_Constrainer_Args::max_cut() const
0122
0123
0124
0125
0126 {
0127 return _max_cut;
0128 }
0129
0130 double Chisq_Constrainer_Args::cutsize() const
0131
0132
0133
0134
0135 {
0136 return _cutsize;
0137 }
0138
0139 double Chisq_Constrainer_Args::min_tot_cutsize() const
0140
0141
0142
0143
0144 {
0145 return _min_tot_cutsize;
0146 }
0147
0148 double Chisq_Constrainer_Args::chisq_test_eps() const
0149
0150
0151
0152
0153 {
0154 return _chisq_test_eps;
0155 }
0156
0157 const Base_Constrainer_Args& Chisq_Constrainer_Args::base_constrainer_args() const
0158
0159
0160
0161
0162 {
0163 return _base_constrainer_args;
0164 }
0165
0166 }
0167
0168
0169
0170 namespace {
0171
0172 using namespace hitfit;
0173
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189
0190 bool solve_linear_system(const Matrix& H,
0191 const Diagonal_Matrix& Y,
0192 const Matrix& By,
0193 const Row_Vector& r,
0194 Column_Vector& alpha,
0195 Column_Vector& d,
0196 Matrix& W,
0197 Matrix& U,
0198 Matrix& V)
0199
0200
0201
0202
0203
0204
0205
0206
0207
0208
0209
0210
0211
0212
0213
0214
0215
0216 {
0217 int nconstraints = H.num_row();
0218 int nbadvars = Y.num_row();
0219
0220
0221 Matrix A(nconstraints + nbadvars, nconstraints + nbadvars);
0222 A.sub(1, 1, -H);
0223 if (nbadvars > 0) {
0224 A.sub(nconstraints + 1, nconstraints + 1, Y);
0225 A.sub(1, nconstraints + 1, By.T());
0226 A.sub(nconstraints + 1, 1, By);
0227 }
0228
0229
0230 Column_Vector yy(nconstraints + nbadvars, 0);
0231 yy.sub(1, r.T());
0232
0233
0234
0235 Matrix Ai;
0236 int ierr = 0;
0237 do {
0238 Ai = A.inverse(ierr);
0239 if (ierr) {
0240 int allzero = 0;
0241 for (int i = 1; i <= nconstraints; i++) {
0242 allzero = 1;
0243 for (int j = 1; j <= nconstraints; j++) {
0244 if (A(i, j) != 0) {
0245 allzero = 0;
0246 break;
0247 }
0248 }
0249 if (allzero) {
0250 A(i, i) = 1;
0251 break;
0252 }
0253 }
0254 if (!allzero)
0255 return false;
0256 }
0257 } while (ierr != 0);
0258
0259
0260 Column_Vector xx = Ai * yy;
0261
0262
0263
0264 W = Ai.sub(1, nconstraints, 1, nconstraints);
0265 if (nbadvars > 0) {
0266 U = Ai.sub(nconstraints + 1, nconstraints + nbadvars, nconstraints + 1, nconstraints + nbadvars);
0267 V = Ai.sub(nconstraints + 1, nconstraints + nbadvars, 1, nconstraints);
0268 d = xx.sub(nconstraints + 1, nconstraints + nbadvars);
0269 }
0270
0271 alpha = xx.sub(1, nconstraints);
0272
0273 return true;
0274 }
0275
0276 }
0277
0278 namespace hitfit {
0279
0280
0281
0282 Chisq_Constrainer::Chisq_Constrainer(const Chisq_Constrainer_Args& args)
0283
0284
0285
0286
0287
0288
0289 : Base_Constrainer(args.base_constrainer_args()), _args(args) {}
0290
0291 double Chisq_Constrainer::fit(Constraint_Calculator& constraint_calculator,
0292 const Column_Vector& xm,
0293 Column_Vector& x,
0294 const Column_Vector& ym,
0295 Column_Vector& y,
0296 const Matrix& G_i,
0297 const Diagonal_Matrix& Y,
0298 Column_Vector& pullx,
0299 Column_Vector& pully,
0300 Matrix& Q,
0301 Matrix& R,
0302 Matrix& S)
0303
0304
0305
0306
0307
0308
0309
0310
0311
0312
0313
0314
0315
0316
0317
0318
0319
0320
0321
0322
0323
0324
0325
0326
0327
0328
0329
0330
0331
0332 {
0333
0334
0335 int nvars = x.num_row();
0336 assert(nvars == G_i.num_col());
0337 assert(nvars == xm.num_row());
0338
0339 int nbadvars = y.num_row();
0340 assert(nbadvars == Y.num_col());
0341 assert(nbadvars == ym.num_row());
0342
0343
0344
0345 Matrix G(nvars, nvars);
0346 if (_args.use_G()) {
0347 int ierr = 0;
0348 G = G_i.inverse(ierr);
0349 assert(!ierr);
0350 }
0351
0352 int nconstraints = constraint_calculator.nconstraints();
0353
0354
0355 Row_Vector F(nconstraints);
0356 Matrix Bx(nvars, nconstraints);
0357 Matrix By(nbadvars, nconstraints);
0358
0359
0360
0361
0362 if (!call_constraint_fcn(constraint_calculator, x, y, F, Bx, By)) {
0363
0364
0365 return -999.0;
0366 }
0367
0368
0369 double constraint_sum_last = -1000;
0370 double chisq_last = -1000;
0371 bool near_convergence = false;
0372 double last_step_cutsize = 1;
0373
0374 unsigned nit = 0;
0375
0376
0377 Column_Vector c = x - xm;
0378 Column_Vector d = y - ym;
0379
0380 Matrix E(nvars, nconstraints);
0381 Matrix W(nconstraints, nconstraints);
0382 Matrix U(nbadvars, nbadvars);
0383 Matrix V(nbadvars, nconstraints);
0384
0385
0386 do {
0387
0388 E = G_i * Bx;
0389 Matrix H = E.T() * Bx;
0390 Row_Vector r = c.T() * Bx + d.T() * By - F;
0391
0392
0393
0394
0395 Column_Vector alpha(nvars);
0396 Column_Vector d1(nbadvars);
0397 if (!solve_linear_system(H, Y, By, r, alpha, d1, W, U, V)) {
0398
0399
0400 return -998.0;
0401 }
0402
0403
0404 Column_Vector c1 = -E * alpha;
0405 double chisq = -scalar(r * alpha);
0406
0407 double psi_cut = 0;
0408
0409
0410 x = c1 + xm;
0411 y = d1 + ym;
0412
0413
0414 Matrix save_By = By;
0415 Row_Vector save_negF = -F;
0416 double this_step_cutsize = 1;
0417 double constraint_sum = -1;
0418 unsigned ncut = 0;
0419
0420
0421
0422
0423
0424
0425
0426
0427
0428
0429 while (!call_constraint_fcn(constraint_calculator, x, y, F, Bx, By) ||
0430 ((constraint_sum = norm_infinity(F)) > _args.constraint_sum_eps() && nit > 0 &&
0431 constraint_sum > constraint_sum_last)) {
0432
0433 if (nit > 0 && _args.printfit() && ncut == 0) {
0434 cout << "(" << chisq << " " << chisq_last << ") ";
0435 }
0436
0437
0438
0439
0440 if (ncut == 0 && abs(chisq - chisq_last) < _args.chisq_diff_eps()) {
0441
0442
0443
0444 if (_args.printfit())
0445 cout << " directed step ";
0446
0447
0448
0449 Column_Vector beta(nconstraints);
0450 Column_Vector delta(nbadvars);
0451 solve_linear_system(H, Y, save_By, save_negF, beta, delta, W, U, V);
0452
0453
0454 Column_Vector gamma = -E * beta;
0455
0456
0457 x = c + xm + gamma;
0458 y = d + ym + delta;
0459
0460
0461
0462 if (call_constraint_fcn(constraint_calculator, x, y, F, Bx, By) && (constraint_sum = norm_infinity(F)) > 0 &&
0463 (constraint_sum < constraint_sum_last)) {
0464
0465
0466 chisq = chisq_last - scalar((-save_negF + r * 2) * beta);
0467 c1 = x - xm;
0468 d1 = y - ym;
0469
0470
0471 break;
0472 }
0473 }
0474
0475
0476
0477 if (ncut == 0)
0478 psi_cut = scalar((save_negF - r) * alpha);
0479
0480
0481 if (++ncut > _args.max_cut()) {
0482
0483
0484 return -997.0;
0485 }
0486
0487
0488
0489
0490
0491
0492 double this_cutsize = _args.cutsize();
0493 if (ncut == 1 && last_step_cutsize < 1) {
0494 this_cutsize = 2 * last_step_cutsize;
0495 if (this_cutsize > _args.cutsize())
0496 this_cutsize = _args.cutsize();
0497 }
0498
0499
0500 this_step_cutsize *= this_cutsize;
0501
0502
0503 if (this_step_cutsize < _args.min_tot_cutsize()) {
0504
0505
0506 return -996.0;
0507 }
0508
0509
0510 double cutleft = 1 - this_cutsize;
0511 c1 = c1 * this_cutsize + c * cutleft;
0512 d1 = d1 * this_cutsize + d * cutleft;
0513
0514
0515 if (chisq_last >= 0) {
0516 chisq = this_cutsize * this_cutsize * chisq + cutleft * cutleft * chisq_last +
0517 2 * this_cutsize * cutleft * psi_cut;
0518 psi_cut = this_cutsize * psi_cut + cutleft * chisq_last;
0519 } else
0520 chisq = chisq_last;
0521
0522
0523 if (_args.printfit()) {
0524 cout << constraint_sum << " cut " << ncut << " size " << setiosflags(ios_base::scientific) << this_cutsize
0525 << " tot size " << this_step_cutsize << resetiosflags(ios_base::scientific) << " " << chisq << "\n";
0526 }
0527
0528
0529 x = c1 + xm;
0530 y = d1 + ym;
0531
0532
0533 }
0534
0535
0536
0537 last_step_cutsize = this_step_cutsize;
0538
0539
0540
0541 double chisq_b = 0;
0542 if (_args.use_G()) {
0543 chisq_b = scalar(c1.T() * G * c1) + scalar(d1.T() * Y * d1);
0544 if (chisq >= 0 && abs((chisq - chisq_b) / chisq) > _args.chisq_test_eps()) {
0545 cout << chisq << " " << chisq_b << "lost precision?\n";
0546 abort();
0547 }
0548 }
0549
0550
0551 if (_args.printfit()) {
0552 cout << chisq << " ";
0553 if (_args.use_G())
0554 cout << chisq_b << " ";
0555 }
0556
0557 double z2 = abs(chisq - chisq_last);
0558
0559 if (_args.printfit()) {
0560 cout << constraint_sum << " " << z2 << "\n";
0561 }
0562
0563 c = c1;
0564 d = d1;
0565 chisq_last = chisq;
0566 constraint_sum_last = constraint_sum;
0567
0568
0569
0570 if (chisq >= 0 && constraint_sum < _args.constraint_sum_eps() && z2 < _args.chisq_diff_eps()) {
0571 if (near_convergence)
0572 break;
0573 near_convergence = true;
0574 } else
0575 near_convergence = false;
0576
0577
0578 if (++nit > _args.maxit()) {
0579
0580
0581 return -995.0;
0582 }
0583
0584 } while (true);
0585
0586
0587
0588
0589 Q = E * W * E.T();
0590 S = -E * V.T();
0591 R = U;
0592
0593
0594 pullx = Column_Vector(nvars);
0595 for (int i = 1; i <= nvars; i++) {
0596 double a = Q(i, i);
0597 if (a < 0)
0598 pullx(i) = c(i) / sqrt(-a);
0599 else {
0600 pullx(i) = 0;
0601
0602 }
0603 }
0604
0605 pully = Column_Vector(nbadvars);
0606 for (int i = 1; i <= nbadvars; i++) {
0607 double a = 1 - Y(i, i) * R(i, i);
0608 if (a > 0)
0609 pully(i) = d(i) * sqrt(Y(i, i) / a);
0610 else {
0611 pully(i) = 0;
0612
0613 }
0614 }
0615
0616
0617 Q = Q + G_i;
0618
0619
0620 return chisq_last;
0621 }
0622
0623 std::ostream& Chisq_Constrainer::print(std::ostream& s) const
0624
0625
0626
0627
0628
0629
0630
0631
0632
0633 {
0634 Base_Constrainer::print(s);
0635 s << " printfit: " << _args.printfit() << " use_G: " << _args.use_G() << "\n";
0636 s << " constraint_sum_eps: " << _args.constraint_sum_eps() << " chisq_diff_eps: " << _args.chisq_diff_eps()
0637 << " chisq_test_eps: " << _args.chisq_test_eps() << "\n";
0638 s << " maxit: " << _args.maxit() << " max_cut: " << _args.max_cut()
0639 << " min_tot_cutsize: " << _args.min_tot_cutsize() << " cutsize: " << _args.cutsize() << "\n";
0640 return s;
0641 }
0642
0643 }