|
|
|
Minimize a 
subject to a set of constraints. Based on the SQUAW algorithm.
More...
#include <Chisq_Constrainer.h>
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. | |
| virtual | ~Chisq_Constrainer () |
Private Attributes | |
| const Chisq_Constrainer_Args | _args |
Minimize a subject to a set of constraints. Based on the SQUAW algorithm.
Definition at line 253 of file Chisq_Constrainer.h.
| hitfit::Chisq_Constrainer::Chisq_Constrainer | ( | const Chisq_Constrainer_Args & | args | ) |
Constructor. Create an instance of Chisq_Constrainer from a Chisq_Constrainer_Args object.
| 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 | ( | ) | [inline, virtual] |
| 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.
| 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. |
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.
| 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;
}
const Chisq_Constrainer_Args hitfit::Chisq_Constrainer::_args [private] |
Parameter settings for this instance of Chisq_Constrainer.
Reimplemented from hitfit::Base_Constrainer.
Definition at line 377 of file Chisq_Constrainer.h.
Referenced by fit().
1.7.3