CMS 3D CMS Logo

Chisq_Constrainer.cc
Go to the documentation of this file.
1 //
2 //
3 // File: src/Chisq_Constrainer.cc
4 // Purpose: Minimize a chisq subject to a set of constraints.
5 // Based on the SQUAW algorithm.
6 // Created: Jul, 2000, sss, based on run 1 mass analysis code.
7 //
8 // CMSSW File : src/Chisq_Constrainer.cc
9 // Original Author : Scott Stuart Snyder <snyder@bnl.gov> for D0
10 // Imported to CMSSW by Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>
11 //
12 
39 #include <cmath>
40 #include <cassert>
41 #include <iostream>
42 #include <iomanip>
43 
44 using std::abs;
45 using std::cout;
46 using std::fixed;
47 using std::ios_base;
48 using std::ostream;
49 using std::resetiosflags;
50 using std::setiosflags;
51 using std::sqrt;
52 
53 //*************************************************************************
54 
55 namespace hitfit {
56 
58  //
59  // Purpose: Constructor.
60  //
61  // Inputs:
62  // defs - The Defaults instance from which to initialize.
63  //
64  : _base_constrainer_args(defs) {
65  _use_G = defs.get_bool("use_G");
66  _printfit = defs.get_bool("printfit");
67  _constraint_sum_eps = defs.get_float("constraint_sum_eps");
68  _chisq_diff_eps = defs.get_float("chisq_diff_eps");
69  _maxit = defs.get_int("maxit");
70  _max_cut = defs.get_int("maxcut");
71  _cutsize = defs.get_float("cutsize");
72  _min_tot_cutsize = defs.get_float("min_tot_cutsize");
73  _chisq_test_eps = defs.get_float("chisq_test_eps");
74  }
75 
77  //
78  // Purpose: Return the printfit parameter.
79  // See the header for documentation.
80  //
81  {
82  return _printfit;
83  }
84 
86  //
87  // Purpose: Return the use_G parameter.
88  // See the header for documentation.
89  //
90  {
91  return _use_G;
92  }
93 
95  //
96  // Purpose: Return the constraint_sum_eps parameter.
97  // See the header for documentation.
98  //
99  {
100  return _constraint_sum_eps;
101  }
102 
104  //
105  // Purpose: Return the chisq_diff_eps parameter.
106  // See the header for documentation.
107  //
108  {
109  return _chisq_diff_eps;
110  }
111 
113  //
114  // Purpose: Return the maxit parameter.
115  // See the header for documentation.
116  //
117  {
118  return _maxit;
119  }
120 
122  //
123  // Purpose: Return the max_cut parameter.
124  // See the header for documentation.
125  //
126  {
127  return _max_cut;
128  }
129 
131  //
132  // Purpose: Return the cutsize parameter.
133  // See the header for documentation.
134  //
135  {
136  return _cutsize;
137  }
138 
140  //
141  // Purpose: Return the min_tot_cutsize parameter.
142  // See the header for documentation.
143  //
144  {
145  return _min_tot_cutsize;
146  }
147 
149  //
150  // Purpose: Return the chisq_test_eps parameter.
151  // See the header for documentation.
152  //
153  {
154  return _chisq_test_eps;
155  }
156 
158  //
159  // Purpose: Return the contained subobject parameters.
160  // See the header for documentation.
161  //
162  {
163  return _base_constrainer_args;
164  }
165 
166 } // namespace hitfit
167 
168 //*************************************************************************
169 
170 namespace {
171 
172  using namespace hitfit;
173 
190  bool solve_linear_system(const Matrix& H,
191  const Diagonal_Matrix& Y,
192  const Matrix& By,
193  const Row_Vector& r,
195  Column_Vector& d,
196  Matrix& W,
197  Matrix& U,
198  Matrix& V)
199  //
200  // Purpose: Solve the system
201  //
202  // [ -H B.t ] [ alpha ] [ r ]
203  // [ ] * [ ] = [ ]
204  // [ B Y ] [ d ] [ 0 ]
205  //
206  // for alpha and d.
207  //
208  // Also returns the inverse matrices:
209  //
210  // [ W V.t ] [ -H B.t ]
211  // [ ] = [ ] ^ -1
212  // [ V U ] [ B Y ]
213  //
214  // Returns true if successful, false if not.
215  //
216  {
217  int nconstraints = H.num_row();
218  int nbadvars = Y.num_row();
219 
220  // Form the matrix on the LHS from H, By, and Y.
221  Matrix A(nconstraints + nbadvars, nconstraints + nbadvars);
222  A.sub(1, 1, -H);
223  if (nbadvars > 0) {
224  A.sub(nconstraints + 1, nconstraints + 1, Y);
225  A.sub(1, nconstraints + 1, By.T());
226  A.sub(nconstraints + 1, 1, By);
227  }
228 
229  // Form the RHS vector from r.
230  Column_Vector yy(nconstraints + nbadvars, 0);
231  yy.sub(1, r.T());
232 
233  // Invert the matrix.
234  // Try to handle singularities correctly.
235  Matrix Ai;
236  int ierr = 0;
237  do {
238  Ai = A.inverse(ierr);
239  if (ierr) {
240  int allzero = 0;
241  for (int i = 1; i <= nconstraints; i++) {
242  allzero = 1;
243  for (int j = 1; j <= nconstraints; j++) {
244  if (A(i, j) != 0) {
245  allzero = 0;
246  break;
247  }
248  }
249  if (allzero) {
250  A(i, i) = 1;
251  break;
252  }
253  }
254  if (!allzero)
255  return false;
256  }
257  } while (ierr != 0);
258 
259  // Solve the system of equations.
260  Column_Vector xx = Ai * yy;
261 
262  // Extract the needed pieces from the inverted matrix
263  // and the solution vector.
264  W = Ai.sub(1, nconstraints, 1, nconstraints);
265  if (nbadvars > 0) {
266  U = Ai.sub(nconstraints + 1, nconstraints + nbadvars, nconstraints + 1, nconstraints + nbadvars);
267  V = Ai.sub(nconstraints + 1, nconstraints + nbadvars, 1, nconstraints);
268  d = xx.sub(nconstraints + 1, nconstraints + nbadvars);
269  }
270 
271  alpha = xx.sub(1, nconstraints);
272 
273  return true;
274  }
275 
276 } // unnamed namespace
277 
278 namespace hitfit {
279 
280  //*************************************************************************
281 
283  //
284  // Purpose: Constructor.
285  //
286  // Inputs:
287  // args - The parameter settings for this instance.
288  //
289  : Base_Constrainer(args.base_constrainer_args()), _args(args) {}
290 
291  double Chisq_Constrainer::fit(Constraint_Calculator& constraint_calculator,
292  const Column_Vector& xm,
293  Column_Vector& x,
294  const Column_Vector& ym,
295  Column_Vector& y,
296  const Matrix& G_i,
297  const Diagonal_Matrix& Y,
298  Column_Vector& pullx,
299  Column_Vector& pully,
300  Matrix& Q,
301  Matrix& R,
302  Matrix& S)
303  //
304  // Purpose: Do a constrained fit.
305  //
306  // Call the number of well-measured variables Nw, the number of
307  // poorly-measured variables Np, and the number of constraints Nc.
308  //
309  // Inputs:
310  // constraint_calculator - The object that will be used to evaluate
311  // the constraints.
312  // xm(Nw) - The measured values of the well-measured variables.
313  // ym(Np) - The measured values of the poorly-measured variables.
314  // x(Nw) - The starting values for the well-measured variables.
315  // y(Np) - The starting values for the poorly-measured variables.
316  // G_i(Nw,Nw) - The error matrix for the well-measured variables.
317  // Y(Np,Np) - The inverse error matrix for the poorly-measured variables.
318  //
319  // Outputs:
320  // x(Nw) - The fit values of the well-measured variables.
321  // y(Np) - The fit values of the poorly-measured variables.
322  // pullx(Nw) - The pull quantities for the well-measured variables.
323  // pully(Nw) - The pull quantities for the poorly-measured variables.
324  // Q(Nw,Nw) - The final error matrix for the well-measured variables.
325  // R(Np,Np) - The final error matrix for the poorly-measured variables.
326  // S(Nw,Np) - The final cross error matrix for the two sets of variables.
327  //
328  // Returns:
329  // The minimum chisq satisfying the constraints.
330  // Returns a value < 0 if the fit failed to converge.
331  //
332  {
333  // Check that the various matrices we've been passed have consistent
334  // dimensionalities.
335  int nvars = x.num_row();
336  assert(nvars == G_i.num_col());
337  assert(nvars == xm.num_row());
338 
339  int nbadvars = y.num_row();
340  assert(nbadvars == Y.num_col());
341  assert(nbadvars == ym.num_row());
342 
343  // If we're going to check the chisq calculation by explicitly using G,
344  // calculate it now from its inverse G_i.
345  Matrix G(nvars, nvars);
346  if (_args.use_G()) {
347  int ierr = 0;
348  G = G_i.inverse(ierr);
349  assert(!ierr);
350  }
351 
352  int nconstraints = constraint_calculator.nconstraints();
353 
354  // Results of the constraint evaluation function.
355  Row_Vector F(nconstraints); // Constraint vector.
356  Matrix Bx(nvars, nconstraints); // Gradients wrt x
357  Matrix By(nbadvars, nconstraints); // Gradients wrt y
358 
359  // (2) Evaluate the constraints at the starting point.
360  // If the starting point is rejected as invalid,
361  // give up and return an error.
362  if (!call_constraint_fcn(constraint_calculator, x, y, F, Bx, By)) {
363  // cout << "Bad initial values!";
364  // return -1000;
365  return -999.0;
366  }
367 
368  // (3) Initialize variables for the fitting loop.
369  double constraint_sum_last = -1000;
370  double chisq_last = -1000;
371  bool near_convergence = false;
372  double last_step_cutsize = 1;
373 
374  unsigned nit = 0;
375 
376  // Initialize the displacement vectors c and d.
377  Column_Vector c = x - xm;
378  Column_Vector d = y - ym;
379 
380  Matrix E(nvars, nconstraints);
381  Matrix W(nconstraints, nconstraints);
382  Matrix U(nbadvars, nbadvars);
383  Matrix V(nbadvars, nconstraints);
384 
385  // (4) Fitting loop:
386  do {
387  // (5) Calculate E, H, and r.
388  E = G_i * Bx;
389  Matrix H = E.T() * Bx;
390  Row_Vector r = c.T() * Bx + d.T() * By - F;
391 
392  // (6) Solve the linearized system for the new values
393  // of the Lagrange multipliers
394  // $\alpha$ and the new value for the displacements d.
395  Column_Vector alpha(nvars);
396  Column_Vector d1(nbadvars);
397  if (!solve_linear_system(H, Y, By, r, alpha, d1, W, U, V)) {
399  // return -1000;
400  return -998.0;
401  }
402 
403  // (7) Compute the new values for the displacements c and the chisq.
404  Column_Vector c1 = -E * alpha;
405  double chisq = -scalar(r * alpha);
406 
407  double psi_cut = 0;
408 
409  // (8) Find where this step is going to be taking us.
410  x = c1 + xm;
411  y = d1 + ym;
412 
413  // (9) Set up for cutting this step, should we have to.
414  Matrix save_By = By;
415  Row_Vector save_negF = -F;
416  double this_step_cutsize = 1;
417  double constraint_sum = -1;
418  unsigned ncut = 0;
419 
420  // (10) Evaluate the constraints at the new point.
421  // If the point is rejected, we have to try to cut the step.
422  // We accept the step if:
423  // The constraint sum is below the convergence threshold
424  // constraint_sum_eps, or
425  // This is the first iteration, or
426  // The constraint sum has decreased since the last iteration.
427  // Otherwise, the constraints have gotten worse, and we
428  // try to cut the step.
429  while (!call_constraint_fcn(constraint_calculator, x, y, F, Bx, By) ||
430  ((constraint_sum = norm_infinity(F)) > _args.constraint_sum_eps() && nit > 0 &&
431  constraint_sum > constraint_sum_last)) {
432  // Doing step cutting...
433  if (nit > 0 && _args.printfit() && ncut == 0) {
434  cout << "(" << chisq << " " << chisq_last << ") ";
435  }
436 
437  // (10a) If this is the first time we've tried to cut this step,
438  // test to see if the chisq is stationary. If it hasn't changed
439  // since the last iteration, try a directed step.
440  if (ncut == 0 && abs(chisq - chisq_last) < _args.chisq_diff_eps()) {
441  // Trying a directed step now.
442  // Try to make the smallest step which satisfies the
443  // (linearized) constraints.
444  if (_args.printfit())
445  cout << " directed step ";
446 
447  // (10a.i) Solve the linearized system for $\beta$ and
448  // the y-displacement vector $\delta$.
449  Column_Vector beta(nconstraints);
450  Column_Vector delta(nbadvars);
451  solve_linear_system(H, Y, save_By, save_negF, beta, delta, W, U, V);
452 
453  // (10a.ii) Get the x-displacement vector $\gamma$.
454  Column_Vector gamma = -E * beta;
455 
456  // (10a.iii) Find the destination of the directed step.
457  x = c + xm + gamma;
458  y = d + ym + delta;
459 
460  // (10a.iv) Accept this point if it's not rejected by the constraint
461  // function, and the constraints improve.
462  if (call_constraint_fcn(constraint_calculator, x, y, F, Bx, By) && (constraint_sum = norm_infinity(F)) > 0 &&
463  (constraint_sum < constraint_sum_last)) {
464  // Accept this step. Calculate the chisq and new displacement
465  // vectors.
466  chisq = chisq_last - scalar((-save_negF + r * 2) * beta);
467  c1 = x - xm;
468  d1 = y - ym;
469 
470  // Exit from step cutting loop.
471  break;
472  }
473  }
474 
475  // If this is the first time we're cutting the step,
476  // initialize $\psi$.
477  if (ncut == 0)
478  psi_cut = scalar((save_negF - r) * alpha);
479 
480  // (10b) Give up if we've tried to cut this step too many times.
481  if (++ncut > _args.max_cut()) {
482  // cout << " Too many cut steps ";
483  // return -1000;
484  return -997.0;
485  }
486 
487  // (10c) Set up the size by which we're going to cut this step.
488  // Normally, this is cutsize. But if this is the first time we're
489  // cutting this step and the last step was also cut, set the cut
490  // size to twice the final cut size from the last step (provided
491  // that it is less than cutsize).
492  double this_cutsize = _args.cutsize();
493  if (ncut == 1 && last_step_cutsize < 1) {
494  this_cutsize = 2 * last_step_cutsize;
495  if (this_cutsize > _args.cutsize())
496  this_cutsize = _args.cutsize();
497  }
498 
499  // (10d) Keep track of the total amount by which we've cut this step.
500  this_step_cutsize *= this_cutsize;
501 
502  // If it falls below min_tot_cutsize, give up.
503  if (this_step_cutsize < _args.min_tot_cutsize()) {
504  // cout << "Cut size underflow ";
505  // return -1000;
506  return -996.0;
507  }
508 
509  // (10e) Cut the step: calculate the new displacement vectors.
510  double cutleft = 1 - this_cutsize;
511  c1 = c1 * this_cutsize + c * cutleft;
512  d1 = d1 * this_cutsize + d * cutleft;
513 
514  // (10f) Calculate the new chisq.
515  if (chisq_last >= 0) {
516  chisq = this_cutsize * this_cutsize * chisq + cutleft * cutleft * chisq_last +
517  2 * this_cutsize * cutleft * psi_cut;
518  psi_cut = this_cutsize * psi_cut + cutleft * chisq_last;
519  } else
520  chisq = chisq_last;
521 
522  // Log what we've done.
523  if (_args.printfit()) {
524  cout << constraint_sum << " cut " << ncut << " size " << setiosflags(ios_base::scientific) << this_cutsize
525  << " tot size " << this_step_cutsize << resetiosflags(ios_base::scientific) << " " << chisq << "\n";
526  }
527 
528  // Find the new step destination.
529  x = c1 + xm;
530  y = d1 + ym;
531 
532  // Now, go and test the step again for acceptability.
533  }
534 
535  // (11) At this point, we have an acceptable step.
536  // Shuffle things around to prepare for the next step.
537  last_step_cutsize = this_step_cutsize;
538 
539  // If requested, calculate the chisq using G to test for
540  // possible loss of precision.
541  double chisq_b = 0;
542  if (_args.use_G()) {
543  chisq_b = scalar(c1.T() * G * c1) + scalar(d1.T() * Y * d1);
544  if (chisq >= 0 && abs((chisq - chisq_b) / chisq) > _args.chisq_test_eps()) {
545  cout << chisq << " " << chisq_b << "lost precision?\n";
546  abort();
547  }
548  }
549 
550  // Log what we're doing.
551  if (_args.printfit()) {
552  cout << chisq << " ";
553  if (_args.use_G())
554  cout << chisq_b << " ";
555  }
556 
557  double z2 = abs(chisq - chisq_last);
558 
559  if (_args.printfit()) {
560  cout << constraint_sum << " " << z2 << "\n";
561  }
562 
563  c = c1;
564  d = d1;
565  chisq_last = chisq;
566  constraint_sum_last = constraint_sum;
567 
568  // (12) Test for convergence. The conditions must be satisfied
569  // for two iterations in a row.
570  if (chisq >= 0 && constraint_sum < _args.constraint_sum_eps() && z2 < _args.chisq_diff_eps()) {
571  if (near_convergence)
572  break; // Converged! Exit loop.
573  near_convergence = true;
574  } else
575  near_convergence = false;
576 
577  // (13) Give up if we've done this too many times.
578  if (++nit > _args.maxit()) {
579  // cout << "too many iterations";
580  // return -1000;
581  return -995.0;
582  }
583 
584  } while (true);
585 
586  // (15) Ok, we have a successful fit!
587 
588  // Calculate the error matrices.
589  Q = E * W * E.T();
590  S = -E * V.T();
591  R = U;
592 
593  // And the vectors of pull functions.
594  pullx = Column_Vector(nvars);
595  for (int i = 1; i <= nvars; i++) {
596  double a = Q(i, i);
597  if (a < 0)
598  pullx(i) = c(i) / sqrt(-a);
599  else {
600  pullx(i) = 0;
601  // cout << " bad pull fcn for var " << i << " (" << a << ") ";
602  }
603  }
604 
605  pully = Column_Vector(nbadvars);
606  for (int i = 1; i <= nbadvars; i++) {
607  double a = 1 - Y(i, i) * R(i, i);
608  if (a > 0)
609  pully(i) = d(i) * sqrt(Y(i, i) / a);
610  else {
611  pully(i) = 0;
612  // cout << " bad pull fcn for badvar " << i << " ";
613  }
614  }
615 
616  // Finish calculation of Q.
617  Q = Q + G_i;
618 
619  // Return the final chisq.
620  return chisq_last;
621  }
622 
623  std::ostream& Chisq_Constrainer::print(std::ostream& s) const
624  //
625  // Purpose: Print our state.
626  //
627  // Inputs:
628  // s - The stream to which to write.
629  //
630  // Returns:
631  // The stream S.
632  //
633  {
635  s << " printfit: " << _args.printfit() << " use_G: " << _args.use_G() << "\n";
636  s << " constraint_sum_eps: " << _args.constraint_sum_eps() << " chisq_diff_eps: " << _args.chisq_diff_eps()
637  << " chisq_test_eps: " << _args.chisq_test_eps() << "\n";
638  s << " maxit: " << _args.maxit() << " max_cut: " << _args.max_cut()
639  << " min_tot_cutsize: " << _args.min_tot_cutsize() << " cutsize: " << _args.cutsize() << "\n";
640  return s;
641  }
642 
643 } // namespace hitfit
const Base_Constrainer_Args & base_constrainer_args() const
Define an abstract interface for getting parameter settings.
CLHEP::HepVector Column_Vector
Definition: matutil.h:63
Abstract base class for evaluating constraints. Users derive from this class and implement the eval()...
CLHEP::HepMatrix Matrix
Definition: matutil.h:62
Chisq_Constrainer(const Chisq_Constrainer_Args &args)
Base class for constrained fitter.
Chisq_Constrainer_Args(const Defaults &defs)
Constructor, creates an instance of Chisq_Constrainer_Args from a Defaults object.
T sqrt(T t)
Definition: SSEVec.h:19
const Base_Constrainer_Args _base_constrainer_args
Minimize a subject to a set of constraints, based on the SQUAW algorithm.
double fit(Constraint_Calculator &constraint_calculator, const Column_Vector &xm, Column_Vector &x, const Column_Vector &ym, Column_Vector &y, const Matrix &G_i, const Diagonal_Matrix &Y, Column_Vector &pullx, Column_Vector &pully, Matrix &Q, Matrix &R, Matrix &S) override
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
bool call_constraint_fcn(Constraint_Calculator &constraint_calculator, const Column_Vector &x, const Column_Vector &y, Row_Vector &F, Matrix &Bx, Matrix &By) const
Helper function to evaluate constraints. This takes care of checking what the user function returns a...
Hold on to parameters for the Chisq_Constrainer class.
d
Definition: ztail.py:151
std::ostream & print(std::ostream &s) const override
Print the state of this instance of Chisq_Constrainer.
Hold on to parameters for the Base_Constrainer class.
CLHEP::HepDiagMatrix Diagonal_Matrix
Definition: matutil.h:64
virtual double get_float(std::string name) const =0
virtual bool get_bool(std::string name) const =0
Define an interface for getting parameter settings.
Definition: Defaults.h:57
double a
Definition: hdecay.h:119
const Chisq_Constrainer_Args _args
alpha
zGenParticlesMatch = cms.InputTag(""),
virtual int get_int(std::string name) const =0
static uInt32 F(BLOWFISH_CTX *ctx, uInt32 x)
Definition: blowfish.cc:163
virtual std::ostream & print(std::ostream &s) const
Print out internal state to output stream.
Row-vector class. CLHEP doesn&#39;t have a row-vector class, so HitFit uses its own. This is only a simpl...
Definition: matutil.h:79
double scalar(const CLHEP::HepGenMatrix &m)
Return the matrix as a scalar. Raise an assertion if the matris is not .
Definition: matutil.cc:166