FunctionController.cxx

Go to the documentation of this file.
00001 
00012 #ifdef _MSC_VER
00013 #include "msdevstudio/MSconfig.h"
00014 #endif
00015 
00016 #include "FunctionController.h"
00017 
00018 #include "Gammaq.h"
00019 #include "DisplayController.h"
00020 
00021 #include "datareps/CompositeFunctionRep.h"
00022 #include "datareps/FunctionParameter.h"
00023 #include "datareps/FunctionRep1.h"
00024 #include "datareps/FunctionRep2.h"
00025 
00026 #include "datasrcs/DataSourceController.h"
00027 #include "datasrcs/NTuple.h"
00028 #include "datasrcs/TupleCut.h"
00029 
00030 #include "functions/FunctionBase.h"
00031 #include "functions/FunctionFactory.h"
00032 
00033 #include "minimizers/Fitter.h"
00034 #include "minimizers/FitterFactory.h"
00035 #include "minimizers/NumLinAlg.h"
00036 
00037 #include "pattern/string_convert.h"
00038 #include "plotters/PlotterBase.h"
00039 #include "plotters/CompositePlotter.h"
00040 #include "projectors/BinningProjector.h"
00041 #include "projectors/NTupleProjector.h"
00042 
00043 #include <algorithm>
00044 #include <functional>
00045 #include <iterator>
00046 #include <stdexcept>
00047 
00048 #include <cmath>
00049 #include <cassert>
00050 
00051 #ifdef ITERATOR_MEMBER_DEFECT
00052 using namespace std;
00053 #else
00054 using std::distance;
00055 using std::find;
00056 using std::find_if;
00057 using std::list;
00058 using std::mem_fun;
00059 using std::string;
00060 using std::vector;
00061 using std::sqrt;
00062 using std::cos;
00063 using std::sin;
00064 using std::atan;
00065 using std::min;
00066 using std::max;
00067 using std::abs;
00068 #endif
00069 
00070 using namespace hippodraw;
00071 
00072 using namespace Numeric;
00073 
00074 FunctionController * FunctionController::s_instance = 0;
00075 
00076 FunctionController::FunctionController ( )
00077   : m_plotter ( 0 ), 
00078     m_x ( 0 ), 
00079     m_y ( 0 ),
00080     m_confid_count ( 0 )
00081 {
00082   m_deltaXSq.resize( 6 );
00083 
00084   m_deltaXSq[ 0 ] = 15.1;
00085   m_deltaXSq[ 1 ] = 18.4;
00086   m_deltaXSq[ 2 ] = 21.1;
00087   m_deltaXSq[ 3 ] = 23.5;
00088   m_deltaXSq[ 4 ] = 25.7;
00089   m_deltaXSq[ 5 ] = 27.8;
00090 }
00091 
00092 FunctionController::~FunctionController ( )
00093 {
00094 }
00095 
00096 FunctionController * FunctionController::instance ( )
00097 {
00098   if ( s_instance == 0 ) {
00099     s_instance = new FunctionController ( );
00100   }
00101   return s_instance;
00102 }
00103 
00104 const vector < string > &
00105 FunctionController::
00106 getFitterNames () const
00107 {
00108   FitterFactory * factory = FitterFactory::instance ();
00109 
00110   return factory -> names ();
00111 }
00112 
00113 const std::string &
00114 FunctionController::
00115 getDefaultFitter () const
00116 {
00117   FitterFactory * factory = FitterFactory::instance ();
00118 
00119   return factory -> getDefault ();
00120 }
00121 
00122 int FunctionController::
00123 getUniqueNonFunctionIndex ( const PlotterBase * plotter ) const
00124 {
00125   int index = -1;
00126   int count = 0;
00127 
00128   int number = plotter->getNumDataReps ();
00129   for ( int i = 0; i < number; i++ ) {
00130     DataRep * r = plotter->getDataRep ( i );
00131     FunctionRep * fr = dynamic_cast < FunctionRep * > ( r );
00132     if ( fr != 0 ) continue;
00133     index = i;
00134     count++;
00135   }
00136   if ( count != 1 ) index = -1;
00137 
00138   return index;
00139 }
00140 
00141 void
00142 FunctionController::
00143 findFunctions ( const PlotterBase * plotter ) const
00144 {
00145   m_func_reps.clear ();
00146   int number = plotter->getNumDataReps ();
00147 
00148   for ( int i = 0; i < number; i++ ) {
00149     DataRep * rep = plotter->getDataRep ( i );
00150     FunctionRep * frep = dynamic_cast < FunctionRep * > ( rep );
00151     if ( frep == 0 ) continue;
00152     m_func_reps.push_back ( frep );
00153   }
00154 }
00155 
00156 FunctionRep * 
00157 FunctionController::
00158 createFunctionRep ( const std::string & name,
00159                     DataRep * rep )
00160 {
00161   FunctionFactory * factory = FunctionFactory::instance ();
00162   FunctionBase * function = factory -> create ( name );
00163   if ( function == 0 ) {
00164     string what ( "FunctionController: Unable to create function of type `" );
00165     what += name;
00166     what += "'";
00167     throw std::runtime_error ( what );
00168   }
00169 
00170   return createFunctionRep ( function, rep );
00171 }
00172 
00173 FunctionRep * 
00174 FunctionController::
00175 createFunctionRep ( FunctionBase * function, 
00176                     DataRep * rep )
00177 {
00178   FunctionRep * frep = 0;
00179 
00180   if ( function -> isComposite () ) {
00181     frep = new CompositeFunctionRep ( function, rep );
00182   }
00183   else {
00184     unsigned int dims = function -> dimensions ();
00185     if ( dims == 2 ) {
00186       frep = new FunctionRep2 ( function, rep );
00187     }
00188     else {
00189       frep = new FunctionRep1 ( function, rep );
00190     }
00191   }
00192 
00193   if ( rep != 0 ) {// only makes sense to add a fitter if function has a target
00194     FitterFactory * factory = FitterFactory::instance ();
00195     const string & name = factory -> getDefault ();
00196     bool ok = setFitter ( frep, name );
00197     if ( ok == false ) {
00198       delete frep;
00199       frep = 0;
00200     }
00201   }
00202 
00203   return frep;
00204 }
00205 
00209 FunctionBase *
00210 FunctionController::
00211 addFunction ( PlotterBase * plotter, const std::string & name )
00212 {
00213   DataRep * rep = plotter -> getDataRep ( 0 );
00214   FunctionRep * frep = getFunctionRep ( plotter );// might be null
00215   addFunction ( plotter, name, frep, rep );
00216   
00217   frep = getFunctionRep ( plotter ); // now should be LinearSum
00218 
00219   return frep -> getFunction ();
00220 }
00221 
00225 FunctionRep *
00226 FunctionController::
00227 addFunction ( PlotterBase * plotter,
00228               const std::string & name,
00229               FunctionRep * frep,
00230               DataRep * rep )
00231 {
00232   FunctionRep * func_rep = createFunctionRep ( name, rep );
00233   func_rep -> initializeWith ( rep );
00234 
00235   return addFunctionRep ( plotter, rep, frep, func_rep );
00236 }
00237 
00238 void
00239 FunctionController::
00240 addFunction ( PlotterBase * plotter, 
00241               FunctionRep * func_rep )
00242 {
00243   int index = plotter->activePlotIndex ();
00244 
00245   if ( index < 0 ) {
00246     index = getUniqueNonFunctionIndex ( plotter );
00247   }
00248   if ( index >= 0 ) {
00249     plotter->setActivePlot ( index, false );
00250     DataRep * drep = plotter->getDataRep ( index );
00251     func_rep -> initializeWith ( drep );
00252     fillFunctionReps ( m_func_reps, plotter, drep );
00253     FunctionRep * frep = 0;
00254     if ( m_func_reps.empty () == false ) {
00255       frep = m_func_reps.front ();
00256     }
00257     addFunctionRep ( plotter, drep, frep, func_rep );
00258   }
00259 }
00260 
00261 
00262 void
00263 FunctionController::
00264 addDataRep ( PlotterBase * plotter, DataRep * rep )
00265 {
00266   FunctionRep * frep = dynamic_cast < FunctionRep * > ( rep );
00267   if ( frep == 0 ) { // not a function rep
00268     DisplayController * controller = DisplayController::instance ();
00269     controller -> addDataRep ( plotter, rep );
00270   }
00271   else { // a function rep
00272     if ( frep -> isComposite () == false ) {
00273       DataRep * target = frep -> getTarget ();
00274       addFunctionRep ( plotter, target, 0, frep );
00275     }
00276   }
00277 }
00278 
00279 void
00280 FunctionController::
00281 setTupleCut ( FunctionRep * rep )
00282 {
00283   rep -> setTupleCut ();
00284   rep -> setCutRange ( true );
00285 }
00286 
00287 void
00288 FunctionController::
00289 setTupleCut ( const PlotterBase * plotter, DataRep * data_rep )
00290 {
00291   data_rep -> addCut ();
00292   FunctionRep * func_rep = getFunctionRep ( plotter, data_rep );
00293   if ( func_rep != 0 ) {
00294     setTupleCut ( func_rep );
00295   }
00296 }
00297 
00298 void
00299 FunctionController::
00300 removeTupleCut ( const PlotterBase * plotter, DataRep * data_rep )
00301 {
00302   FunctionRep * func_rep = getFunctionRep ( plotter, data_rep );
00303   if ( func_rep != 0 ) {
00304     func_rep -> removeCut ();
00305   }
00306 }
00307 
00308 FunctionRep *
00309 FunctionController::
00310 addFunctionRep ( PlotterBase * plotter,
00311                  DataRep * drep,
00312                  FunctionRep * frep,
00313                  FunctionRep * func_rep )
00314 {
00315   FunctionRep * return_rep = frep;
00316 
00317   plotter->addDataRep ( func_rep );
00318 
00319   fillFunctionReps ( m_func_reps, plotter, drep );
00320 
00321   FitterFactory * factory = FitterFactory::instance ();
00322   const string & name = factory -> getDefault ();
00323   setFitter ( func_rep, name );
00324 
00325   if ( frep != 0 ) {
00326     if ( frep -> isComposite () ) {
00327       frep -> addToComposite ( func_rep );
00328     }
00329     else {
00330       FunctionRep * composite = createFunctionRep ( "Linear Sum", drep );
00331       composite -> initializeWith ( drep );
00332       plotter -> addDataRep ( composite );
00333       setFitter ( composite, name );
00334       setTupleCut ( func_rep );
00335       composite -> addToComposite ( frep );
00336       composite -> addToComposite ( func_rep );
00337       return_rep = composite;
00338     }
00339   }
00340   func_rep->notifyObservers ();
00341   setTupleCut ( func_rep );
00342   if ( return_rep == 0 ) return_rep = func_rep;
00343 
00344   return return_rep;
00345 }
00346 
00347 void
00348 FunctionController::
00349 removeFunction ( PlotterBase * plotter,
00350                  FunctionRep * frep )
00351 {
00352   if ( frep -> isComposite () ) {
00353     CompositeFunctionRep * composite
00354       = dynamic_cast < CompositeFunctionRep * > ( frep );
00355     const vector < FunctionRep * > & freps = composite -> getFunctionReps ();
00356     unsigned int size = freps.size();
00357 
00358 // remove from end so we don't destroy validty of reference vector
00359     for ( int i = size -1; i >= 0; i-- ) {
00360       FunctionRep * rep = freps[i];
00361       plotter -> removeDataRep ( rep );
00362       composite -> removeFromComposite ( rep );
00363       delete rep;
00364     }
00365     plotter -> removeDataRep ( frep );
00366     delete frep;
00367   }
00368   else {
00369     fillFunctionReps ( m_func_reps, plotter, 0 ); // get all
00370     vector < FunctionRep * > :: iterator it 
00371       = find ( m_func_reps.begin(), m_func_reps.end (),
00372                frep );
00373     assert ( it != m_func_reps.end () );
00374     m_func_reps.erase ( it );
00375  
00376     plotter -> removeDataRep ( frep );
00377     for ( unsigned int i = 0; i < m_func_reps.size (); i++ ) {
00378       FunctionRep * rep = m_func_reps[i];
00379       if ( rep -> isComposite () ) {
00380         rep -> removeFromComposite ( frep );
00381       }
00382     }
00383     // Remove all composites who now have one or none functions within them.
00384   
00385     vector < FunctionRep * >::iterator iter = m_func_reps.begin ();
00386   
00387     // This while loop is tricky because modifying the vector that we
00388     // are iterating over.   
00389     while ( iter != m_func_reps.end() ){
00390       FunctionRep * rep = *iter;
00391       if ( rep -> isComposite () ) {
00392         CompositeFunctionRep * composite
00393           = dynamic_cast < CompositeFunctionRep * > ( rep );
00394         if ( composite -> count() < 2 ) {
00395           plotter->removeDataRep ( composite );
00396           const CompositeFunctionRep::FunctionRepList_t & reps 
00397             = composite -> getFunctionReps ();
00398           FunctionRep * crep = reps.front ();
00399           crep -> setInComposite ( false );
00400           delete composite;
00401           iter = m_func_reps.erase ( iter ); // updates iter 
00402         }
00403         else {
00404           iter ++;
00405         }
00406       }
00407       else {
00408         iter ++;
00409       }
00410     }
00411     delete frep; // the functionrep
00412   }
00413 }
00414 
00415 bool
00416 FunctionController::
00417 hasFunction ( const PlotterBase * plotter, const DataRep * rep )
00418 {
00419   bool yes = false;
00420   assert ( plotter != 0 );
00421 
00422   fillFunctionReps ( m_func_reps, plotter, rep );
00423   yes = m_func_reps.empty () == false;
00424 
00425   return yes;
00426 }
00427 
00428 const vector< string > & 
00429 FunctionController::
00430 functionNames ( PlotterBase * plotter, DataRep * rep )
00431 {
00432   fillFunctionReps ( m_func_reps, plotter, rep );
00433 
00434   m_f_names.clear ();
00435   vector< FunctionRep * >:: const_iterator it = m_func_reps.begin ();
00436 
00437   for ( ; it != m_func_reps.end (); ++it ) {
00438     FunctionRep * fp = *it;
00439     if ( rep != 0 && fp->getTarget() != rep ) continue;
00440     FunctionBase * function = fp->getFunction ();
00441     const string & name = function->name ();
00442     m_f_names.push_back ( name );
00443   }
00444 
00445   return m_f_names;
00446 }
00447 
00448 void
00449 FunctionController::
00450 fillFunctionReps ( std::vector < FunctionRep * > & freps,
00451                    const PlotterBase * plotter,
00452                    const DataRep * drep ) const
00453 {
00454   freps.clear ();
00455   int number = plotter -> getNumDataReps ();
00456   for ( int i = 0; i < number; i++ ) {
00457     DataRep * rep = plotter -> getDataRep ( i );
00458     FunctionRep * frep = dynamic_cast < FunctionRep * > ( rep );
00459     if ( frep != 0 ) {
00460       if ( drep != 0 ) {
00461         if ( frep -> getTarget () == drep ) {
00462           freps.push_back ( frep );
00463         }
00464       }
00465       else { // any target
00466         freps.push_back ( frep );
00467       }
00468     }
00469   }
00470 }
00471 
00472 void
00473 FunctionController::
00474 fillTopLevelFunctionReps ( std::vector < FunctionRep * > & freps,
00475                            const PlotterBase * plotter,
00476                            const DataRep * drep ) const
00477 {
00478   freps.clear ();
00479   fillFunctionReps ( m_func_reps, plotter, drep );
00480   unsigned int size = m_func_reps.size ();
00481   for ( unsigned int i = 0; i < size; i++ ) {
00482     FunctionRep * rep = m_func_reps[i];
00483     if ( rep -> isInComposite () == false ) {
00484       freps.push_back ( rep );
00485     }
00486   }
00487 }
00488 
00489 void
00490 FunctionController::
00491 setErrorsFromComposite ( const PlotterBase * plotter,
00492                          const FunctionRep * composite )
00493 {
00494   const vector < double > & errors = composite -> principleErrors();
00495   vector < double >::const_iterator begin = errors.begin();
00496 
00497   DataRep * target = composite -> getTarget ();
00498   fillFunctionReps ( m_func_reps, plotter, target );
00499 
00500   vector < FunctionRep * >::iterator first = m_func_reps.begin();
00501 
00502   while ( first != m_func_reps.end() ) {
00503     FunctionRep * rep = *first++;
00504     if ( rep -> isComposite() ) continue;
00505     const vector < double > & t = rep -> principleErrors ();
00506     unsigned int number = t.size();
00507     rep->setPrincipleErrors ( begin, 
00508                               begin + number );
00509     begin += number;
00510   }
00511 
00512 }
00513 
00514 bool
00515 FunctionController::
00516 fitFunction ( PlotterBase * plotter, unsigned int )
00517 {
00518   FunctionRep * rep = getFunctionRep ( plotter );
00519 
00520   return fitFunction ( plotter, rep );
00521 }
00522 
00523 FunctionRep *
00524 FunctionController::
00525 getComposite ( const PlotterBase * plotter, FunctionRep * rep )
00526 {
00527   const DataRep * target = rep -> getTarget ();
00528   fillFunctionReps ( m_func_reps, plotter, target );
00529   for ( unsigned int i = 0; i < m_func_reps.size (); i++ ) {
00530     FunctionRep * fr = m_func_reps[i];
00531     if ( fr -> isComposite () ) {
00532       CompositeFunctionRep * composite 
00533         = dynamic_cast < CompositeFunctionRep * > ( fr );
00534       bool yes = composite -> isMember ( rep );
00535       if ( yes ) {
00536         rep = composite;
00537         break;
00538       }
00539     }
00540   }
00541 
00542   return rep;
00543 }
00544 
00545 bool
00546 FunctionController::
00547 fitFunction ( PlotterBase * plotter, FunctionRep * rep )
00548 {
00549   rep = getComposite ( plotter, rep );
00550   bool ok = rep -> fitFunction ();
00551   if ( rep -> isComposite () ) {
00552     setErrorsFromComposite ( plotter, rep );
00553   }
00554 
00555   return ok;
00556 }
00557 
00558 bool
00559 FunctionController::
00560 tryFitFunction ( PlotterBase * plotter, FunctionRep * func_rep )
00561 {
00562   saveParameters ( plotter );
00563   double chi_sq = func_rep -> objectiveValue ();
00564   bool ok = fitFunction ( plotter, func_rep );
00565   if ( ok ) {
00566     double new_chi_sq = func_rep -> objectiveValue ();
00567     if ( new_chi_sq > chi_sq ) {
00568       restoreParameters ( plotter );
00569     }
00570   } else {
00571     restoreParameters ( plotter );
00572   }
00573 
00574   return ok;
00575 }
00576 
00577 void
00578 FunctionController::
00579 setFitRange ( PlotterBase * plotter, const Range & range )
00580 {
00581   if ( hasFunction ( plotter, 0 ) ) {
00582     FunctionRep * frep = getFunctionRep ( plotter );
00583     frep -> setCutRange ( range );
00584     plotter -> update ();
00585   }
00586 }
00587 
00588 void
00589 FunctionController::
00590 setFitRange ( PlotterBase * plotter, double low, double high )
00591 {
00592   const Range range ( low, high );
00593   setFitRange ( plotter, range );
00594 }
00595 
00596 void FunctionController::saveParameters ( PlotterBase * plotter )
00597 {
00598   findFunctions ( plotter );
00599   vector< FunctionRep * >::iterator it = m_func_reps.begin ();
00600 
00601   for ( ; it != m_func_reps.end (); ++it ) {
00602     FunctionRep * frep = *it;
00603     FunctionBase * function = frep->getFunction ();
00604     if ( function->isComposite () ) continue;
00605     frep->saveParameters ();
00606   }
00607 }
00608 
00609 void FunctionController::restoreParameters ( PlotterBase * plotter )
00610 {
00611   findFunctions ( plotter );
00612   vector< FunctionRep * >::iterator it = m_func_reps.begin ();
00613 
00614   for ( ; it != m_func_reps.end (); ++it ) {
00615     FunctionRep * frep = *it;
00616     FunctionBase * function = frep->getFunction ();
00617     if ( function->isComposite () ) continue;
00618     frep->restoreParameters ();
00619   }
00620 
00621 }
00622 
00623 FunctionRep *
00624 FunctionController::
00625 getFunctionRep ( const PlotterBase * plotter ) const
00626 {
00627   FunctionRep * frep = 0;
00628   DataRep * rep = 0;
00629   int index = plotter->activePlotIndex ();
00630 
00631   if ( index >= 0 ) {
00632     rep = plotter -> getDataRep ( index );
00633   }
00634   else {
00635     rep = plotter -> getTarget ();
00636   }
00637 
00638   FunctionRep * test = dynamic_cast < FunctionRep * > ( rep );
00639   if ( test != 0 ) { // if ctrl-clicked, could be the function
00640     frep = test;
00641   }
00642   else {
00643     frep = getFunctionRep ( plotter, rep );
00644   }
00645 
00646   return frep;
00647 }
00648 
00652 FunctionRep *
00653 FunctionController::
00654 getFunctionRep ( const PlotterBase * plotter, const DataRep * datarep ) const
00655 {
00656   FunctionRep * frep = 0;
00657 
00658   fillFunctionReps ( m_func_reps, plotter, datarep );
00659 
00660   for ( unsigned int i = 0; i < m_func_reps.size(); i++ ) {
00661     FunctionRep * rep = m_func_reps[i];
00662     const DataRep * drep = rep -> getTarget ();
00663     if ( drep != datarep ) continue;
00664     frep = rep;
00665     if ( frep -> isComposite () ) break; // use composite if found
00666   }
00667 
00668   return frep;
00669 }
00670 
00671 ViewBase * FunctionController::
00672 createFuncView ( const ViewFactory * factory,
00673                  PlotterBase * plotter,
00674                  const std::string & name )
00675 {
00676   FunctionRep * frep = getFunctionRep ( plotter );
00677   assert ( frep != 0 );
00678 
00679   DisplayController * controller = DisplayController::instance ();
00680   string nullstring ("");
00681 
00682   return controller->createTextView ( factory, frep, name, nullstring );
00683 }
00684 
00685 bool
00686 FunctionController::
00687 setFitter ( const PlotterBase * plotter, const std::string & name )
00688 {
00689   FunctionRep * frep = getFunctionRep ( plotter );
00690   bool ok = false;
00691   if ( frep != 0 ) {
00692     ok = setFitter ( frep, name );
00693   }
00694 
00695   return ok;
00696 }
00697 
00698 const string &
00699 FunctionController::
00700 getFitterName ( const PlotterBase * plotter )
00701 {
00702   FunctionRep * frep = getFunctionRep ( plotter );
00703   Fitter * fitter = frep -> getFitter ();
00704 
00705   return fitter -> name ();
00706 }
00707 
00708 Fitter *
00709 FunctionController::
00710 getFitter ( const PlotterBase * plotter )
00711 {
00712   FunctionRep * frep = getFunctionRep ( plotter );
00713   Fitter * fitter = frep -> getFitter ();
00714 
00715   return fitter;
00716 }
00717 
00718 
00719 bool
00720 FunctionController::
00721 setFitter ( FunctionRep * frep, const std::string & name )
00722 {
00723   FitterFactory * factory = FitterFactory::instance ();
00724   Fitter * fitter = factory -> create ( name );
00725 
00726   return frep -> setFitter ( fitter );
00727 }
00728 
00729 void
00730 FunctionController::
00731 setDefaultFitter ( const std::string & name )
00732 {
00733   FitterFactory * factory = FitterFactory::instance ();
00734   factory -> setDefault ( name );
00735 }
00736 
00737 bool
00738 FunctionController::
00739 changeFitter ( const PlotterBase * plotter,
00740                const DataRep * drep,
00741                const std::string & name )
00742 {
00743   FitterFactory * factory = FitterFactory::instance ();
00744   Fitter * proto = factory -> prototype ( name );
00745 
00746   FunctionRep * frep = getFunctionRep ( plotter, drep );
00747   const FunctionBase * func = frep -> getFunction ();
00748   bool yes = proto -> isCompatible ( func );
00749   if ( yes ) {
00750     const Fitter * old_fitter = frep -> getFitter ();
00751     Fitter * new_fitter = factory -> create ( name );
00752     new_fitter -> copyFrom ( old_fitter );
00753     frep -> setFitter ( new_fitter );
00754   }
00755   return yes;
00756 }
00757 
00758 double
00759 FunctionController::
00760 getObjectiveValue ( const PlotterBase * plotter, const DataRep * datarep )
00761 {
00762   FunctionRep * rep = getFunctionRep ( plotter, datarep );
00763 
00764   return rep -> objectiveValue ();
00765 }
00766 
00767 const vector < vector < double > > &
00768 FunctionController::
00769 getCovarianceMatrix ( const PlotterBase * plotter )
00770 {
00771   FunctionRep * rep = getFunctionRep ( plotter );
00772 
00773   return rep -> covarianceMatrix ();
00774 }
00775 
00780 double
00781 FunctionController::
00782 getChiSquared ( const PlotterBase * plotter )
00783 {
00784   const DataRep * rep = plotter -> getTarget ();
00785 
00786   return getObjectiveValue ( plotter, rep );
00787 }
00788 
00789 int
00790 FunctionController::
00791 getDegreesOfFreedom ( const PlotterBase * plotter )
00792 {
00793   FunctionRep * rep = getFunctionRep ( plotter );
00794 
00795   return rep -> degreesOfFreedom ();
00796 }
00797 
00798 NTuple *
00799 FunctionController::
00800 createNTuple ( const PlotterBase * plotter, const FunctionRep * rep )
00801 {
00802   NTuple * ntuple = 0;
00803   if ( rep == 0 ) {
00804     int old_index = plotter -> activePlotIndex ();
00805     int index = getUniqueNonFunctionIndex ( plotter );
00806     if ( index < 0 ) return 0;
00807 
00808     int saved_index = plotter -> activePlotIndex ();
00809     PlotterBase * pb = const_cast < PlotterBase * > ( plotter );
00810     pb -> setActivePlot ( index, false ); // no redraw
00811     ntuple = plotter -> createNTuple ();
00812     pb -> setActivePlot ( saved_index, false );
00813     rep = getFunctionRep ( plotter );
00814     pb -> setActivePlot ( old_index, true );
00815   }
00816   else {
00817     const DataRep * target = rep -> getTarget ();
00818     ntuple = target -> createNTuple ();
00819   }
00820   if ( rep != 0 ) {
00821     const FunctionBase * function = rep -> getFunction ();
00822     vector < double > & x = ntuple -> getColumn ( 0 ); // X coordinate
00823     vector < double > & y = ntuple -> getColumn ( 1 ); // Y coordinate
00824     unsigned int size = ntuple -> rows ();
00825     vector < double > values ( size );
00826     vector < double > residuals ( size );
00827   
00828     for ( unsigned int i = 0; i < x.size(); i++ ) {
00829       values [i] = function -> operator () ( x[i] );
00830       residuals [i] =  y[i] - values [i];
00831     }
00832 
00833     ntuple -> addColumn ( function->name (), values );
00834     ntuple -> addColumn ( "Residuals", residuals );
00835   }
00836 
00837   return ntuple;
00838 }
00839 
00840 PlotterBase *
00841 FunctionController::
00842 createResidualsDisplay ( PlotterBase * plotter, const FunctionRep * frep )
00843 {
00844   NTuple * ntuple = createNTuple ( plotter, frep );
00845   ntuple -> setTitle ( plotter -> getTitle () );
00846   const string old_label ( "Residuals" );
00847   DataSourceController::instance () -> registerNTuple ( ntuple );
00848   vector < string > bindings ( 4 );
00849   bindings[0] = ntuple -> getLabelAt ( 0 );
00850   bindings[1] = old_label;
00851   bindings[2] = "nil";
00852   bindings[3] = ntuple -> getLabelAt ( 3 );
00853   bool axis_scaled = plotter -> isAxisScaled ( Axes::Y );
00854   if ( axis_scaled ) {
00855     int i = ntuple -> indexOf ( old_label );
00856     const string new_label ("Residuals (scaled)"); 
00857     ntuple -> setLabelAt ( new_label , i );
00858     bindings[1] = new_label;
00859   }
00860   DisplayController * controller = DisplayController::instance ();
00861 
00862 // Make scaling of x-axis match that of the original plot.
00863   PlotterBase * new_plotter = controller -> createDisplay ( "XY Plot",
00864                                                             *ntuple,
00865                                                             bindings );
00866   if ( axis_scaled ) {
00867     double factor = plotter -> getScaleFactor ( Axes::Y );
00868     new_plotter -> setScaleFactor ( Axes::Y, factor );
00869   }
00870   controller->setLog( new_plotter, "x", 
00871                       controller -> getLog( plotter, "x" ) );
00872   
00873   return new_plotter;
00874 }
00875 
00876 PlotterBase *
00877 FunctionController::
00878 createNewEllipsoidDisplay ( PlotterBase * masterPlot, FunctionRep * rep )
00879 {
00880   int n = 100;
00881   double xmin, xmax, ymin, ymax;
00882   
00883   // Set up the NTuple
00884   NTuple* ntuple = new( NTuple );
00885   ellipsoidNTuple ( masterPlot, rep, ntuple, n, xmin, xmax, ymin, ymax );
00886   ntuple -> setTitle ( masterPlot -> getTitle () );
00887   
00888   vector < string > bindings ( 1 );
00889   bindings[0] = ntuple ->  getLabelAt ( 0 );
00890   
00891   DisplayController * dcontroller = DisplayController::instance ();
00892 
00893   // Z plot. X and Y coordinates are merely indices
00894   PlotterBase * ellipsePlot = dcontroller ->
00895     createDisplay ( "Image", *ntuple, bindings );
00896   
00897   // Rescale and translate the Z plot so that  the plot
00898   // becomes contained in the rectangle bound by xmin, xmax
00899   // ymin and ymax.
00900   ellipsePlot -> setBinWidth( "x", ( xmax - xmin ) / ( n - 1.0 ) );
00901   ellipsePlot -> setBinWidth( "y", ( ymax - ymin ) / ( n - 1.0 ) );
00902   
00903   ellipsePlot -> setOffsets ( xmin, ymin );
00904   
00905   ellipsePlot -> setRange ( string( "X" ), xmin, xmax );
00906   ellipsePlot -> setRange ( string( "Y" ), ymin, ymax );
00907   
00908   // Establish the relation between source masterPlot and
00909   // this new plot. This relationship shall be exploited
00910   // when we refresh the error plots.
00911   int index = dcontroller -> activeDataRepIndex( masterPlot );
00912   ellipsePlot -> setParentDataRepIndex( index );
00913   ellipsePlot -> setParentPlotter( masterPlot );
00914   
00915   return ellipsePlot;
00916 }
00917 
00918 void
00919 FunctionController::
00920 ellipsoidNTuple ( PlotterBase* masterPlot, FunctionRep * frep,
00921                   NTuple* nt, int n,
00922                   double & xmin, double & xmax,
00923                   double & ymin, double & ymax )
00924 {
00925   string text ( "Confidence ellipse " );
00926   text += String::convert ( ++m_confid_count );
00927   nt -> setName ( text );
00928 
00929   DataSourceController * controller = DataSourceController::instance ();
00930   controller -> registerNTuple ( nt );
00931 
00932   assert( frep );
00933   
00934   // Get projected covariance i.e. take the sub-matrix
00935   // out of the covariance matrix which corresponds
00936   // to the 2 parameters of interest whose correlation
00937   // we would like to see.
00938   vector< vector< double > > Sigma( 2 );
00939   Sigma[0].resize( 2, 0.0 );
00940   Sigma[1].resize( 2, 0.0 );
00941 
00942   const vector < vector < double > > & covariance 
00943     = frep -> covarianceMatrix ();  
00944 
00945   Sigma[0][0] = covariance [m_x][m_x];
00946   Sigma[0][1] = covariance [m_x][m_y];
00947   Sigma[1][0] = covariance [m_y][m_x];
00948   Sigma[1][1] = covariance [m_y][m_y];
00949 
00950   // Invert the projected covariance 
00951   vector< vector< double > > SigmaInv;
00952   invertMatrix( Sigma, SigmaInv );
00953   
00954   // Decide the center of the ellipse to be parameter
00955   // value of the 2 parameters of interest whose correlation
00956   // we would like to see.
00957   vector< double > xbar( 2 );
00958   
00959   const Fitter * fitter = getFitter ( masterPlot );
00960   vector < double > free_parms;
00961   fitter -> fillFreeParameters ( free_parms );
00962   xbar[ 0 ] = free_parms [ m_x ];
00963   xbar[ 1 ] = free_parms [ m_y ];
00964   
00965   // Get the bounding ellipse and its bounds. Bounding ellipse is the
00966   // 99.99% confidence ellipsoid. For mu = 2 the delta chi-square turns from 
00967   // Numerical Recipes in C as 18.4. etc etc.
00968   unsigned int mu = free_parms.size();
00969   NTuple * boundingTuple = ellipse( xbar, SigmaInv, m_deltaXSq[ mu - 1 ] ); 
00970   
00971   xmin = boundingTuple -> minElement ( 0 );
00972   xmax = boundingTuple -> maxElement ( 0 );
00973   ymin = boundingTuple -> minElement ( 1 );
00974   ymax = boundingTuple -> maxElement ( 1 );
00975   
00976   delete boundingTuple; // Served its purpose
00977   
00978   // Create the NTuple for the contour plot
00979   // This NTuple is to be build column-wise
00980   // keeping in view the way Z-plot works.
00981   vector< double > p( n * n ), a( 2 );
00982   
00983   double dx    = ( xmax - xmin ) / ( n - 1.0 );
00984   double dy    = ( ymax - ymin ) / ( n - 1.0 );
00985   double delta = 0.0;
00986   
00987   for( int i = 0; i < n; i++ )
00988     for( int j = 0; j < n; j++ )
00989       {
00990         a[ 0 ] = xmin + i * dx - xbar[ 0 ] ;
00991         a[ 1 ] = ymin + j * dy - xbar[ 1 ];
00992         delta  = quadraticProduct( SigmaInv, a );
00993         p[ n * i + j ] = 1 - gammq( mu / 2.0, delta / 2.0 ); 
00994       }
00995 
00996   nt -> clear();
00997   nt -> addColumn( "Percent", 100 * p );
00998   
00999   return;
01000 }
01001 
01002 PlotterBase *
01003 FunctionController::
01004 refreshEllipsoidDisplay ( PlotterBase * plot, FunctionRep * frep )
01005 {
01006   int n = 100;
01007   double xmin, xmax, ymin, ymax;
01008   
01009   PlotterBase* masterPlot = plot->getParentPlotter();
01010 
01011   // Set up the NTuple
01012   NTuple* ntuple = new( NTuple );
01013   ellipsoidNTuple ( masterPlot, frep, ntuple, n, xmin, xmax, ymin, ymax );
01014   ntuple -> setTitle ( masterPlot -> getTitle () );
01015   
01016   // First get the selected DataRep from the plotter. The DataRep 
01017   // will have a projector that is a NTupleProjector. It doesn't
01018   // know, but we do, so we downcast and set the NTuple.
01019   DataRep * drep = plot -> selectedDataRep();
01020   NTupleProjector * ntProjector =
01021     dynamic_cast < NTupleProjector * > ( drep -> getProjector() );
01022   ntProjector -> setNTuple ( ntuple );
01023   
01024   NTuple * nt = const_cast < NTuple * > ( ntuple );
01025   nt  -> addObserver ( ntProjector );
01026     
01027   // Rescale and translate the Z plot so that  the plot
01028   // becomes contained in the rectangle bound by xmin, xmax
01029   // ymin and ymax.
01030   plot -> setBinWidth( "x", ( xmax - xmin ) / ( n - 1.0 ) );
01031   plot -> setBinWidth( "y", ( ymax - ymin ) / ( n - 1.0 ) );
01032   
01033   plot -> setOffsets ( xmin, ymin );
01034   
01035   plot -> setRange ( string( "X" ), xmin, xmax );
01036   plot -> setRange ( string( "Y" ), ymin, ymax );
01037     
01038   return plot;
01039 }
01040 
01041 NTuple * 
01042 FunctionController::
01043 ellipse ( const std::vector< double > & xbar,
01044           std::vector< std::vector < double > > & SigmaInv,
01045           double Csquare )
01046 {
01047   // First argument should be a 2 D vector, the center of the ellipse //
01048   assert( xbar.size() == 2 );
01049 
01050   // Second argument should be a 2 x 2 SPD matrix i.e. two rows two cols //
01051   assert( SigmaInv.size() == 2 );
01052   assert( SigmaInv[0].size() == 2 || SigmaInv[1].size() == 2 );
01053   
01054   // Second argument should be a 2 x 2 -- Symmetric --  PD matrix
01055   // Under numerical round-offs it might not be true so we take mean of the
01056   // entries.
01057   double temp = ( SigmaInv[0][1] + SigmaInv[1][0] ) / 2;
01058   SigmaInv[0][1] = temp;
01059   SigmaInv[1][0] = temp;
01060   
01061   // Third argument should be a positive scalar "c^2", Square of "Radius" // 
01062   assert( Csquare > DBL_EPSILON );
01063   
01064   // The eigenvalues of the SigmaInv matrix
01065   double b = -( SigmaInv[0][0] + SigmaInv[1][1] );
01066   double c = SigmaInv[0][0] * SigmaInv[1][1] - SigmaInv[0][1] * SigmaInv[1][0];
01067 
01068   double lambda1 = ( -b + sqrt( b * b - 4 * c ) ) / 2;
01069   double lambda2 = ( -b - sqrt( b * b - 4 * c ) ) / 2;
01070     
01071   // Determining the angles the ellipse axis make with X-axis
01072   double alpha1 =  atan( (SigmaInv[0][0] - lambda1) / SigmaInv[0][1] );
01073   double a      =  sqrt( Csquare / lambda1 );
01074   b             =  sqrt( Csquare / lambda2 );
01075 
01076   // Creating a rotation matrix
01077   double Rot00 =  cos( alpha1 );
01078   double Rot11 =  cos( alpha1 );
01079   double Rot01 = -sin( alpha1 );
01080   double Rot10 =  sin( alpha1 );
01081 
01082   // Generating N-points on the ellipse
01083   int N = 100;
01084   vector< double > x, y;
01085   
01086   x.resize( N, 0.0 );
01087   y.resize( N, 0.0 );
01088   
01089   for( int i = 0; i < N; i++)
01090     {
01091       // Unrotated untranslated point
01092       double theta = ( 2 * M_PI * i ) / ( (double) (N-1) );
01093       double x0 = a * cos( theta );
01094       double x1 = b * sin( theta );
01095       
01096       //x = Rot * x; i.e. force in a rotation
01097       double xrot0 = Rot00 * x0 + Rot01 * x1;
01098       double xrot1 = Rot10 * x0 + Rot11 * x1;
01099       
01100       //x = x + xbar; i.e. force in a translation
01101       x[i] = xrot0 + xbar[0];
01102       y[i] = xrot1 + xbar[1];
01103     }
01104   
01105   // Create an NTuple out of the above vector
01106   NTuple* ntuple = new( NTuple );
01107   ntuple -> addColumn( "X", x );
01108   ntuple -> addColumn( "Y", y );
01109   
01110   return ntuple;
01111 }
01112 
01113 int FunctionController::setEllpsoidParamIndex( hippodraw::Axes::Type axes,
01114                                                int index )
01115 {
01116   assert( axes == Axes::X || axes == Axes::Y );
01117 
01118   if( axes == Axes::X )
01119     m_x = index;
01120   
01121   if( axes == Axes::Y )
01122     m_y = index;
01123       
01124   return EXIT_SUCCESS;
01125 }
01126 
01127 bool
01128 FunctionController::
01129 functionExists ( const std::string & name )
01130 {
01131   FunctionFactory * factory = FunctionFactory::instance ();
01132 
01133   return factory -> exists ( name );
01134 }
01135 
01136 bool
01137 FunctionController::
01138 isCompatible ( const std::string & function,
01139                const std::string & fitter )
01140 {
01141   bool yes = true;
01142 
01143   FunctionFactory * fun_fac = FunctionFactory::instance ();
01144   FunctionBase * fun_proto = fun_fac -> prototype ( function );
01145 
01146   FitterFactory * fit_fac = FitterFactory::instance ();
01147   Fitter * fit_proto = fit_fac -> prototype ( fitter );
01148 
01149   yes = fit_proto -> isCompatible ( fun_proto );
01150 
01151   return yes;
01152 }
01153 
01154 const vector < string > &
01155 FunctionController::
01156 getFunctionNames () const
01157 {
01158   return FunctionFactory::instance () -> names ();
01159 }

Generated for HippoDraw Class Library by doxygen