hitfit::Chisq_Constrainer Class Reference

Minimize a subject to a set of constraints. Based on the SQUAW algorithm. More...

#include <Chisq_Constrainer.h>

Inheritance diagram for hitfit::Chisq_Constrainer:

Public Member Functions
	Chisq_Constrainer (const Chisq_Constrainer_Args &args)
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
std::ostream &	print (std::ostream &s) const override
	Print the state of this instance of Chisq_Constrainer. More...
	~Chisq_Constrainer () override
Public Member Functions inherited from hitfit::Base_Constrainer
	Base_Constrainer (const Base_Constrainer_Args &args)
virtual	~Base_Constrainer ()

Private Attributes
const Chisq_Constrainer_Args	_args

Additional Inherited Members
Protected Member Functions inherited from hitfit::Base_Constrainer
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 againts numerical derivatives, it that was requested. More...

Detailed Description

Minimize a subject to a set of constraints. Based on the SQUAW algorithm.

Definition at line 244 of file Chisq_Constrainer.h.

Constructor & Destructor Documentation

Chisq_Constrainer()

hitfit::Chisq_Constrainer::Chisq_Constrainer

(

const Chisq_Constrainer_Args &

args

)

Constructor. Create an instance of Chisq_Constrainer from a Chisq_Constrainer_Args object.

Parameters

args	The parameter settings for this instance.

Definition at line 282 of file Chisq_Constrainer.cc.

      : Base_Constrainer(args.base_constrainer_args()), _args(args) {}

~Chisq_Constrainer()

hitfit::Chisq_Constrainer::~Chisq_Constrainer ( )

inlineoverride

Destructor.

Definition at line 265 of file Chisq_Constrainer.h.

{}

Member Function Documentation

fit()

double hitfit::Chisq_Constrainer::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  )

overridevirtual

@ brief Do a constrained fit. Call the number of well-measured variables Nw, the number of poorly-measured variables Np, and the number of constraints Nc.

Parameters

constraint_calculator	The object that will be used to evaluate the constraints.
xm	Measured values of the well-measured variables, has dimension Nw.
x	Before the fit: starting value of the well-measured variables for the fit, has dimension Nw. After the fit: final values of the well-measured variables.
ym	Measured values of the poorly-measured variables, has dimension Np.
y	Before the fit: starting value of the poorly-measured variables for the fit, has dimension Np. After the fit: final values of the poorly-measured variables.
G_i	Error matrix for the well-measured variables, has dimension Nw,Nw.
Y	Inverse error matrix for the poorly-measured variables, has dimension Np,Np.
pullx	Pull quantities for the well-measured variables, has dimension Nw.
pully	Pull quantities for the poorly-measured variables, has dimension Np.
Q	Error matrix for the well-measured variables, has dimension Nw,Nw.
R	Error matrix for the poorly-measured variables, has dimension Np,Np.
S	Error matrix for the correlation between well-measured variables and poorly-measured variables, has dimension Nw,Np.

Input:

constraint_calculator, xm, x, ym, y, G_i, y..

Output:

x, y, pullx, pully, Q, R, S.

Return:

of the fit. Should returns a negative value if the fit does not converge.

cout << "singular matrix!";

Implements hitfit::Base_Constrainer.

Definition at line 291 of file Chisq_Constrainer.cc.

References _args, a, funct::abs(), alpha, cms::cuda::assert(), HLT_2022v15_cff::beta, HltBtagPostValidation_cff::c, alignmentValidation::c1, hitfit::Base_Constrainer::call_constraint_fcn(), hitfit::Chisq_Constrainer_Args::chisq_diff_eps(), hitfit::Chisq_Constrainer_Args::chisq_test_eps(), hitfit::Chisq_Constrainer_Args::constraint_sum_eps(), gather_cfg::cout, hitfit::Chisq_Constrainer_Args::cutsize(), ztail::d, d1, dumpMFGeometry_cfg::delta, F(), cmssw_cycle_finder::G, CustomPhysics_cfi::gamma, data-class-funcs::H, mps_fire::i, hitfit::Chisq_Constrainer_Args::max_cut(), hitfit::Chisq_Constrainer_Args::maxit(), hitfit::Chisq_Constrainer_Args::min_tot_cutsize(), hitfit::Constraint_Calculator::nconstraints(), hitfit::Chisq_Constrainer_Args::printfit(), dttmaxenums::R, hitfit::scalar(), mathSSE::sqrt(), mitigatedMETSequence_cff::U, hitfit::Chisq_Constrainer_Args::use_G(), cms::cuda::V, x, BeamSpotPI::Y, and testProducerWithPsetDescEmpty_cfi::z2.

Referenced by trackingPlots.Iteration::modules().

  {
    // Check that the various matrices we've been passed have consistent
    // dimensionalities.
    int nvars = x.num_row();
    assert(nvars == G_i.num_col());
    assert(nvars == xm.num_row());

    int nbadvars = y.num_row();
    assert(nbadvars == Y.num_col());
    assert(nbadvars == ym.num_row());

    // If we're going to check the chisq calculation by explicitly using G,
    // calculate it now from its inverse G_i.
    Matrix G(nvars, nvars);
    if (_args.use_G()) {
      int ierr = 0;
      G = G_i.inverse(ierr);
      assert(!ierr);
    }

    int nconstraints = constraint_calculator.nconstraints();

    // Results of the constraint evaluation function.
    Row_Vector F(nconstraints);         // Constraint vector.
    Matrix Bx(nvars, nconstraints);     // Gradients wrt x
    Matrix By(nbadvars, nconstraints);  // Gradients wrt y

    // (2) Evaluate the constraints at the starting point.
    // If the starting point is rejected as invalid,
    // give up and return an error.
    if (!call_constraint_fcn(constraint_calculator, x, y, F, Bx, By)) {
      //    cout << "Bad initial values!";
      //    return -1000;
      return -999.0;
    }

    // (3) Initialize variables for the fitting loop.
    double constraint_sum_last = -1000;
    double chisq_last = -1000;
    bool near_convergence = false;
    double last_step_cutsize = 1;

    unsigned nit = 0;

    // Initialize the displacement vectors c and d.
    Column_Vector c = x - xm;
    Column_Vector d = y - ym;

    Matrix E(nvars, nconstraints);
    Matrix W(nconstraints, nconstraints);
    Matrix U(nbadvars, nbadvars);
    Matrix V(nbadvars, nconstraints);

    // (4) Fitting loop:
    do {
      // (5) Calculate E, H, and r.
      E = G_i * Bx;
      Matrix H = E.T() * Bx;
      Row_Vector r = c.T() * Bx + d.T() * By - F;

      // (6) Solve the linearized system for the new values
      // of the Lagrange multipliers
      // $\alpha$ and the new value for the displacements d.
      Column_Vector alpha(nvars);
      Column_Vector d1(nbadvars);
      if (!solve_linear_system(H, Y, By, r, alpha, d1, W, U, V)) {
        //      return -1000;
        return -998.0;
      }

      // (7) Compute the new values for the displacements c and the chisq.
      Column_Vector c1 = -E * alpha;
      double chisq = -scalar(r * alpha);

      double psi_cut = 0;

      // (8) Find where this step is going to be taking us.
      x = c1 + xm;
      y = d1 + ym;

      // (9) Set up for cutting this step, should we have to.
      Matrix save_By = By;
      Row_Vector save_negF = -F;
      double this_step_cutsize = 1;
      double constraint_sum = -1;
      unsigned ncut = 0;

      // (10) Evaluate the constraints at the new point.
      // If the point is rejected, we have to try to cut the step.
      // We accept the step if:
      //  The constraint sum is below the convergence threshold
      //    constraint_sum_eps, or
      //  This is the first iteration, or
      //  The constraint sum has decreased since the last iteration.
      // Otherwise, the constraints have gotten worse, and we
      // try to cut the step.
      while (!call_constraint_fcn(constraint_calculator, x, y, F, Bx, By) ||
             ((constraint_sum = norm_infinity(F)) > _args.constraint_sum_eps() && nit > 0 &&
              constraint_sum > constraint_sum_last)) {
        // Doing step cutting...
        if (nit > 0 && _args.printfit() && ncut == 0) {
          cout << "(" << chisq << " " << chisq_last << ") ";
        }

        // (10a) If this is the first time we've tried to cut this step,
        // test to see if the chisq is stationary.  If it hasn't changed
        // since the last iteration, try a directed step.
        if (ncut == 0 && abs(chisq - chisq_last) < _args.chisq_diff_eps()) {
          // Trying a directed step now.
          // Try to make the smallest step which satisfies the
          // (linearized) constraints.
          if (_args.printfit())
            cout << " directed step ";

          // (10a.i) Solve the linearized system for $\beta$ and
          // the y-displacement vector $\delta$.
          Column_Vector beta(nconstraints);
          Column_Vector delta(nbadvars);
          solve_linear_system(H, Y, save_By, save_negF, beta, delta, W, U, V);

          // (10a.ii) Get the x-displacement vector $\gamma$.
          Column_Vector gamma = -E * beta;

          // (10a.iii) Find the destination of the directed step.
          x = c + xm + gamma;
          y = d + ym + delta;

          // (10a.iv) Accept this point if it's not rejected by the constraint
          // function, and the constraints improve.
          if (call_constraint_fcn(constraint_calculator, x, y, F, Bx, By) && (constraint_sum = norm_infinity(F)) > 0 &&
              (constraint_sum < constraint_sum_last)) {
            // Accept this step.  Calculate the chisq and new displacement
            // vectors.
            chisq = chisq_last - scalar((-save_negF + r * 2) * beta);
            c1 = x - xm;
            d1 = y - ym;

            // Exit from step cutting loop.
            break;
          }
        }

        // If this is the first time we're cutting the step,
        // initialize $\psi$.
        if (ncut == 0)
          psi_cut = scalar((save_negF - r) * alpha);

        // (10b) Give up if we've tried to cut this step too many times.
        if (++ncut > _args.max_cut()) {
          //    cout << " Too many cut steps ";
          //    return -1000;
          return -997.0;
        }

        // (10c) Set up the size by which we're going to cut this step.
        // Normally, this is cutsize.  But if this is the first time we're
        // cutting this step and the last step was also cut, set the cut
        // size to twice the final cut size from the last step (provided
        // that it is less than cutsize).
        double this_cutsize = _args.cutsize();
        if (ncut == 1 && last_step_cutsize < 1) {
          this_cutsize = 2 * last_step_cutsize;
          if (this_cutsize > _args.cutsize())
            this_cutsize = _args.cutsize();
        }

        // (10d) Keep track of the total amount by which we've cut this step.
        this_step_cutsize *= this_cutsize;

        // If it falls below min_tot_cutsize, give up.
        if (this_step_cutsize < _args.min_tot_cutsize()) {
          //    cout << "Cut size underflow ";
          //    return -1000;
          return -996.0;
        }

        // (10e) Cut the step: calculate the new displacement vectors.
        double cutleft = 1 - this_cutsize;
        c1 = c1 * this_cutsize + c * cutleft;
        d1 = d1 * this_cutsize + d * cutleft;

        // (10f) Calculate the new chisq.
        if (chisq_last >= 0) {
          chisq = this_cutsize * this_cutsize * chisq + cutleft * cutleft * chisq_last +
                  2 * this_cutsize * cutleft * psi_cut;
          psi_cut = this_cutsize * psi_cut + cutleft * chisq_last;
        } else
          chisq = chisq_last;

        // Log what we've done.
        if (_args.printfit()) {
          cout << constraint_sum << " cut " << ncut << " size " << setiosflags(ios_base::scientific) << this_cutsize
               << " tot size " << this_step_cutsize << resetiosflags(ios_base::scientific) << " " << chisq << "\n";
        }

        // Find the new step destination.
        x = c1 + xm;
        y = d1 + ym;

        // Now, go and test the step again for acceptability.
      }

      // (11) At this point, we have an acceptable step.
      // Shuffle things around to prepare for the next step.
      last_step_cutsize = this_step_cutsize;

      // If requested, calculate the chisq using G to test for
      // possible loss of precision.
      double chisq_b = 0;
      if (_args.use_G()) {
        chisq_b = scalar(c1.T() * G * c1) + scalar(d1.T() * Y * d1);
        if (chisq >= 0 && abs((chisq - chisq_b) / chisq) > _args.chisq_test_eps()) {
          cout << chisq << " " << chisq_b << "lost precision?\n";
          abort();
        }
      }

      // Log what we're doing.
      if (_args.printfit()) {
        cout << chisq << " ";
        if (_args.use_G())
          cout << chisq_b << " ";
      }

      double z2 = abs(chisq - chisq_last);

      if (_args.printfit()) {
        cout << constraint_sum << " " << z2 << "\n";
      }

      c = c1;
      d = d1;
      chisq_last = chisq;
      constraint_sum_last = constraint_sum;

      // (12) Test for convergence.  The conditions must be satisfied
      // for two iterations in a row.
      if (chisq >= 0 && constraint_sum < _args.constraint_sum_eps() && z2 < _args.chisq_diff_eps()) {
        if (near_convergence)
          break;  // Converged!  Exit loop.
        near_convergence = true;
      } else
        near_convergence = false;

      // (13) Give up if we've done this too many times.
      if (++nit > _args.maxit()) {
        //      cout << "too many iterations";
        //      return -1000;
        return -995.0;
      }

    } while (true);

    // (15) Ok, we have a successful fit!

    // Calculate the error matrices.
    Q = E * W * E.T();
    S = -E * V.T();
    R = U;

    // And the vectors of pull functions.
    pullx = Column_Vector(nvars);
    for (int i = 1; i <= nvars; i++) {
      double a = Q(i, i);
      if (a < 0)
        pullx(i) = c(i) / sqrt(-a);
      else {
        pullx(i) = 0;
        //      cout << " bad pull fcn for var " << i << " (" << a << ") ";
      }
    }

    pully = Column_Vector(nbadvars);
    for (int i = 1; i <= nbadvars; i++) {
      double a = 1 - Y(i, i) * R(i, i);
      if (a > 0)
        pully(i) = d(i) * sqrt(Y(i, i) / a);
      else {
        pully(i) = 0;
        //      cout << " bad pull fcn for badvar " << i << " ";
      }
    }

    // Finish calculation of Q.
    Q = Q + G_i;

    // Return the final chisq.
    return chisq_last;
  }

print()

std::ostream & hitfit::Chisq_Constrainer::print ( std::ostream &  s ) const

overridevirtual

Print the state of this instance of Chisq_Constrainer.

Parameters

s	The output stream to which the output is sent.

Reimplemented from hitfit::Base_Constrainer.

Definition at line 623 of file Chisq_Constrainer.cc.

References hitfit::Base_Constrainer::print(), and alignCSCRings::s.

  {
    Base_Constrainer::print(s);
    s << " printfit: " << _args.printfit() << "  use_G: " << _args.use_G() << "\n";
    s << " constraint_sum_eps: " << _args.constraint_sum_eps() << "  chisq_diff_eps: " << _args.chisq_diff_eps()
      << "  chisq_test_eps: " << _args.chisq_test_eps() << "\n";
    s << " maxit: " << _args.maxit() << "  max_cut: " << _args.max_cut()
      << "  min_tot_cutsize: " << _args.min_tot_cutsize() << "  cutsize: " << _args.cutsize() << "\n";
    return s;
  }

Member Data Documentation

_args

const Chisq_Constrainer_Args hitfit::Chisq_Constrainer::_args

private

Parameter settings for this instance of Chisq_Constrainer.

Definition at line 366 of file Chisq_Constrainer.h.

Referenced by fit().