Minimize a $\chi^{2}$ 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 $\chi^{2}$ subject to a set of constraints. Based on the SQUAW algorithm.

Definition at line 252 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 307 of file Chisq_Constrainer.cc.

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

virtual hitfit::Chisq_Constrainer::~Chisq_Constrainer ( )

inlinevirtual

Destructor.

Definition at line 274 of file Chisq_Constrainer.h.

274 {}

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:: $\chi^{2}$ of the fit. Should returns a negative value if the fit does not converge.

cout << "singular matrix!";

Implements hitfit::Base_Constrainer.

Definition at line 320 of file Chisq_Constrainer.cc.

 {
   // 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 673 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 376 of file Chisq_Constrainer.h.

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

Public Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation