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

Generated for HippoDraw Class Library by doxygen