49 using std::resetiosflags;
50 using std::setiosflags;
64 : _base_constrainer_args(defs) {
190 bool solve_linear_system(
const Matrix&
H,
217 int nconstraints = H.num_row();
218 int nbadvars = Y.num_row();
221 Matrix A(nconstraints + nbadvars, nconstraints + nbadvars);
224 A.sub(nconstraints + 1, nconstraints + 1, Y);
225 A.sub(1, nconstraints + 1, By.T());
226 A.sub(nconstraints + 1, 1, By);
238 Ai = A.inverse(ierr);
241 for (
int i = 1;
i <= nconstraints;
i++) {
243 for (
int j = 1;
j <= nconstraints;
j++) {
264 W = Ai.sub(1, nconstraints, 1, nconstraints);
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);
271 alpha = xx.sub(1, nconstraints);
335 int nvars = x.num_row();
336 assert(nvars == G_i.num_col());
337 assert(nvars == xm.num_row());
339 int nbadvars = y.num_row();
340 assert(nbadvars == Y.num_col());
341 assert(nbadvars == ym.num_row());
348 G = G_i.inverse(ierr);
352 int nconstraints = constraint_calculator.
nconstraints();
356 Matrix Bx(nvars, nconstraints);
357 Matrix By(nbadvars, nconstraints);
369 double constraint_sum_last = -1000;
370 double chisq_last = -1000;
371 bool near_convergence =
false;
372 double last_step_cutsize = 1;
380 Matrix E(nvars, nconstraints);
381 Matrix W(nconstraints, nconstraints);
383 Matrix V(nbadvars, nconstraints);
397 if (!solve_linear_system(H, Y, By, r, alpha, d1, W, U, V)) {
405 double chisq = -
scalar(r * alpha);
416 double this_step_cutsize = 1;
417 double constraint_sum = -1;
431 constraint_sum > constraint_sum_last)) {
434 cout <<
"(" << chisq <<
" " << chisq_last <<
") ";
445 cout <<
" directed step ";
451 solve_linear_system(H, Y, save_By, save_negF, beta, delta, W, U, V);
462 if (
call_constraint_fcn(constraint_calculator, x, y, F, Bx, By) && (constraint_sum = norm_infinity(F)) > 0 &&
463 (constraint_sum < constraint_sum_last)) {
466 chisq = chisq_last -
scalar((-save_negF + r * 2) * beta);
478 psi_cut =
scalar((save_negF - r) * alpha);
493 if (ncut == 1 && last_step_cutsize < 1) {
494 this_cutsize = 2 * last_step_cutsize;
500 this_step_cutsize *= this_cutsize;
510 double cutleft = 1 - this_cutsize;
511 c1 = c1 * this_cutsize + c * cutleft;
512 d1 = d1 * this_cutsize + d * cutleft;
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;
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";
537 last_step_cutsize = this_step_cutsize;
545 cout << chisq <<
" " << chisq_b <<
"lost precision?\n";
552 cout << chisq <<
" ";
554 cout << chisq_b <<
" ";
557 double z2 =
abs(chisq - chisq_last);
560 cout << constraint_sum <<
" " << z2 <<
"\n";
566 constraint_sum_last = constraint_sum;
571 if (near_convergence)
573 near_convergence =
true;
575 near_convergence =
false;
595 for (
int i = 1;
i <= nvars;
i++) {
606 for (
int i = 1;
i <= nbadvars;
i++) {
607 double a = 1 -
Y(
i,
i) *
R(
i,
i);
const Base_Constrainer_Args & base_constrainer_args() const
Define an abstract interface for getting parameter settings.
double min_tot_cutsize() const
CLHEP::HepVector Column_Vector
Abstract base class for evaluating constraints. Users derive from this class and implement the eval()...
Chisq_Constrainer(const Chisq_Constrainer_Args &args)
double constraint_sum_eps() const
Base class for constrained fitter.
Chisq_Constrainer_Args(const Defaults &defs)
Constructor, creates an instance of Chisq_Constrainer_Args from a Defaults object.
double chisq_diff_eps() const
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)
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.
std::ostream & print(std::ostream &s) const override
Print the state of this instance of Chisq_Constrainer.
double chisq_test_eps() const
Hold on to parameters for the Base_Constrainer class.
CLHEP::HepDiagMatrix Diagonal_Matrix
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.
const Chisq_Constrainer_Args _args
alpha
zGenParticlesMatch = cms.InputTag(""),
virtual int get_int(std::string name) const =0
double _constraint_sum_eps
static uInt32 F(BLOWFISH_CTX *ctx, uInt32 x)
virtual std::ostream & print(std::ostream &s) const
Print out internal state to output stream.
Row-vector class. CLHEP doesn't have a row-vector class, so HitFit uses its own. This is only a simpl...
double scalar(const CLHEP::HepGenMatrix &m)
Return the matrix as a scalar. Raise an assertion if the matris is not .