FunctionController.cxx
Go to the documentation of this file.
1 
12 #ifdef _MSC_VER
13 #include "msdevstudio/MSconfig.h"
14 #endif
15 
16 #include "FunctionController.h"
17 
18 #include "Gammaq.h"
19 #include "DisplayController.h"
20 
23 #include "datareps/FunctionRep1.h"
24 #include "datareps/FunctionRep2.h"
25 
27 #include "datasrcs/NTuple.h"
28 #include "datasrcs/TupleCut.h"
29 
30 #include "functions/FunctionBase.h"
32 
33 #include "minimizers/Fitter.h"
35 #include "minimizers/NumLinAlg.h"
36 
37 #include "pattern/string_convert.h"
38 #include "plotters/PlotterBase.h"
42 
43 #include <algorithm>
44 #include <functional>
45 #include <iterator>
46 #include <stdexcept>
47 
48 #include <cmath>
49 #include <cassert>
50 
51 #ifdef ITERATOR_MEMBER_DEFECT
52 using namespace std;
53 #else
54 using std::distance;
55 using std::find;
56 using std::find_if;
57 using std::list;
58 using std::mem_fun;
59 using std::string;
60 using std::vector;
61 using std::sqrt;
62 using std::cos;
63 using std::sin;
64 using std::atan;
65 using std::min;
66 using std::max;
67 using std::abs;
68 #endif
69 
70 using namespace hippodraw;
71 
72 using namespace Numeric;
73 
74 FunctionController * FunctionController::s_instance = 0;
75 
76 FunctionController::FunctionController ( )
77  : m_plotter ( 0 ),
78  m_x ( 0 ),
79  m_y ( 0 ),
80  m_confid_count ( 0 )
81 {
82  m_deltaXSq.resize( 6 );
83 
84  m_deltaXSq[ 0 ] = 15.1;
85  m_deltaXSq[ 1 ] = 18.4;
86  m_deltaXSq[ 2 ] = 21.1;
87  m_deltaXSq[ 3 ] = 23.5;
88  m_deltaXSq[ 4 ] = 25.7;
89  m_deltaXSq[ 5 ] = 27.8;
90 }
91 
93 {
94 }
95 
97 {
98  if ( s_instance == 0 ) {
100  }
101  return s_instance;
102 }
103 
104 const vector < string > &
107 {
108  FitterFactory * factory = FitterFactory::instance ();
109 
110  return factory -> names ();
111 }
112 
113 const std::string &
116 {
117  FitterFactory * factory = FitterFactory::instance ();
118 
119  return factory -> getDefault ();
120 }
121 
123 getUniqueNonFunctionIndex ( const PlotterBase * plotter ) const
124 {
125  int index = -1;
126  int count = 0;
127 
128  int number = plotter->getNumDataReps ();
129  for ( int i = 0; i < number; i++ ) {
130  DataRep * r = plotter->getDataRep ( i );
131  FunctionRep * fr = dynamic_cast < FunctionRep * > ( r );
132  if ( fr != 0 ) continue;
133  index = i;
134  count++;
135  }
136  if ( count != 1 ) index = -1;
137 
138  return index;
139 }
140 
141 void
143 findFunctions ( const PlotterBase * plotter ) const
144 {
145  m_func_reps.clear ();
146  int number = plotter->getNumDataReps ();
147 
148  for ( int i = 0; i < number; i++ ) {
149  DataRep * rep = plotter->getDataRep ( i );
150  FunctionRep * frep = dynamic_cast < FunctionRep * > ( rep );
151  if ( frep == 0 ) continue;
152  m_func_reps.push_back ( frep );
153  }
154 }
155 
156 FunctionRep *
158 createFunctionRep ( const std::string & name,
159  DataRep * rep )
160 {
162  FunctionBase * function = factory -> create ( name );
163  if ( function == 0 ) {
164  string what ( "FunctionController: Unable to create function of type `" );
165  what += name;
166  what += "'";
167  throw std::runtime_error ( what );
168  }
169 
170  return createFunctionRep ( function, rep );
171 }
172 
173 FunctionRep *
176  DataRep * rep )
177 {
178  FunctionRep * frep = 0;
179 
180  if ( function -> isComposite () ) {
181  frep = new CompositeFunctionRep ( function, rep );
182  }
183  else {
184  unsigned int dims = function -> dimensions ();
185  if ( dims == 2 ) {
186  frep = new FunctionRep2 ( function, rep );
187  }
188  else {
189  frep = new FunctionRep1 ( function, rep );
190  }
191  }
192 
193  if ( rep != 0 ) {// only makes sense to add a fitter if function has a target
194  FitterFactory * factory = FitterFactory::instance ();
195  const string & name = factory -> getDefault ();
196  bool ok = setFitter ( frep, name );
197  if ( ok == false ) {
198  delete frep;
199  frep = 0;
200  }
201  }
202 
203  return frep;
204 }
205 
209 FunctionBase *
211 addFunction ( PlotterBase * plotter, const std::string & name )
212 {
213  DataRep * rep = plotter -> getDataRep ( 0 );
214  FunctionRep * frep = getFunctionRep ( plotter );// might be null
215  addFunction ( plotter, name, frep, rep );
216 
217  frep = getFunctionRep ( plotter ); // now should be LinearSum
218 
219  return frep -> getFunction ();
220 }
221 
225 FunctionRep *
228  const std::string & name,
229  FunctionRep * frep,
230  DataRep * rep )
231 {
232  FunctionRep * func_rep = createFunctionRep ( name, rep );
233  func_rep -> initializeWith ( rep );
234 
235  return addFunctionRep ( plotter, rep, frep, func_rep );
236 }
237 
238 void
241  FunctionRep * func_rep )
242 {
243  int index = plotter->activePlotIndex ();
244 
245  if ( index < 0 ) {
246  index = getUniqueNonFunctionIndex ( plotter );
247  }
248  if ( index >= 0 ) {
249  plotter->setActivePlot ( index, false );
250  DataRep * drep = plotter->getDataRep ( index );
251  func_rep -> initializeWith ( drep );
252  fillFunctionReps ( m_func_reps, plotter, drep );
253  FunctionRep * frep = 0;
254  if ( m_func_reps.empty () == false ) {
255  frep = m_func_reps.front ();
256  }
257  addFunctionRep ( plotter, drep, frep, func_rep );
258  }
259 }
260 
261 
262 void
264 addDataRep ( PlotterBase * plotter, DataRep * rep )
265 {
266  FunctionRep * frep = dynamic_cast < FunctionRep * > ( rep );
267  if ( frep == 0 ) { // not a function rep
269  controller -> addDataRep ( plotter, rep );
270  }
271  else { // a function rep
272  if ( frep -> isComposite () == false ) {
273  DataRep * target = frep -> getTarget ();
274  addFunctionRep ( plotter, target, 0, frep );
275  }
276  }
277 }
278 
279 void
282 {
283  rep -> setTupleCut ();
284  rep -> setCutRange ( true );
285 }
286 
287 void
289 setTupleCut ( const PlotterBase * plotter, DataRep * data_rep )
290 {
291  data_rep -> addCut ();
292  FunctionRep * func_rep = getFunctionRep ( plotter, data_rep );
293  if ( func_rep != 0 ) {
294  setTupleCut ( func_rep );
295  }
296 }
297 
298 void
300 removeTupleCut ( const PlotterBase * plotter, DataRep * data_rep )
301 {
302  FunctionRep * func_rep = getFunctionRep ( plotter, data_rep );
303  if ( func_rep != 0 ) {
304  func_rep -> removeCut ();
305  }
306 }
307 
308 FunctionRep *
311  DataRep * drep,
312  FunctionRep * frep,
313  FunctionRep * func_rep )
314 {
315  FunctionRep * return_rep = frep;
316 
317  plotter->addDataRep ( func_rep );
318 
319  fillFunctionReps ( m_func_reps, plotter, drep );
320 
321  FitterFactory * factory = FitterFactory::instance ();
322  const string & name = factory -> getDefault ();
323  setFitter ( func_rep, name );
324 
325  if ( frep != 0 ) {
326  if ( frep -> isComposite () ) {
327  frep -> addToComposite ( func_rep );
328  }
329  else {
330  FunctionRep * composite = createFunctionRep ( "Linear Sum", drep );
331  composite -> initializeWith ( drep );
332  plotter -> addDataRep ( composite );
333  setFitter ( composite, name );
334  setTupleCut ( func_rep );
335  composite -> addToComposite ( frep );
336  composite -> addToComposite ( func_rep );
337  return_rep = composite;
338  }
339  }
340  func_rep->notifyObservers ();
341  setTupleCut ( func_rep );
342  if ( return_rep == 0 ) return_rep = func_rep;
343 
344  return return_rep;
345 }
346 
347 void
350  FunctionRep * frep )
351 {
352  if ( frep -> isComposite () ) {
353  CompositeFunctionRep * composite
354  = dynamic_cast < CompositeFunctionRep * > ( frep );
355  const vector < FunctionRep * > & freps = composite -> getFunctionReps ();
356  unsigned int size = freps.size();
357 
358 // remove from end so we don't destroy validty of reference vector
359  for ( int i = size -1; i >= 0; i-- ) {
360  FunctionRep * rep = freps[i];
361  plotter -> removeDataRep ( rep );
362  composite -> removeFromComposite ( rep );
363  delete rep;
364  }
365  plotter -> removeDataRep ( frep );
366  delete frep;
367  }
368  else {
369  fillFunctionReps ( m_func_reps, plotter, 0 ); // get all
370  vector < FunctionRep * > :: iterator it
371  = find ( m_func_reps.begin(), m_func_reps.end (),
372  frep );
373  assert ( it != m_func_reps.end () );
374  m_func_reps.erase ( it );
375 
376  plotter -> removeDataRep ( frep );
377  for ( unsigned int i = 0; i < m_func_reps.size (); i++ ) {
378  FunctionRep * rep = m_func_reps[i];
379  if ( rep -> isComposite () ) {
380  rep -> removeFromComposite ( frep );
381  }
382  }
383  // Remove all composites who now have one or none functions within them.
384 
385  vector < FunctionRep * >::iterator iter = m_func_reps.begin ();
386 
387  // This while loop is tricky because modifying the vector that we
388  // are iterating over.
389  while ( iter != m_func_reps.end() ){
390  FunctionRep * rep = *iter;
391  if ( rep -> isComposite () ) {
392  CompositeFunctionRep * composite
393  = dynamic_cast < CompositeFunctionRep * > ( rep );
394  if ( composite -> count() < 2 ) {
395  plotter->removeDataRep ( composite );
397  = composite -> getFunctionReps ();
398  FunctionRep * crep = reps.front ();
399  crep -> setInComposite ( false );
400  delete composite;
401  iter = m_func_reps.erase ( iter ); // updates iter
402  }
403  else {
404  iter ++;
405  }
406  }
407  else {
408  iter ++;
409  }
410  }
411  delete frep; // the functionrep
412  }
413 }
414 
415 bool
417 hasFunction ( const PlotterBase * plotter, const DataRep * rep )
418 {
419  bool yes = false;
420  assert ( plotter != 0 );
421 
422  fillFunctionReps ( m_func_reps, plotter, rep );
423  yes = m_func_reps.empty () == false;
424 
425  return yes;
426 }
427 
428 const vector< string > &
430 functionNames ( PlotterBase * plotter, DataRep * rep )
431 {
432  fillFunctionReps ( m_func_reps, plotter, rep );
433 
434  m_f_names.clear ();
435  vector< FunctionRep * >:: const_iterator it = m_func_reps.begin ();
436 
437  for ( ; it != m_func_reps.end (); ++it ) {
438  FunctionRep * fp = *it;
439  if ( rep != 0 && fp->getTarget() != rep ) continue;
440  FunctionBase * function = fp->getFunction ();
441  const string & name = function->name ();
442  m_f_names.push_back ( name );
443  }
444 
445  return m_f_names;
446 }
447 
448 void
450 fillFunctionReps ( std::vector < FunctionRep * > & freps,
451  const PlotterBase * plotter,
452  const DataRep * drep ) const
453 {
454  freps.clear ();
455  int number = plotter -> getNumDataReps ();
456  for ( int i = 0; i < number; i++ ) {
457  DataRep * rep = plotter -> getDataRep ( i );
458  FunctionRep * frep = dynamic_cast < FunctionRep * > ( rep );
459  if ( frep != 0 ) {
460  if ( drep != 0 ) {
461  if ( frep -> getTarget () == drep ) {
462  freps.push_back ( frep );
463  }
464  }
465  else { // any target
466  freps.push_back ( frep );
467  }
468  }
469  }
470 }
471 
472 void
474 fillTopLevelFunctionReps ( std::vector < FunctionRep * > & freps,
475  const PlotterBase * plotter,
476  const DataRep * drep ) const
477 {
478  freps.clear ();
479  fillFunctionReps ( m_func_reps, plotter, drep );
480  unsigned int size = m_func_reps.size ();
481  for ( unsigned int i = 0; i < size; i++ ) {
482  FunctionRep * rep = m_func_reps[i];
483  if ( rep -> isInComposite () == false ) {
484  freps.push_back ( rep );
485  }
486  }
487 }
488 
489 void
492  const FunctionRep * composite )
493 {
494  const vector < double > & errors = composite -> principleErrors();
495  vector < double >::const_iterator begin = errors.begin();
496 
497  DataRep * target = composite -> getTarget ();
498  fillFunctionReps ( m_func_reps, plotter, target );
499 
500  vector < FunctionRep * >::iterator first = m_func_reps.begin();
501 
502  while ( first != m_func_reps.end() ) {
503  FunctionRep * rep = *first++;
504  if ( rep -> isComposite() ) continue;
505  const vector < double > & t = rep -> principleErrors ();
506  unsigned int number = t.size();
507  rep->setPrincipleErrors ( begin,
508  begin + number );
509  begin += number;
510  }
511 
512 }
513 
514 bool
516 fitFunction ( PlotterBase * plotter, unsigned int )
517 {
518  FunctionRep * rep = getFunctionRep ( plotter );
519 
520  return fitFunction ( plotter, rep );
521 }
522 
523 FunctionRep *
525 getComposite ( const PlotterBase * plotter, FunctionRep * rep )
526 {
527  const DataRep * target = rep -> getTarget ();
528  fillFunctionReps ( m_func_reps, plotter, target );
529  for ( unsigned int i = 0; i < m_func_reps.size (); i++ ) {
530  FunctionRep * fr = m_func_reps[i];
531  if ( fr -> isComposite () ) {
532  CompositeFunctionRep * composite
533  = dynamic_cast < CompositeFunctionRep * > ( fr );
534  bool yes = composite -> isMember ( rep );
535  if ( yes ) {
536  rep = composite;
537  break;
538  }
539  }
540  }
541 
542  return rep;
543 }
544 
545 bool
547 fitFunction ( PlotterBase * plotter, FunctionRep * rep )
548 {
549  rep = getComposite ( plotter, rep );
550  bool ok = rep -> fitFunction ();
551  if ( rep -> isComposite () ) {
552  setErrorsFromComposite ( plotter, rep );
553  }
554 
555  return ok;
556 }
557 
558 bool
560 tryFitFunction ( PlotterBase * plotter, FunctionRep * func_rep )
561 {
562  saveParameters ( plotter );
563  double chi_sq = func_rep -> objectiveValue ();
564  bool ok = fitFunction ( plotter, func_rep );
565  if ( ok ) {
566  double new_chi_sq = func_rep -> objectiveValue ();
567  if ( new_chi_sq > chi_sq ) {
568  restoreParameters ( plotter );
569  }
570  } else {
571  restoreParameters ( plotter );
572  }
573 
574  return ok;
575 }
576 
577 void
579 setFitRange ( PlotterBase * plotter, const Range & range )
580 {
581  if ( hasFunction ( plotter, 0 ) ) {
582  FunctionRep * frep = getFunctionRep ( plotter );
583  frep -> setCutRange ( range );
584  plotter -> update ();
585  }
586 }
587 
588 void
590 setFitRange ( PlotterBase * plotter, double low, double high )
591 {
592  const Range range ( low, high );
593  setFitRange ( plotter, range );
594 }
595 
597 {
598  findFunctions ( plotter );
599  vector< FunctionRep * >::iterator it = m_func_reps.begin ();
600 
601  for ( ; it != m_func_reps.end (); ++it ) {
602  FunctionRep * frep = *it;
603  FunctionBase * function = frep->getFunction ();
604  if ( function->isComposite () ) continue;
605  frep->saveParameters ();
606  }
607 }
608 
610 {
611  findFunctions ( plotter );
612  vector< FunctionRep * >::iterator it = m_func_reps.begin ();
613 
614  for ( ; it != m_func_reps.end (); ++it ) {
615  FunctionRep * frep = *it;
616  FunctionBase * function = frep->getFunction ();
617  if ( function->isComposite () ) continue;
618  frep->restoreParameters ();
619  }
620 
621 }
622 
623 FunctionRep *
625 getFunctionRep ( const PlotterBase * plotter ) const
626 {
627  FunctionRep * frep = 0;
628  DataRep * rep = 0;
629  int index = plotter->activePlotIndex ();
630 
631  if ( index >= 0 ) {
632  rep = plotter -> getDataRep ( index );
633  }
634  else {
635  rep = plotter -> getTarget ();
636  }
637 
638  FunctionRep * test = dynamic_cast < FunctionRep * > ( rep );
639  if ( test != 0 ) { // if ctrl-clicked, could be the function
640  frep = test;
641  }
642  else {
643  frep = getFunctionRep ( plotter, rep );
644  }
645 
646  return frep;
647 }
648 
652 FunctionRep *
654 getFunctionRep ( const PlotterBase * plotter, const DataRep * datarep ) const
655 {
656  FunctionRep * frep = 0;
657 
658  fillFunctionReps ( m_func_reps, plotter, datarep );
659 
660  for ( unsigned int i = 0; i < m_func_reps.size(); i++ ) {
661  FunctionRep * rep = m_func_reps[i];
662  const DataRep * drep = rep -> getTarget ();
663  if ( drep != datarep ) continue;
664  frep = rep;
665  if ( frep -> isComposite () ) break; // use composite if found
666  }
667 
668  return frep;
669 }
670 
672 createFuncView ( const ViewFactory * factory,
673  PlotterBase * plotter,
674  const std::string & name )
675 {
676  FunctionRep * frep = getFunctionRep ( plotter );
677  assert ( frep != 0 );
678 
680  string nullstring ("");
681 
682  return controller->createTextView ( factory, frep, name, nullstring );
683 }
684 
685 bool
687 setFitter ( const PlotterBase * plotter, const std::string & name )
688 {
689  FunctionRep * frep = getFunctionRep ( plotter );
690  bool ok = false;
691  if ( frep != 0 ) {
692  ok = setFitter ( frep, name );
693  }
694 
695  return ok;
696 }
697 
698 const string &
700 getFitterName ( const PlotterBase * plotter )
701 {
702  FunctionRep * frep = getFunctionRep ( plotter );
703  Fitter * fitter = frep -> getFitter ();
704 
705  return fitter -> name ();
706 }
707 
708 Fitter *
710 getFitter ( const PlotterBase * plotter )
711 {
712  FunctionRep * frep = getFunctionRep ( plotter );
713  Fitter * fitter = frep -> getFitter ();
714 
715  return fitter;
716 }
717 
718 
719 bool
721 setFitter ( FunctionRep * frep, const std::string & name )
722 {
723  FitterFactory * factory = FitterFactory::instance ();
724  Fitter * fitter = factory -> create ( name );
725 
726  return frep -> setFitter ( fitter );
727 }
728 
729 void
731 setDefaultFitter ( const std::string & name )
732 {
733  FitterFactory * factory = FitterFactory::instance ();
734  factory -> setDefault ( name );
735 }
736 
737 bool
739 changeFitter ( const PlotterBase * plotter,
740  const DataRep * drep,
741  const std::string & name )
742 {
743  FitterFactory * factory = FitterFactory::instance ();
744  Fitter * proto = factory -> prototype ( name );
745 
746  FunctionRep * frep = getFunctionRep ( plotter, drep );
747  const FunctionBase * func = frep -> getFunction ();
748  bool yes = proto -> isCompatible ( func );
749  if ( yes ) {
750  const Fitter * old_fitter = frep -> getFitter ();
751  Fitter * new_fitter = factory -> create ( name );
752  new_fitter -> copyFrom ( old_fitter );
753  frep -> setFitter ( new_fitter );
754  }
755  return yes;
756 }
757 
758 double
760 getObjectiveValue ( const PlotterBase * plotter, const DataRep * datarep )
761 {
762  FunctionRep * rep = getFunctionRep ( plotter, datarep );
763 
764  return rep -> objectiveValue ();
765 }
766 
767 const vector < vector < double > > &
769 getCovarianceMatrix ( const PlotterBase * plotter )
770 {
771  FunctionRep * rep = getFunctionRep ( plotter );
772 
773  return rep -> covarianceMatrix ();
774 }
775 
780 double
782 getChiSquared ( const PlotterBase * plotter )
783 {
784  const DataRep * rep = plotter -> getTarget ();
785 
786  return getObjectiveValue ( plotter, rep );
787 }
788 
789 int
791 getDegreesOfFreedom ( const PlotterBase * plotter )
792 {
793  FunctionRep * rep = getFunctionRep ( plotter );
794 
795  return rep -> degreesOfFreedom ();
796 }
797 
798 NTuple *
800 createNTuple ( const PlotterBase * plotter, const FunctionRep * rep )
801 {
802  NTuple * ntuple = 0;
803  if ( rep == 0 ) {
804  int old_index = plotter -> activePlotIndex ();
805  int index = getUniqueNonFunctionIndex ( plotter );
806  if ( index < 0 ) return 0;
807 
808  int saved_index = plotter -> activePlotIndex ();
809  PlotterBase * pb = const_cast < PlotterBase * > ( plotter );
810  pb -> setActivePlot ( index, false ); // no redraw
811  ntuple = plotter -> createNTuple ();
812  pb -> setActivePlot ( saved_index, false );
813  rep = getFunctionRep ( plotter );
814  pb -> setActivePlot ( old_index, true );
815  }
816  else {
817  const DataRep * target = rep -> getTarget ();
818  ntuple = target -> createNTuple ();
819  }
820  if ( rep != 0 ) {
821  const FunctionBase * function = rep -> getFunction ();
822  vector < double > & x = ntuple -> getColumn ( 0 ); // X coordinate
823  vector < double > & y = ntuple -> getColumn ( 1 ); // Y coordinate
824  unsigned int size = ntuple -> rows ();
825  vector < double > values ( size );
826  vector < double > residuals ( size );
827 
828  for ( unsigned int i = 0; i < x.size(); i++ ) {
829  values [i] = function -> operator () ( x[i] );
830  residuals [i] = y[i] - values [i];
831  }
832 
833  ntuple -> addColumn ( function->name (), values );
834  ntuple -> addColumn ( "Residuals", residuals );
835  }
836 
837  return ntuple;
838 }
839 
840 PlotterBase *
842 createResidualsDisplay ( PlotterBase * plotter, const FunctionRep * frep )
843 {
844  NTuple * ntuple = createNTuple ( plotter, frep );
845  ntuple -> setTitle ( plotter -> getTitle () );
846  const string old_label ( "Residuals" );
847  DataSourceController::instance () -> registerNTuple ( ntuple );
848  vector < string > bindings ( 4 );
849  bindings[0] = ntuple -> getLabelAt ( 0 );
850  bindings[1] = old_label;
851  bindings[2] = "nil";
852  bindings[3] = ntuple -> getLabelAt ( 3 );
853  bool axis_scaled = plotter -> isAxisScaled ( Axes::Y );
854  if ( axis_scaled ) {
855  int i = ntuple -> indexOf ( old_label );
856  const string new_label ("Residuals (scaled)");
857  ntuple -> setLabelAt ( new_label , i );
858  bindings[1] = new_label;
859  }
861 
862 // Make scaling of x-axis match that of the original plot.
863  PlotterBase * new_plotter = controller -> createDisplay ( "XY Plot",
864  *ntuple,
865  bindings );
866  if ( axis_scaled ) {
867  double factor = plotter -> getScaleFactor ( Axes::Y );
868  new_plotter -> setScaleFactor ( Axes::Y, factor );
869  }
870  controller->setLog( new_plotter, "x",
871  controller -> getLog( plotter, "x" ) );
872 
873  return new_plotter;
874 }
875 
876 PlotterBase *
879 {
880  int n = 100;
881  double xmin, xmax, ymin, ymax;
882 
883  // Set up the NTuple
884  NTuple* ntuple = new( NTuple );
885  ellipsoidNTuple ( masterPlot, rep, ntuple, n, xmin, xmax, ymin, ymax );
886  ntuple -> setTitle ( masterPlot -> getTitle () );
887 
888  vector < string > bindings ( 1 );
889  bindings[0] = ntuple -> getLabelAt ( 0 );
890 
892 
893  // Z plot. X and Y coordinates are merely indices
894  PlotterBase * ellipsePlot = dcontroller ->
895  createDisplay ( "Image", *ntuple, bindings );
896 
897  // Rescale and translate the Z plot so that the plot
898  // becomes contained in the rectangle bound by xmin, xmax
899  // ymin and ymax.
900  ellipsePlot -> setBinWidth( "x", ( xmax - xmin ) / ( n - 1.0 ) );
901  ellipsePlot -> setBinWidth( "y", ( ymax - ymin ) / ( n - 1.0 ) );
902 
903  ellipsePlot -> setOffsets ( xmin, ymin );
904 
905  ellipsePlot -> setRange ( string( "X" ), xmin, xmax );
906  ellipsePlot -> setRange ( string( "Y" ), ymin, ymax );
907 
908  // Establish the relation between source masterPlot and
909  // this new plot. This relationship shall be exploited
910  // when we refresh the error plots.
911  int index = dcontroller -> activeDataRepIndex( masterPlot );
912  ellipsePlot -> setParentDataRepIndex( index );
913  ellipsePlot -> setParentPlotter( masterPlot );
914 
915  return ellipsePlot;
916 }
917 
918 void
921  NTuple* nt, int n,
922  double & xmin, double & xmax,
923  double & ymin, double & ymax )
924 {
925  string text ( "Confidence ellipse " );
926  text += String::convert ( ++m_confid_count );
927  nt -> setName ( text );
928 
930  controller -> registerNTuple ( nt );
931 
932  assert( frep );
933 
934  // Get projected covariance i.e. take the sub-matrix
935  // out of the covariance matrix which corresponds
936  // to the 2 parameters of interest whose correlation
937  // we would like to see.
938  vector< vector< double > > Sigma( 2 );
939  Sigma[0].resize( 2, 0.0 );
940  Sigma[1].resize( 2, 0.0 );
941 
942  const vector < vector < double > > & covariance
943  = frep -> covarianceMatrix ();
944 
945  Sigma[0][0] = covariance [m_x][m_x];
946  Sigma[0][1] = covariance [m_x][m_y];
947  Sigma[1][0] = covariance [m_y][m_x];
948  Sigma[1][1] = covariance [m_y][m_y];
949 
950  // Invert the projected covariance
951  vector< vector< double > > SigmaInv;
952  invertMatrix( Sigma, SigmaInv );
953 
954  // Decide the center of the ellipse to be parameter
955  // value of the 2 parameters of interest whose correlation
956  // we would like to see.
957  vector< double > xbar( 2 );
958 
959  const Fitter * fitter = getFitter ( masterPlot );
960  vector < double > free_parms;
961  fitter -> fillFreeParameters ( free_parms );
962  xbar[ 0 ] = free_parms [ m_x ];
963  xbar[ 1 ] = free_parms [ m_y ];
964 
965  // Get the bounding ellipse and its bounds. Bounding ellipse is the
966  // 99.99% confidence ellipsoid. For mu = 2 the delta chi-square turns from
967  // Numerical Recipes in C as 18.4. etc etc.
968  unsigned int mu = free_parms.size();
969  NTuple * boundingTuple = ellipse( xbar, SigmaInv, m_deltaXSq[ mu - 1 ] );
970 
971  xmin = boundingTuple -> minElement ( 0 );
972  xmax = boundingTuple -> maxElement ( 0 );
973  ymin = boundingTuple -> minElement ( 1 );
974  ymax = boundingTuple -> maxElement ( 1 );
975 
976  delete boundingTuple; // Served its purpose
977 
978  // Create the NTuple for the contour plot
979  // This NTuple is to be build column-wise
980  // keeping in view the way Z-plot works.
981  vector< double > p( n * n ), a( 2 );
982 
983  double dx = ( xmax - xmin ) / ( n - 1.0 );
984  double dy = ( ymax - ymin ) / ( n - 1.0 );
985  double delta = 0.0;
986 
987  for( int i = 0; i < n; i++ )
988  for( int j = 0; j < n; j++ )
989  {
990  a[ 0 ] = xmin + i * dx - xbar[ 0 ] ;
991  a[ 1 ] = ymin + j * dy - xbar[ 1 ];
992  delta = quadraticProduct( SigmaInv, a );
993  p[ n * i + j ] = 1 - gammq( mu / 2.0, delta / 2.0 );
994  }
995 
996  nt -> clear();
997  nt -> addColumn( "Percent", 100 * p );
998 
999  return;
1000 }
1001 
1002 PlotterBase *
1005 {
1006  int n = 100;
1007  double xmin, xmax, ymin, ymax;
1008 
1009  PlotterBase* masterPlot = plot->getParentPlotter();
1010 
1011  // Set up the NTuple
1012  NTuple* ntuple = new( NTuple );
1013  ellipsoidNTuple ( masterPlot, frep, ntuple, n, xmin, xmax, ymin, ymax );
1014  ntuple -> setTitle ( masterPlot -> getTitle () );
1015 
1016  // First get the selected DataRep from the plotter. The DataRep
1017  // will have a projector that is a NTupleProjector. It doesn't
1018  // know, but we do, so we downcast and set the NTuple.
1019  DataRep * drep = plot -> selectedDataRep();
1020  NTupleProjector * ntProjector =
1021  dynamic_cast < NTupleProjector * > ( drep -> getProjector() );
1022  ntProjector -> setNTuple ( ntuple );
1023 
1024  NTuple * nt = const_cast < NTuple * > ( ntuple );
1025  nt -> addObserver ( ntProjector );
1026 
1027  // Rescale and translate the Z plot so that the plot
1028  // becomes contained in the rectangle bound by xmin, xmax
1029  // ymin and ymax.
1030  plot -> setBinWidth( "x", ( xmax - xmin ) / ( n - 1.0 ) );
1031  plot -> setBinWidth( "y", ( ymax - ymin ) / ( n - 1.0 ) );
1032 
1033  plot -> setOffsets ( xmin, ymin );
1034 
1035  plot -> setRange ( string( "X" ), xmin, xmax );
1036  plot -> setRange ( string( "Y" ), ymin, ymax );
1037 
1038  return plot;
1039 }
1040 
1041 NTuple *
1043 ellipse ( const std::vector< double > & xbar,
1044  std::vector< std::vector < double > > & SigmaInv,
1045  double Csquare )
1046 {
1047  // First argument should be a 2 D vector, the center of the ellipse //
1048  assert( xbar.size() == 2 );
1049 
1050  // Second argument should be a 2 x 2 SPD matrix i.e. two rows two cols //
1051  assert( SigmaInv.size() == 2 );
1052  assert( SigmaInv[0].size() == 2 || SigmaInv[1].size() == 2 );
1053 
1054  // Second argument should be a 2 x 2 -- Symmetric -- PD matrix
1055  // Under numerical round-offs it might not be true so we take mean of the
1056  // entries.
1057  double temp = ( SigmaInv[0][1] + SigmaInv[1][0] ) / 2;
1058  SigmaInv[0][1] = temp;
1059  SigmaInv[1][0] = temp;
1060 
1061  // Third argument should be a positive scalar "c^2", Square of "Radius" //
1062  assert( Csquare > DBL_EPSILON );
1063 
1064  // The eigenvalues of the SigmaInv matrix
1065  double b = -( SigmaInv[0][0] + SigmaInv[1][1] );
1066  double c = SigmaInv[0][0] * SigmaInv[1][1] - SigmaInv[0][1] * SigmaInv[1][0];
1067 
1068  double lambda1 = ( -b + sqrt( b * b - 4 * c ) ) / 2;
1069  double lambda2 = ( -b - sqrt( b * b - 4 * c ) ) / 2;
1070 
1071  // Determining the angles the ellipse axis make with X-axis
1072  double alpha1 = atan( (SigmaInv[0][0] - lambda1) / SigmaInv[0][1] );
1073  double a = sqrt( Csquare / lambda1 );
1074  b = sqrt( Csquare / lambda2 );
1075 
1076  // Creating a rotation matrix
1077  double Rot00 = cos( alpha1 );
1078  double Rot11 = cos( alpha1 );
1079  double Rot01 = -sin( alpha1 );
1080  double Rot10 = sin( alpha1 );
1081 
1082  // Generating N-points on the ellipse
1083  int N = 100;
1084  vector< double > x, y;
1085 
1086  x.resize( N, 0.0 );
1087  y.resize( N, 0.0 );
1088 
1089  for( int i = 0; i < N; i++)
1090  {
1091  // Unrotated untranslated point
1092  double theta = ( 2 * M_PI * i ) / ( (double) (N-1) );
1093  double x0 = a * cos( theta );
1094  double x1 = b * sin( theta );
1095 
1096  //x = Rot * x; i.e. force in a rotation
1097  double xrot0 = Rot00 * x0 + Rot01 * x1;
1098  double xrot1 = Rot10 * x0 + Rot11 * x1;
1099 
1100  //x = x + xbar; i.e. force in a translation
1101  x[i] = xrot0 + xbar[0];
1102  y[i] = xrot1 + xbar[1];
1103  }
1104 
1105  // Create an NTuple out of the above vector
1106  NTuple* ntuple = new( NTuple );
1107  ntuple -> addColumn( "X", x );
1108  ntuple -> addColumn( "Y", y );
1109 
1110  return ntuple;
1111 }
1112 
1114  int index )
1115 {
1116  assert( axes == Axes::X || axes == Axes::Y );
1117 
1118  if( axes == Axes::X )
1119  m_x = index;
1120 
1121  if( axes == Axes::Y )
1122  m_y = index;
1123 
1124  return EXIT_SUCCESS;
1125 }
1126 
1127 bool
1129 functionExists ( const std::string & name )
1130 {
1132 
1133  return factory -> exists ( name );
1134 }
1135 
1136 bool
1138 isCompatible ( const std::string & function,
1139  const std::string & fitter )
1140 {
1141  bool yes = true;
1142 
1144  FunctionBase * fun_proto = fun_fac -> prototype ( function );
1145 
1146  FitterFactory * fit_fac = FitterFactory::instance ();
1147  Fitter * fit_proto = fit_fac -> prototype ( fitter );
1148 
1149  yes = fit_proto -> isCompatible ( fun_proto );
1150 
1151  return yes;
1152 }
1153 
1154 const vector < string > &
1157 {
1158  return FunctionFactory::instance () -> names ();
1159 }

Generated for HippoDraw Class Library by doxygen