FunctionController.cxx

Go to the documentation of this file.

#ifdef HAVE_CONFIG_H
#include "config.h"
#else
#ifdef _MSC_VER
#include "msdevstudio/MSconfig.h"
#endif
#endif

#include "FunctionController.h"
#include "Gammaq.h"
#include "DisplayController.h"

#include "datasrcs/DataSourceController.h"
#include "datareps/FunctionRep.h"
#include "datasrcs/NTuple.h"
#include "datasrcs/TupleCut.h"

#include "functions/FunctionBase.h"
#include "functions/FunctionFactory.h"

#include "minimizers/Fitter.h"
#include "minimizers/FitterFactory.h"
#include "minimizers/NTupleFCN.h"
#include "minimizers/NumLinAlg.h"

#include "plotters/PlotterBase.h"
#include "plotters/CompositePlotter.h"
#include "projectors/BinningProjector.h"
#include "projectors/FunctionProjector.h"
#include "projectors/NTupleProjector.h"

#include <algorithm>
#include <functional>
#include <iterator>

#ifdef SSTREAM_DEFECT
#include <strstream>
using std::ostrstream;
#else
#include <sstream>
using std::ostringstream;
#endif

#include <cfloat>

using namespace hippodraw;
using namespace hippodraw::Numeric;

#ifdef ITERATOR_MEMBER_DEFECT
using namespace std;
#else
using std::distance;
using std::find;
using std::find_if;
using std::list;
using std::mem_fun;
using std::string;
using std::vector;
using std::sqrt;
using std::cos;
using std::sin;
using std::atan;
using std::min;
using std::max;
using std::abs;
#endif

FunctionController * FunctionController::s_instance = 0;

FunctionController::FunctionController ( )
  : m_plotter ( 0 ), 
    m_x ( 0 ), 
    m_y ( 0 ),
    m_confid_count ( 0 )
{
  m_deltaXSq.resize( 6 );

  m_deltaXSq[ 0 ] = 15.1;
  m_deltaXSq[ 1 ] = 18.4;
  m_deltaXSq[ 2 ] = 21.1;
  m_deltaXSq[ 3 ] = 23.5;
  m_deltaXSq[ 4 ] = 25.7;
  m_deltaXSq[ 5 ] = 27.8;
}

FunctionController::~FunctionController ( )
{
}

FunctionController * FunctionController::instance ( )
{
  if ( s_instance == 0 ) {
    s_instance = new FunctionController ( );
  }
  return s_instance;
}

const vector < string > &
FunctionController::
getFitterNames () const
{
  FitterFactory * factory = FitterFactory::instance ();

  return factory -> names ();
}

int FunctionController::
getUniqueNonFunctionIndex ( const PlotterBase * plotter ) const
{
  int index = -1;
  int count = 0;

  int number = plotter->getNumDataReps ();
  for ( int i = 0; i < number; i++ ) {
    DataRep * r = plotter->getDataRep ( i );
    FunctionRep * fr = dynamic_cast < FunctionRep * > ( r );
    if ( fr != 0 ) continue;
    index = i;
    count++;
  }
  if ( count != 1 ) index = -1;

  return index;
}

void
FunctionController::
findFunctions ( const PlotterBase * plotter ) const
{
  m_func_reps.clear ();
  int number = plotter->getNumDataReps ();

  for ( int i = 0; i < number; i++ ) {
    DataRep * rep = plotter->getDataRep ( i );
    FunctionRep * frep = dynamic_cast < FunctionRep * > ( rep );
    if ( frep == 0 ) continue;
    m_func_reps.push_back ( frep );
  }
}

FunctionRep * 
FunctionController::
createFunctionRep ( const std::string & name,
                    DataRep * rep )
{
  FunctionRep * frep = new FunctionRep ( name, rep );

  return frep;
}

FunctionRep * 
FunctionController::
createFunctionRep ( FunctionBase * function, 
                    DataRep * rep )
{
  FunctionRep * frep = new FunctionRep ( function, rep );

  const vector < string > & fitters = getFitterNames ();
  setFitter ( frep, fitters[0] );

  return frep;
}

FunctionBase *
FunctionController::
addFunction ( PlotterBase * plotter, const std::string & name )
{
  DataRep * rep = plotter -> getDataRep ( 0 );
  addFunction ( plotter, name, rep );
  FunctionRep * frep = getFunctionRep ( plotter );

  return frep -> getFunction ();
}

int FunctionController::addFunction ( PlotterBase * plotter,
                                      const std::string & name, 
                                      DataRep * rep )
{
  FunctionRep * func_rep = createFunctionRep ( name, rep );
  func_rep -> initializeWith ( rep );
  return addFunctionRep ( plotter, rep, func_rep );
}

int FunctionController::addFunction ( PlotterBase * plotter, 
                                      FunctionRep * func_rep )
{
  int index = plotter->activePlotIndex ();

  if ( index < 0 ) {
    index = getUniqueNonFunctionIndex ( plotter );
  }
  if ( index < 0 ) return -1;
  plotter->setActivePlot ( index, false );
  DataRep * drep = plotter->getDataRep ( index );
  func_rep -> initializeWith ( drep );
  return addFunctionRep ( plotter, drep, func_rep );
}


void
FunctionController::
addDataRep ( PlotterBase * plotter, DataRep * rep )
{
  FunctionRep * frep = dynamic_cast < FunctionRep * > ( rep );
  if ( frep == 0 ) { // not a function rep
    DisplayController * controller = DisplayController::instance ();
    controller -> addDataRep ( plotter, rep );
  }
  else { // a function rep
    if ( frep -> isComposite () == false ) {
      DataRep * target = frep -> getTarget ();
      addFunctionRep ( plotter, target, frep );
    }
  }
}

int
FunctionController::
addFunctionRep ( PlotterBase * plotter,
                 DataRep * drep,
                 FunctionRep * func_rep )
{
  plotter->addDataRep ( func_rep );
  findFunctions ( plotter ); // to re-cache

  FunctionBase * function = func_rep->getFunction ();
  if ( function->isComposite () ) {
    buildComposite ( function );
  }
  func_rep->notifyObservers ();
  const vector < string > & fitters = getFitterNames ();
  setFitter ( func_rep, fitters[0] );

  // Check the list of functions. If there are 2 or more functionReps, and
  // all are non LinearSum, and the function you added right now is not a
  // linearSum, then add a linearSum functionrep.
  
  // Allow at most one Linear Sum in the list of all functionreps on a
  // plotter. If a Linear sum is already present, then add the current
  // function to that composite.

  if ( function->isComposite () == false ) {
    int count = 0;
    vector< FunctionRep * >:: const_iterator it = m_func_reps.begin (); 
    for ( ; it != m_func_reps.end (); ++it ) {
      FunctionRep * localfp = *it;
      FunctionBase * localfunction = localfp->getFunction ();
      const string & name = localfunction->name ();

      if ( name != "Linear Sum" ){
        count ++;
      }
      else{
        // The current function is a Linear Sum. So add the new function
        // to its list of functions.
        localfunction->addToComposite ( function );
        count = 0;
        break;
      }
    }
    
    if ( count > 1 ){
      addFunction ( plotter, "Linear Sum", drep );
    }
  
  }

  int retval = findIndex ( func_rep );
  return retval;

}

void FunctionController::buildComposite ( const PlotterBase * plotter )
{
  FunctionRep * composite_rep = getComposite ( plotter );
  if ( composite_rep == 0 ) return;

  FunctionBase * composite = composite_rep->getFunction ();
  buildComposite ( composite );
}

void FunctionController::buildComposite ( FunctionBase * composite )
{
  vector< FunctionRep * > :: const_iterator it = m_func_reps.begin ();

  for ( ; it != m_func_reps.end (); ++it ) {
    FunctionBase * function = (*it)->getFunction ();
    if ( function != composite ) {
      composite->addToComposite ( function );
    }
  }
}

void FunctionController::removeFunction ( PlotterBase * plotter,
                                          unsigned int index )
{
  findFunctions ( plotter );
  assert ( index < m_func_reps.size () );

  vector < FunctionRep * > :: iterator it = m_func_reps.begin () + index;
  FunctionRep * frep = *it;
  m_func_reps.erase ( it );

  plotter->removeDataRep ( frep );

  FunctionBase * function = frep -> getFunction();  

  for ( unsigned int i = 0; i < m_func_reps.size(); i++ ){
    FunctionBase * curFunction = m_func_reps[i]->getFunction();
    if ( curFunction->isComposite() ){
      curFunction->removeFromComposite ( function );
      m_func_reps[i]->setDirty ( true );
    }
  }

  // Remove all composites who now have one or none functions within them.
  
  vector < FunctionRep * >::iterator iter = m_func_reps.begin ();
  
  // This while loop is tricky because modifying the vector that we
  // are iterating over.   
  while ( iter != m_func_reps.end() ){
    FunctionBase * func = ( *iter )->getFunction();
    if ( func->isComposite() ) {
      if ( func->count() < 2 ) {
        plotter->removeDataRep ( *iter );
        delete *iter;
        iter = m_func_reps.erase ( iter ); // updates iter 
      }
      else {
        iter ++;
      }
    }
    else {
      iter ++;
    }
  }
  delete frep; // the functionrep
}

int FunctionController::findIndex ( const FunctionRep * fp )
{
  vector< FunctionRep * >::iterator it
    = find ( m_func_reps.begin (), m_func_reps.end (), fp );

  if ( it == m_func_reps.end () ) return -1;

#ifdef DISTANCE_DEFECT
  return it - m_func_reps.begin ();
#else
  return distance ( m_func_reps.begin (), it );
#endif
}

bool
FunctionController::
hasFunction ( const PlotterBase * plotter, const DataRep * rep )
{
  bool yes = false;
  assert ( plotter != 0 );

  findFunctions ( plotter );

  if ( m_func_reps.empty () == false ) {
    if ( rep != 0 ) {
      vector< FunctionRep * >:: const_iterator it = m_func_reps.begin ();

      for ( ; it != m_func_reps.end (); ++it ) {
        FunctionRep * fp = *it;
        if ( fp->getTarget() == rep ) {
          yes = true;
          break;
        }
      }
    } 
    else {
      yes = true;
    }
  }

  return yes;
}

const vector< string > & 
FunctionController::
functionNames ( PlotterBase * plotter, DataRep * rep )
{
  findFunctions ( plotter );

  m_f_names.clear ();
  vector< FunctionRep * >:: const_iterator it = m_func_reps.begin ();

  for ( ; it != m_func_reps.end (); ++it ) {
    FunctionRep * fp = *it;
    if ( rep != 0 && fp->getTarget() != rep ) continue;
    FunctionBase * function = fp->getFunction ();
    const string & name = function->name ();
    m_f_names.push_back ( name );
  }

  return m_f_names;
}

const std::vector < FunctionRep * > & 
FunctionController::
getFunctionReps ( const PlotterBase * plotter ) const
{
  FunctionController * tmp = const_cast< FunctionController * > ( this );
  tmp->findFunctions ( plotter );
  
  return m_func_reps;
}

FunctionRep *
FunctionController::
getComposite ( const PlotterBase * plotter ) const
{
  FunctionController * tmp = const_cast< FunctionController * > ( this );
  tmp->findFunctions ( plotter );

#ifdef MEMFUN_DEFECT
  vector < FunctionRep * >::iterator first = m_func_reps.begin();
  while ( first != m_func_reps.end() ) {
          FunctionRep * frep = *first++;
          if ( frep->isComposite () ) break;
  }
#else
  vector < FunctionRep * >::const_iterator first
    = find_if ( m_func_reps.begin(), m_func_reps.end(),
                mem_fun ( &FunctionRep::isComposite ) );
#endif

  if ( first == m_func_reps.end() ) return 0;

  return *first;
}

const vector < double > &
FunctionController::
getErrors ( const PlotterBase * plotter )
{
  FunctionRep * rep = getFunctionRep ( plotter );

  return rep -> principleErrors ();
}

void
FunctionController::
setErrorsFromComposite ( const PlotterBase * plotter,
                         const FunctionRep * composite )
{
  const vector < double > & errors = composite -> principleErrors();
  vector < double >::const_iterator begin = errors.begin();

  getFunctionReps ( plotter );
  vector < FunctionRep * >::iterator first = m_func_reps.begin();

  while ( first != m_func_reps.end() ) {
    FunctionRep * rep = *first++;
    if ( rep -> isComposite() ) continue;
    ProjectorBase * pbase = rep -> getProjector ();
    FunctionProjector * projector 
      = dynamic_cast < FunctionProjector * > ( pbase );
    assert ( projector );
    const vector < double > & t = projector -> principleErrors ();
    unsigned int number = t.size();
    projector->setPrincipleErrors ( begin, 
                                    begin + number );
    begin += number;
  }

}

bool
FunctionController::
fitFunction ( PlotterBase * plotter, unsigned int )
{
  FunctionRep * rep = getFunctionRep ( plotter );

  return fitFunction ( plotter, rep );
}

bool
FunctionController::
fitFunction ( PlotterBase * plotter, FunctionRep * rep )
{
  FunctionRep * composite = getComposite ( plotter );
  bool ok = false;
  if ( composite != 0 ) {
    ok = composite -> fitFunction ();
    setErrorsFromComposite ( plotter, composite );
  }
  else {
    ok = rep -> fitFunction ();
  }

  return ok;
}

void
FunctionController::
tryFitFunction ( PlotterBase * plotter, const DataRep * datarep )
{
  saveParameters ( plotter );
  double chi_sq = getObjectiveValue ( plotter, datarep );
  FunctionRep * func_rep = getFunctionRep ( plotter, datarep );
  bool ok = fitFunction ( plotter, func_rep );
  if ( ok ) {
    double new_chi_sq = getObjectiveValue ( plotter, datarep );
    if ( new_chi_sq > chi_sq ) {
      restoreParameters ( plotter );
    }
  } else {
    restoreParameters ( plotter );
  }
}

void
FunctionController::
fitFunctions ( PlotterBase * plotter )
{
  if ( hasFunction ( plotter, 0 ) ) {
    buildComposite ( plotter );
    FunctionRep * frep = getFunctionRep ( plotter );
    fitFunction ( plotter, frep );
  }}

void
FunctionController::
setFitRange ( PlotterBase * plotter, const Range & range )
{
  if ( hasFunction ( plotter, 0 ) ) {
    FunctionRep * frep = getFunctionRep ( plotter );
    frep -> setCutRange ( range );
    plotter -> update ();
  }
}

void
FunctionController::
setFitRange ( PlotterBase * plotter, double low, double high )
{
  const Range range ( low, high );
  setFitRange ( plotter, range );
}

const std::vector < double >
FunctionController::
parameters ( const PlotterBase * plotter, unsigned int )
{
  FunctionRep * rep = getFunctionRep ( plotter );

  return rep -> parameters ();
}

const std::vector < double >
FunctionController::
principleErrors ( const PlotterBase * plotter, unsigned int )
{
  FunctionRep * rep = getFunctionRep ( plotter );

  return rep -> principleErrors ();
}

void FunctionController::saveParameters ( PlotterBase * plotter )
{
  findFunctions ( plotter );
  vector< FunctionRep * >::iterator it = m_func_reps.begin ();

  for ( ; it != m_func_reps.end (); ++it ) {
    FunctionRep * frep = *it;
    FunctionBase * function = frep->getFunction ();
    if ( function->isComposite () ) continue;
    frep->saveParameters ();
  }
}

void FunctionController::restoreParameters ( PlotterBase * plotter )
{
  findFunctions ( plotter );
  vector< FunctionRep * >::iterator it = m_func_reps.begin ();

  for ( ; it != m_func_reps.end (); ++it ) {
    FunctionRep * frep = *it;
    FunctionBase * function = frep->getFunction ();
    if ( function->isComposite () ) continue;
    frep->restoreParameters ();
  }

}

FunctionRep *
FunctionController::
getFunctionRep ( const PlotterBase * plotter ) const
{
  FunctionRep * frep = 0;
  DataRep * rep = 0;
  int index = plotter->activePlotIndex ();

  if ( index >= 0 ) {
    rep = plotter -> getDataRep ( index );
  }
  else {
    rep = plotter -> getTarget ();
  }

  FunctionRep * test = dynamic_cast < FunctionRep * > ( rep );
  if ( test != 0 ) { // if ctrl-clicked, could be the function
    frep = test;
  }
  else {
    frep = getFunctionRep ( plotter, rep );
  }

  return frep;
}

FunctionRep *
FunctionController::
getFunctionRep ( const PlotterBase * plotter, const DataRep * datarep ) const
{
  FunctionRep * frep = 0;

  findFunctions ( plotter ); 

  for ( unsigned int i = 0; i < m_func_reps.size(); i++ ) {
    FunctionRep * rep = m_func_reps[i];
    const DataRep * drep = rep -> getTarget ();
    if ( drep != datarep ) continue;
    frep = rep;
    if ( frep -> isComposite () ) break; // use composite if found
  }

  return frep;
}

ViewBase * FunctionController::
createFuncView ( const ViewFactory * factory,
                 PlotterBase * plotter,
                 const std::string & name )
{
  FunctionRep * frep = getFunctionRep ( plotter );
  assert ( frep != 0 );

  DisplayController * controller = DisplayController::instance ();
  string nullstring ("");

  return controller->createTextView ( factory, frep, name, nullstring );
}

bool
FunctionController::
setFitter ( const PlotterBase * plotter, const std::string & name )
{
  FunctionRep * frep = getFunctionRep ( plotter );
  bool ok = false;
  if ( frep != 0 ) {
    ok = setFitter ( frep, name );
  }

  return ok;
}

const string &
FunctionController::
getFitterName ( const PlotterBase * plotter )
{
  FunctionRep * frep = getFunctionRep ( plotter );
  Fitter * fitter = frep -> getFitter ();

  return fitter -> name ();
}

Fitter *
FunctionController::
getFitter ( const PlotterBase * plotter )
{
  FunctionRep * frep = getFunctionRep ( plotter );
  Fitter * fitter = frep -> getFitter ();

  return fitter;
}


bool
FunctionController::
setFitter ( FunctionRep * frep, const std::string & name )
{
  FitterFactory * factory = FitterFactory::instance ();
  Fitter * fitter = factory -> create ( name );

  return setFitter ( frep, fitter );
}

bool
FunctionController::
setFitter ( FunctionRep * frep, Fitter * fitter )
{
  DataRep * drep = frep -> getTarget ();
  FunctionBase * function = frep -> getFunction ();
  StatedFCN * fcn = fitter -> getFCN ();

  if ( fcn -> needsIntegrated () == true ) {
    bool yes = drep -> isAxisBinned ( Axes::X );
    if ( ! yes ) return false;
  }

  fcn -> setFunction ( function );

  NTupleFCN * ntfcn = dynamic_cast < NTupleFCN * > ( fcn );
  assert ( ntfcn != 0 );

  const DataSource * target = drep -> getProjectedValues ();
  ntfcn -> setDataSource ( target );
  ntfcn -> setUseErrors ( );

  frep -> setFitter ( fitter );

  return true;
}

double
FunctionController::
getObjectiveValue ( const PlotterBase * plotter, const DataRep * datarep )
{
  FunctionRep * rep = getFunctionRep ( plotter, datarep );

  ProjectorBase * pbase = rep->getProjector ( );
  FunctionProjector * projector 
    = dynamic_cast < FunctionProjector * > ( pbase );

  return projector -> objectiveValue ();
}

const vector < vector < double > > &
FunctionController::
getCovarianceMatrix ( const PlotterBase * plotter )
{
  FunctionRep * rep = getFunctionRep ( plotter );

  ProjectorBase * pbase = rep->getProjector ( );
  FunctionProjector * projector 
    = dynamic_cast < FunctionProjector * > ( pbase );

  return projector -> covariance ();
}

double
FunctionController::
getChiSquared ( const PlotterBase * plotter )
{
  const DataRep * rep = plotter -> getTarget ();

  return getObjectiveValue ( plotter, rep );
}

int
FunctionController::
getDegreesOfFreedom ( const PlotterBase * plotter )
{
  FunctionRep * rep = getFunctionRep ( plotter );

  ProjectorBase * pbase = rep->getProjector ( );
  FunctionProjector * projector 
    = dynamic_cast < FunctionProjector * > ( pbase );

  return projector -> degreesOfFreedom ();
}

NTuple *
FunctionController::
createNTuple ( PlotterBase * plotter )
{
  int old_index = plotter -> activePlotIndex ();

  int index = getUniqueNonFunctionIndex ( plotter );
  if ( index < 0 ) return 0;

  int saved_index = plotter -> activePlotIndex ();
  plotter -> setActivePlot ( index, false ); // no redraw
  NTuple * ntuple = plotter -> createNTuple ();
  plotter -> setActivePlot ( saved_index, false );

  FunctionRep * rep = getFunctionRep ( plotter );
  if ( rep == 0 ) return ntuple;

  ProjectorBase * pbase = rep->getProjector ( );
  FunctionProjector * projector 
    = dynamic_cast < FunctionProjector * > ( pbase );
  FunctionBase * function = projector -> function();

  vector < double > & x = ntuple -> getColumn ( 0 ); // X coordinate
  vector < double > & y = ntuple -> getColumn ( 1 ); // Y coordingate
  unsigned int size = ntuple -> rows ();
  vector < double > values ( size );
  vector < double > residuals ( size );
  
  for ( unsigned int i = 0; i < x.size(); i++ ) {
    values [i] = function -> operator () ( x[i] );
    residuals [i] =  y[i] - values [i];
  }

  ntuple -> addColumn ( function->name (), values );
  ntuple -> addColumn ( "Residuals", residuals );

  plotter -> setActivePlot ( old_index, true );

  return ntuple;
}
 
PlotterBase *
FunctionController::
createResidualsDisplay ( PlotterBase * plotter )
{
  NTuple * ntuple = createNTuple ( plotter );
  ntuple -> setTitle ( plotter -> getTitle () );
  DataSourceController::instance () -> registerNTuple ( ntuple );
  vector < string > bindings ( 4 );
  bindings[0] = ntuple -> getLabelAt ( 0 );
  bindings[1] = "Residuals";
  bindings[2] = "nil";
  bindings[3] = ntuple -> getLabelAt ( 3 );

  DisplayController * controller = DisplayController::instance ();

// Make scaling of x-axis match that of the original plot.
  PlotterBase * new_plotter = controller -> createDisplay ( "XY Plot",
                                                            *ntuple,
                                                            bindings );
  controller->setLog( new_plotter, "x", 
                      controller -> getLog( plotter, "x" ) );
  
  return new_plotter;
}

PlotterBase *
FunctionController::
createNewEllipsoidDisplay ( PlotterBase * masterPlot )
{
  int n = 100;
  double xmin, xmax, ymin, ymax;
  
  // Set up the NTuple
  NTuple* ntuple = new( NTuple );
  ellipsoidNTuple ( masterPlot, ntuple, n, xmin, xmax, ymin, ymax );
  ntuple -> setTitle ( masterPlot -> getTitle () );
  
  vector < string > bindings ( 1 );
  bindings[0] = ntuple ->  getLabelAt ( 0 );
  
  DisplayController * dcontroller = DisplayController::instance ();

  // Z plot. X and Y coordinates are merely inidices
  PlotterBase * ellipsePlot = dcontroller ->
    createDisplay ( "Z Plot", *ntuple, bindings );
  
  // Rescale and translate the Z plot so that  the plot
  // becomes contained in the rectangle bound by xmin, xmax
  // ymin and ymax.
  ellipsePlot -> setBinWidth( "x", ( xmax - xmin ) / ( n - 1.0 ) );
  ellipsePlot -> setBinWidth( "y", ( ymax - ymin ) / ( n - 1.0 ) );
  
  ellipsePlot -> setOffsets ( xmin, ymin );
  
  ellipsePlot -> setRange ( string( "X" ), xmin, xmax );
  ellipsePlot -> setRange ( string( "Y" ), ymin, ymax );
  
  // Establish the relation between source masterPlot and
  // this new plot. This relationship shall be exploited
  // when we refresh the error plots.
  int index = dcontroller -> activeDataRepIndex( masterPlot );
  ellipsePlot -> setParentDataRepIndex( index );
  ellipsePlot -> setParentPlotter( masterPlot );
  
  return ellipsePlot;
}


void
FunctionController::
ellipsoidNTuple ( PlotterBase* masterPlot,
                  NTuple* nt, int n,
                  double & xmin, double & xmax,
                  double & ymin, double & ymax )
{
#ifdef SSTREAM_DEFECT
  ostrstream strm_text;
#else
  ostringstream strm_text;
#endif
  strm_text << "Confidence ellipse ";
  m_confid_count++;
  strm_text << m_confid_count;
  const string str_text ( strm_text.str() );
  nt -> setName ( str_text );

  DataSourceController * controller = DataSourceController::instance ();
  controller -> registerNTuple ( nt );

  // Function Rep associated with the selected datarep
  FunctionRep * frep = getFunctionRep ( masterPlot );
  assert( frep );
  
  ProjectorBase * pbase = frep -> getProjector ( );
  FunctionProjector * fprojector 
    = dynamic_cast < FunctionProjector * > ( pbase );
  
  // Get projected covarience i.e. take the submatrix
  // out of the covarience matrix which corrosponds
  // to the 2 parameters of intrest whose correlation
  // we would like to see.
  vector< vector< double > > Sigma( 2 );
  Sigma[0].resize( 2, 0.0 );
  Sigma[1].resize( 2, 0.0 );

  const vector < vector < double > > & covariance 
    = fprojector -> covariance ();  

  Sigma[0][0] = covariance [m_x][m_x];
  Sigma[0][1] = covariance [m_x][m_y];
  Sigma[1][0] = covariance [m_y][m_x];
  Sigma[1][1] = covariance [m_y][m_y];

  // Invert the projected covarience 
  vector< vector< double > > SigmaInv;
  invertMatrix( Sigma, SigmaInv );
  
  // Decide the center of the ellipse to be parameter
  // value of the 2 parameters of intrest whose correlation
  // we would like to see.
  vector< double > xbar( 2 );
  
  const Fitter * fitter = getFitter ( masterPlot );
  vector < double > free_parms;
  fitter -> fillFreeParameters ( free_parms );
  xbar[ 0 ] = free_parms [ m_x ];
  xbar[ 1 ] = free_parms [ m_y ];
  
  // Get the bounding ellipse and its bounds. Bounding ellipse is the
  // 99.99% confidence ellipsoid. For mu = 2 the delta chi-square turns from 
  // Numerical Recipies in C as 18.4. etc etc.
  unsigned int mu = free_parms.size();
  NTuple * boundingTuple = ellipse( xbar, SigmaInv, m_deltaXSq[ mu - 1 ] ); 
  
  xmin = boundingTuple -> minElement ( 0 );
  xmax = boundingTuple -> maxElement ( 0 );
  ymin = boundingTuple -> minElement ( 1 );
  ymax = boundingTuple -> maxElement ( 1 );
  
  delete boundingTuple; // Served its purpose
  
  // Create the ntuple for the contour plot
  // This ntuple is to be build column-wise
  // keeping in view the way Z-plot works.
  vector< double > p( n * n ), a( 2 );
  
  double dx    = ( xmax - xmin ) / ( n - 1.0 );
  double dy    = ( ymax - ymin ) / ( n - 1.0 );
  double delta = 0.0;
  
  for( int i = 0; i < n; i++ )
    for( int j = 0; j < n; j++ )
      {
        a[ 0 ] = xmin + i * dx - xbar[ 0 ] ;
        a[ 1 ] = ymin + j * dy - xbar[ 1 ];
        delta  = quadraticProduct( SigmaInv, a );
        p[ n * i + j ] = 1 - gammq( mu / 2.0, delta / 2.0 ); 
      }

  nt -> clear();
  nt -> addColumn( "Percent", 100 * p );
  
  return;
}

PlotterBase *
FunctionController::
refreshEllipsoidDisplay ( PlotterBase * plot )
{
  int n = 100;
  double xmin, xmax, ymin, ymax;
  
  PlotterBase* masterPlot = plot->getParentPlotter();

  // Set up the NTuple
  NTuple* ntuple = new( NTuple );
  ellipsoidNTuple ( masterPlot, ntuple, n, xmin, xmax, ymin, ymax );
  ntuple -> setTitle ( masterPlot -> getTitle () );
  
  // First get the selected datarep from the plotter. The DataRep 
  // will have a projector that is a NTupleProjector. It doesn't
  // know, but we do, so we downcast and set the ntuple.
  DataRep * drep = plot -> selectedDataRep();
  NTupleProjector * ntProjector =
    dynamic_cast < NTupleProjector * > ( drep -> getProjector() );
  ntProjector -> setNTuple ( ntuple );
  
  NTuple * nt = const_cast < NTuple * > ( ntuple );
  nt  -> addObserver ( ntProjector );
    
  // Rescale and translate the Z plot so that  the plot
  // becomes contained in the rectangle bound by xmin, xmax
  // ymin and ymax.
  plot -> setBinWidth( "x", ( xmax - xmin ) / ( n - 1.0 ) );
  plot -> setBinWidth( "y", ( ymax - ymin ) / ( n - 1.0 ) );
  
  plot -> setOffsets ( xmin, ymin );
  
  plot -> setRange ( string( "X" ), xmin, xmax );
  plot -> setRange ( string( "Y" ), ymin, ymax );
    
  return plot;
}

NTuple * 
FunctionController::
ellipse ( const std::vector< double > & xbar,
          std::vector< std::vector < double > > & SigmaInv,
          double Csquare )
{
  // First agrument should be a 2 D vector, the center of the ellipse //
  assert( xbar.size() == 2 );

  // Second argument should be a 2 x 2 SPD matrix i.e. two rows two cols //
  assert( SigmaInv.size() == 2 );
  assert( SigmaInv[0].size() == 2 || SigmaInv[1].size() == 2 );
  
  // Second argument should be a 2 x 2 -- Symmetric --  PD matrix
  // Under numerical roundoffs it might not be true so we take mean of the
  // enteries.
  double temp = ( SigmaInv[0][1] + SigmaInv[1][0] ) / 2;
  SigmaInv[0][1] = temp;
  SigmaInv[1][0] = temp;
  
  // Third argument should be a positive scalar "c^2", Square of "Radius" // 
  assert( Csquare > DBL_EPSILON );
  
  // The eigenvalues of the SigmaInv matrix
  double b = -( SigmaInv[0][0] + SigmaInv[1][1] );
  double c = SigmaInv[0][0] * SigmaInv[1][1] - SigmaInv[0][1] * SigmaInv[1][0];

  double lambda1 = ( -b + sqrt( b * b - 4 * c ) ) / 2;
  double lambda2 = ( -b - sqrt( b * b - 4 * c ) ) / 2;
    
  // Determinig the angles the ellipse axis make with X-axis
  double alpha1 =  atan( (SigmaInv[0][0] - lambda1) / SigmaInv[0][1] );
  double a      =  sqrt( Csquare / lambda1 );
  b             =  sqrt( Csquare / lambda2 );

  // Creating a rotation matrix
  double Rot00 =  cos( alpha1 );
  double Rot11 =  cos( alpha1 );
  double Rot01 = -sin( alpha1 );
  double Rot10 =  sin( alpha1 );

  // Generating N-points on the ellipse
  int N = 100;
  vector< double > x, y;
  
  x.resize( N, 0.0 );
  y.resize( N, 0.0 );
  
  for( int i = 0; i < N; i++)
    {
      // Unrotated untranslated point
      double theta = ( 2 * M_PI * i ) / ( (double) (N-1) );
      double x0 = a * cos( theta );
      double x1 = b * sin( theta );
      
      //x = Rot * x; i.e. force in a rotation
      double xrot0 = Rot00 * x0 + Rot01 * x1;
      double xrot1 = Rot10 * x0 + Rot11 * x1;
      
      //x = x + xbar; i.e. force in a translation
      x[i] = xrot0 + xbar[0];
      y[i] = xrot1 + xbar[1];
    }
  
  // Create an Ntuple out of the above vector
  NTuple* ntuple = new( NTuple );
  ntuple -> addColumn( "X", x );
  ntuple -> addColumn( "Y", y );
  
  return ntuple;
}

int FunctionController::setEllpsoidParamIndex( hippodraw::Axes::Type axes,
                                               int index )
{
  assert( axes == Axes::X || axes == Axes::Y );

  if( axes == Axes::X )
    m_x = index;
  
  if( axes == Axes::Y )
    m_y = index;
      
  return EXIT_SUCCESS;
}

bool
FunctionController::
functionExists ( const std::string & name )
{
  FunctionFactory * factory = FunctionFactory::instance ();

  return factory -> exists ( name );
}

bool
FunctionController::
isCompatible ( const std::string & function,
               const std::string & fitter )
{
  bool yes = true;

  FunctionFactory * fun_fac = FunctionFactory::instance ();
  FunctionBase * fun_proto = fun_fac -> prototype ( function );

  FitterFactory * fit_fac = FitterFactory::instance ();
  Fitter * fit_proto = fit_fac -> prototype ( fitter );

  if ( fit_proto -> needsDerivatives () == true  &&
       fun_proto -> hasDerivatives () == false ) {
    yes = false;
  }

  return yes;
}

---