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)
virtual 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)
virtual std::ostream &	print (std::ostream &s) const
	Print the state of this instance of Chisq_Constrainer. More...
virtual	~Chisq_Constrainer ()
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 253 of file Chisq_Constrainer.h.

Constructor & Destructor Documentation

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 308 of file Chisq_Constrainer.cc.

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

virtual hitfit::Chisq_Constrainer::~Chisq_Constrainer ( )

inlinevirtual

Destructor.

Definition at line 275 of file Chisq_Constrainer.h.

{}

Member Function Documentation

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  )

virtual

@ 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 321 of file Chisq_Constrainer.cc.

References _args, a, abs, alpha, beta, trackerHits::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(), delta, F(), 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(), and hitfit::Chisq_Constrainer_Args::use_G().

{
  // 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 (1);

  // (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;
}

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

virtual

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 674 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

const Chisq_Constrainer_Args hitfit::Chisq_Constrainer::_args

private

Parameter settings for this instance of Chisq_Constrainer.

Definition at line 377 of file Chisq_Constrainer.h.

Referenced by Vispa.Main.Application.Application::_readCommandLineAttributes(), and fit().