Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:31:16

0001 //
0002 //
0003 // File: src/Base_Constrainer.cc
0004 // Purpose: Abstract base for the chisq fitter classes.
0005 //          This allows for different algorithms to be used.
0006 // Created: Jul, 2000, sss, based on run 1 mass analysis code.
0007 //
0008 // CMSSW File      : src/Base_Constrainer.cc
0009 // Original Author : Scott Stuart Snyder <snyder@bnl.gov> for D0
0010 // Imported to CMSSW by Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>
0011 //
0012 
0013 /**
0014 
0015     @file Base_Constrainer.cc
0016 
0017     @brief Abstract base classes for the \f$\chi^{2}\f$
0018     fitter classes.  Includes helper function(s).
0019     See the documentation for header file Base_Constrainer.h for details.
0020 
0021     @par Creation date:
0022     July 2000.
0023 
0024     @author
0025     Scott Stuart Snyder <snyder@bnl.gov>.
0026 
0027     @par Modification History:
0028     Apr 2009: Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>:
0029     Imported to CMSSW.<br>
0030     Oct 2009: Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>:
0031     Added Doxygen tags for automatic generation of documentation.
0032 
0033     @par Terms of Usage:
0034     With consent from the original author (Scott Snyder).
0035 
0036  */
0037 
0038 #include "TopQuarkAnalysis/TopHitFit/interface/Base_Constrainer.h"
0039 #include "TopQuarkAnalysis/TopHitFit/interface/matutil.h"
0040 #include "TopQuarkAnalysis/TopHitFit/interface/Defaults.h"
0041 #include <iostream>
0042 #include <cmath>
0043 #include <cstdlib>
0044 
0045 using std::abort;
0046 using std::abs;
0047 using std::cout;
0048 using std::ostream;
0049 
0050 //*************************************************************************
0051 // Helper function for doing gradient testing.
0052 //
0053 
0054 namespace {
0055 
0056   /**
0057    @brief Test if <i>a</i> is significantly different than <i>b</i>.
0058    <i>c</i> sets the scale for the comparison, <i>eps</i> gives by
0059    how much they may differ.
0060 
0061    @param a First number to be compared.
0062    @param b Second number to be compared.
0063    @param c Scale of the comparison.
0064    @param eps How much the two numbers may be different.
0065    @par Return:
0066    Returns true if one of the two numbers is significantly larger than
0067    the other one, otherwise false.
0068  */
0069   bool test_different(double a, double b, double c, double eps)
0070   //
0071   // Purpose: Test if A is significantly different than B.
0072   //          C sets the scale for the comparison; EPS gives
0073   //          by how much they may differ.
0074   //
0075   {
0076     double scale = eps * (abs(a) + abs(b) + abs(c));
0077     if (abs(a) != 0 && abs(b) / abs(a) < 0.1 && abs(a) > scale)
0078       scale = abs(a) * .5;
0079     if (scale == 0)
0080       return false;
0081     if (scale < eps)
0082       scale = eps;
0083     return abs(a - b) > scale;
0084   }
0085 
0086 }  // unnamed namespace
0087 
0088 //*************************************************************************
0089 
0090 namespace hitfit {
0091 
0092   //*************************************************************************
0093 
0094   Base_Constrainer_Args::Base_Constrainer_Args(const Defaults& defs)
0095       //
0096       // Purpose: Constructor.
0097       //
0098       // Inputs:
0099       //   defs -        The Defaults instance from which to initialize.
0100       //
0101       : _test_gradient(defs.get_bool("test_gradient")),
0102         _test_step(defs.get_float("test_step")),
0103         _test_eps(defs.get_float("test_eps")) {}
0104 
0105   bool Base_Constrainer_Args::test_gradient() const
0106   //
0107   // Purpose: Return the test_gradient parameter.
0108   //          See the header for documentation.
0109   //
0110   {
0111     return _test_gradient;
0112   }
0113 
0114   double Base_Constrainer_Args::test_step() const
0115   //
0116   // Purpose: Return the test_step parameter.
0117   //          See the header for documentation.
0118   //
0119   {
0120     return _test_step;
0121   }
0122 
0123   double Base_Constrainer_Args::test_eps() const
0124   //
0125   // Purpose: Return the test_eps parameter.
0126   //          See the header for documentation.
0127   //
0128   {
0129     return _test_eps;
0130   }
0131 
0132   //*************************************************************************
0133 
0134   Constraint_Calculator::Constraint_Calculator(int nconstraints)
0135       //
0136       // Purpose: Constructor.
0137       //
0138       // Inputs:
0139       //   nconstraints- The number of constraint functions.
0140       //
0141       : _nconstraints(nconstraints) {}
0142 
0143   int Constraint_Calculator::nconstraints() const
0144   //
0145   // Purpose: Return the number of constraint functions.
0146   //
0147   // Returns:
0148   //   The number of constraint functions.
0149   //
0150   {
0151     return _nconstraints;
0152   }
0153 
0154   //*************************************************************************
0155 
0156   Base_Constrainer::Base_Constrainer(const Base_Constrainer_Args& args)
0157       //
0158       // Purpose: Constructor.
0159       //
0160       // Inputs:
0161       //   args -        The parameter settings for this instance.
0162       //
0163       : _args(args) {}
0164 
0165   std::ostream& Base_Constrainer::print(std::ostream& s) const
0166   //
0167   // Purpose: Print our state.
0168   //
0169   // Inputs:
0170   //   s -           The stream to which to write.
0171   //
0172   // Returns:
0173   //   The stream S.
0174   //
0175   {
0176     s << "Base_Constrainer parameters:\n";
0177     s << " test_gradient: " << _args.test_gradient() << " test_step: " << _args.test_step()
0178       << " test_eps: " << _args.test_eps() << "\n";
0179     return s;
0180   }
0181 
0182   /**
0183     @brief Output stream operator, print the content of this Base_Constrainer
0184     to an output stream.
0185 
0186     @param s The output stream to which to write.
0187 
0188     @param f The instance of Base_Constrainer to be printed.
0189 
0190 */
0191   std::ostream& operator<<(std::ostream& s, const Base_Constrainer& f)
0192   //
0193   // Purpose: Print our state.
0194   //
0195   // Inputs:
0196   //   s -           The stream to which to write.
0197   //   f -           The instance to dump.
0198   //
0199   // Returns:
0200   //   The stream S.
0201   //
0202   {
0203     return f.print(s);
0204   }
0205 
0206   bool Base_Constrainer::call_constraint_fcn(Constraint_Calculator& constraint_calculator,
0207                                              const Column_Vector& x,
0208                                              const Column_Vector& y,
0209                                              Row_Vector& F,
0210                                              Matrix& Bx,
0211                                              Matrix& By) const
0212   //
0213   // Purpose: Call the constraint function for the point x, y.
0214   //          Return F, Bx, By, and a flag saying if the
0215   //          point is acceptable.
0216   //
0217   //          If test_gradient is on, we verify the gradients returned
0218   //          by also computing them numerically.
0219   //
0220   // Inputs:
0221   //   constraints - The user-supplied object to evaluate the constraints.
0222   //   x(Nw) -       Vector of well-measured quantities where we evaluate
0223   //                 the constraints.
0224   //   y(Np) -       Vector of poorly-measured quantities where we evaluate
0225   //                 the constraints.
0226   //
0227   // Outputs:
0228   //   F(Nc) -       The results of the constraint functions.
0229   //   Bx(Nw,Nc) -   Gradients of F with respect to x.
0230   //   By(Np,Nc) -   Gradients of F with respect to y.
0231   //
0232   // Returns:
0233   //   True if the point is accepted, false if it was rejected.
0234   {
0235     // Call the user's function.
0236     bool val = constraint_calculator.eval(x, y, F, Bx, By);
0237 
0238     // If we're not doing gradients numerically, we're done.
0239     if (!_args.test_gradient())
0240       return val;
0241 
0242     // Bail if the point was rejected.
0243     if (!val)
0244       return false;
0245 
0246     int Nw = x.num_row();
0247     int Np = y.num_row();
0248     int Nc = F.num_col();
0249 
0250     // Numerically check Bx.
0251     for (int i = 1; i <= Nc; i++) {
0252       // Step a little along variable I.
0253       Column_Vector step_x(Nw, 0);
0254       step_x(i) = _args.test_step();
0255       Column_Vector new_x = x + step_x;
0256 
0257       // Evaluate the constraints at the new point.
0258       Matrix new_Bx(Nw, Nc);
0259       Matrix new_By(Np, Nc);
0260       Row_Vector new_F(Nc);
0261       if (!constraint_calculator.eval(new_x, y, new_F, new_Bx, new_By))
0262         return false;
0263 
0264       // Calculate what we expect the constraints to be at this point,
0265       // given the user's gradients.
0266       Row_Vector test_F = F + step_x.T() * Bx;
0267 
0268       // Check the results.
0269       for (int j = 1; j <= Nc; j++) {
0270         if (test_different(test_F(j), new_F(j), F(j), _args.test_eps())) {
0271           cout << "bad gradient x " << i << " " << j << "\n";
0272           cout << x;
0273           cout << y;
0274           cout << new_x;
0275           cout << F;
0276           cout << new_F;
0277           cout << Bx;
0278           cout << (test_F - new_F);
0279           abort();
0280         }
0281       }
0282     }
0283 
0284     // Numerically check By.
0285     for (int i = 1; i <= Np; i++) {
0286       // Step a little along variable I.
0287       Column_Vector step_y(Np, 0);
0288       step_y(i) = _args.test_step();
0289       Column_Vector new_y = y + step_y;
0290 
0291       // Evaluate the constraints at the new point.
0292       Matrix new_Bx(Nw, Nc);
0293       Matrix new_By(Np, Nc);
0294       Row_Vector new_F(Nc);
0295       if (!constraint_calculator.eval(x, new_y, new_F, new_Bx, new_By))
0296         return false;
0297 
0298       // Calculate what we expect the constraints to be at this point,
0299       // given the user's gradients.
0300       Row_Vector test_F = F + step_y.T() * By;
0301 
0302       // Check the results.
0303       for (int j = 1; j <= Nc; j++) {
0304         if (test_different(test_F(j), new_F(j), F(j), _args.test_eps())) {
0305           cout << "bad gradient y " << i << " " << j << "\n";
0306           cout << x;
0307           cout << y;
0308           cout << new_y;
0309           cout << F;
0310           cout << new_F;
0311           cout << Bx;
0312           cout << (test_F - new_F);
0313           abort();
0314         }
0315       }
0316     }
0317 
0318     // Done!
0319     return true;
0320   }
0321 
0322 }  // namespace hitfit