Base_Constrainer.cc

Go to the documentation of this file.

//
//
// File: src/Base_Constrainer.cc
// Purpose: Abstract base for the chisq fitter classes.
//          This allows for different algorithms to be used.
// Created: Jul, 2000, sss, based on run 1 mass analysis code.
//
// CMSSW File      : src/Base_Constrainer.cc
// Original Author : Scott Stuart Snyder <snyder@bnl.gov> for D0
// Imported to CMSSW by Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>
//

#include "TopQuarkAnalysis/TopHitFit/interface/Base_Constrainer.h"
#include "TopQuarkAnalysis/TopHitFit/interface/matutil.h"
#include "TopQuarkAnalysis/TopHitFit/interface/Defaults.h"
#include <iostream>
#include <cmath>
#include <cstdlib>

using std::abort;
using std::abs;
using std::cout;
using std::ostream;

//*************************************************************************
// Helper function for doing gradient testing.
//

namespace {

  bool test_different(double a, double b, double c, double eps)
  //
  // Purpose: Test if A is significantly different than B.
  //          C sets the scale for the comparison; EPS gives
  //          by how much they may differ.
  //
  {
    double scale = eps * (abs(a) + abs(b) + abs(c));
    if (abs(a) != 0 && abs(b) / abs(a) < 0.1 && abs(a) > scale)
      scale = abs(a) * .5;
    if (scale == 0)
      return false;
    if (scale < eps)
      scale = eps;
    return abs(a - b) > scale;
  }

}  // unnamed namespace

//*************************************************************************

namespace hitfit {

  //*************************************************************************

  Base_Constrainer_Args::Base_Constrainer_Args(const Defaults& defs)
      //
      // Purpose: Constructor.
      //
      // Inputs:
      //   defs -        The Defaults instance from which to initialize.
      //
      : _test_gradient(defs.get_bool("test_gradient")),
        _test_step(defs.get_float("test_step")),
        _test_eps(defs.get_float("test_eps")) {}

  bool Base_Constrainer_Args::test_gradient() const
  //
  // Purpose: Return the test_gradient parameter.
  //          See the header for documentation.
  //
  {
    return _test_gradient;
  }

  double Base_Constrainer_Args::test_step() const
  //
  // Purpose: Return the test_step parameter.
  //          See the header for documentation.
  //
  {
    return _test_step;
  }

  double Base_Constrainer_Args::test_eps() const
  //
  // Purpose: Return the test_eps parameter.
  //          See the header for documentation.
  //
  {
    return _test_eps;
  }

  //*************************************************************************

  Constraint_Calculator::Constraint_Calculator(int nconstraints)
      //
      // Purpose: Constructor.
      //
      // Inputs:
      //   nconstraints- The number of constraint functions.
      //
      : _nconstraints(nconstraints) {}

  int Constraint_Calculator::nconstraints() const
  //
  // Purpose: Return the number of constraint functions.
  //
  // Returns:
  //   The number of constraint functions.
  //
  {
    return _nconstraints;
  }

  //*************************************************************************

  Base_Constrainer::Base_Constrainer(const Base_Constrainer_Args& args)
      //
      // Purpose: Constructor.
      //
      // Inputs:
      //   args -        The parameter settings for this instance.
      //
      : _args(args) {}

  std::ostream& Base_Constrainer::print(std::ostream& s) const
  //
  // Purpose: Print our state.
  //
  // Inputs:
  //   s -           The stream to which to write.
  //
  // Returns:
  //   The stream S.
  //
  {
    s << "Base_Constrainer parameters:\n";
    s << " test_gradient: " << _args.test_gradient() << " test_step: " << _args.test_step()
      << " test_eps: " << _args.test_eps() << "\n";
    return s;
  }

  std::ostream& operator<<(std::ostream& s, const Base_Constrainer& f)
  //
  // Purpose: Print our state.
  //
  // Inputs:
  //   s -           The stream to which to write.
  //   f -           The instance to dump.
  //
  // Returns:
  //   The stream S.
  //
  {
    return f.print(s);
  }

  bool Base_Constrainer::call_constraint_fcn(Constraint_Calculator& constraint_calculator,
                                             const Column_Vector& x,
                                             const Column_Vector& y,
                                             Row_Vector& F,
                                             Matrix& Bx,
                                             Matrix& By) const
  //
  // Purpose: Call the constraint function for the point x, y.
  //          Return F, Bx, By, and a flag saying if the
  //          point is acceptable.
  //
  //          If test_gradient is on, we verify the gradients returned
  //          by also computing them numerically.
  //
  // Inputs:
  //   constraints - The user-supplied object to evaluate the constraints.
  //   x(Nw) -       Vector of well-measured quantities where we evaluate
  //                 the constraints.
  //   y(Np) -       Vector of poorly-measured quantities where we evaluate
  //                 the constraints.
  //
  // Outputs:
  //   F(Nc) -       The results of the constraint functions.
  //   Bx(Nw,Nc) -   Gradients of F with respect to x.
  //   By(Np,Nc) -   Gradients of F with respect to y.
  //
  // Returns:
  //   True if the point is accepted, false if it was rejected.
  {
    // Call the user's function.
    bool val = constraint_calculator.eval(x, y, F, Bx, By);

    // If we're not doing gradients numerically, we're done.
    if (!_args.test_gradient())
      return val;

    // Bail if the point was rejected.
    if (!val)
      return false;

    int Nw = x.num_row();
    int Np = y.num_row();
    int Nc = F.num_col();

    // Numerically check Bx.
    for (int i = 1; i <= Nc; i++) {
      // Step a little along variable I.
      Column_Vector step_x(Nw, 0);
      step_x(i) = _args.test_step();
      Column_Vector new_x = x + step_x;

      // Evaluate the constraints at the new point.
      Matrix new_Bx(Nw, Nc);
      Matrix new_By(Np, Nc);
      Row_Vector new_F(Nc);
      if (!constraint_calculator.eval(new_x, y, new_F, new_Bx, new_By))
        return false;

      // Calculate what we expect the constraints to be at this point,
      // given the user's gradients.
      Row_Vector test_F = F + step_x.T() * Bx;

      // Check the results.
      for (int j = 1; j <= Nc; j++) {
        if (test_different(test_F(j), new_F(j), F(j), _args.test_eps())) {
          cout << "bad gradient x " << i << " " << j << "\n";
          cout << x;
          cout << y;
          cout << new_x;
          cout << F;
          cout << new_F;
          cout << Bx;
          cout << (test_F - new_F);
          abort();
        }
      }
    }

    // Numerically check By.
    for (int i = 1; i <= Np; i++) {
      // Step a little along variable I.
      Column_Vector step_y(Np, 0);
      step_y(i) = _args.test_step();
      Column_Vector new_y = y + step_y;

      // Evaluate the constraints at the new point.
      Matrix new_Bx(Nw, Nc);
      Matrix new_By(Np, Nc);
      Row_Vector new_F(Nc);
      if (!constraint_calculator.eval(x, new_y, new_F, new_Bx, new_By))
        return false;

      // Calculate what we expect the constraints to be at this point,
      // given the user's gradients.
      Row_Vector test_F = F + step_y.T() * By;

      // Check the results.
      for (int j = 1; j <= Nc; j++) {
        if (test_different(test_F(j), new_F(j), F(j), _args.test_eps())) {
          cout << "bad gradient y " << i << " " << j << "\n";
          cout << x;
          cout << y;
          cout << new_y;
          cout << F;
          cout << new_F;
          cout << Bx;
          cout << (test_F - new_F);
          abort();
        }
      }
    }

    // Done!
    return true;
  }

}  // namespace hitfit