13 #include "msdevstudio/MSconfig.h"
51 #ifdef ITERATOR_MEMBER_DEFECT
70 using namespace hippodraw;
72 using namespace Numeric;
76 FunctionController::FunctionController ( )
104 const vector < string > &
110 return factory -> names ();
119 return factory -> getDefault ();
129 for (
int i = 0; i < number; i++ ) {
132 if ( fr != 0 )
continue;
136 if ( count != 1 ) index = -1;
148 for (
int i = 0; i < number; i++ ) {
151 if ( frep == 0 )
continue;
163 if (
function == 0 ) {
164 string what (
"FunctionController: Unable to create function of type `" );
167 throw std::runtime_error ( what );
180 if (
function -> isComposite () ) {
184 unsigned int dims =
function -> dimensions ();
195 const string & name = factory -> getDefault ();
213 DataRep * rep = plotter -> getDataRep ( 0 );
219 return frep -> getFunction ();
228 const std::string & name,
233 func_rep -> initializeWith ( rep );
251 func_rep -> initializeWith ( drep );
272 if ( frep -> isComposite () ==
false ) {
273 DataRep * target = frep -> getTarget ();
284 rep -> setCutRange (
true );
291 data_rep -> addCut ();
293 if ( func_rep != 0 ) {
303 if ( func_rep != 0 ) {
304 func_rep -> removeCut ();
322 const string & name = factory -> getDefault ();
326 if ( frep -> isComposite () ) {
327 frep -> addToComposite ( func_rep );
331 composite -> initializeWith ( drep );
335 composite -> addToComposite ( frep );
336 composite -> addToComposite ( func_rep );
337 return_rep = composite;
342 if ( return_rep == 0 ) return_rep = func_rep;
352 if ( frep -> isComposite () ) {
355 const vector < FunctionRep * > & freps = composite -> getFunctionReps ();
356 unsigned int size = freps.size();
359 for (
int i = size -1; i >= 0; i-- ) {
361 plotter -> removeDataRep ( rep );
362 composite -> removeFromComposite ( rep );
365 plotter -> removeDataRep ( frep );
370 vector < FunctionRep * > :: iterator it
376 plotter -> removeDataRep ( frep );
377 for (
unsigned int i = 0; i <
m_func_reps.size (); i++ ) {
379 if ( rep -> isComposite () ) {
380 rep -> removeFromComposite ( frep );
385 vector < FunctionRep * >::iterator iter =
m_func_reps.begin ();
391 if ( rep -> isComposite () ) {
394 if ( composite -> count() < 2 ) {
397 = composite -> getFunctionReps ();
399 crep -> setInComposite (
false );
420 assert ( plotter != 0 );
428 const vector< string > &
435 vector< FunctionRep * >:: const_iterator it =
m_func_reps.begin ();
439 if ( rep != 0 && fp->
getTarget() != rep )
continue;
441 const string & name =
function->
name ();
455 int number = plotter -> getNumDataReps ();
456 for (
int i = 0; i < number; i++ ) {
457 DataRep * rep = plotter -> getDataRep ( i );
461 if ( frep -> getTarget () == drep ) {
462 freps.push_back ( frep );
466 freps.push_back ( frep );
481 for (
unsigned int i = 0; i <
size; i++ ) {
483 if ( rep -> isInComposite () ==
false ) {
484 freps.push_back ( rep );
494 const vector < double > & errors = composite -> principleErrors();
495 vector < double >::const_iterator begin = errors.begin();
497 DataRep * target = composite -> getTarget ();
500 vector < FunctionRep * >::iterator first =
m_func_reps.begin();
504 if ( rep -> isComposite() )
continue;
505 const vector < double > & t = rep -> principleErrors ();
506 unsigned int number = t.size();
527 const DataRep * target = rep -> getTarget ();
529 for (
unsigned int i = 0; i <
m_func_reps.size (); i++ ) {
531 if ( fr -> isComposite () ) {
534 bool yes = composite -> isMember ( rep );
551 if ( rep -> isComposite () ) {
563 double chi_sq = func_rep -> objectiveValue ();
566 double new_chi_sq = func_rep -> objectiveValue ();
567 if ( new_chi_sq > chi_sq ) {
583 frep -> setCutRange ( range );
584 plotter -> update ();
592 const Range range ( low, high );
599 vector< FunctionRep * >::iterator it =
m_func_reps.begin ();
604 if ( function->isComposite () )
continue;
612 vector< FunctionRep * >::iterator it =
m_func_reps.begin ();
617 if ( function->isComposite () )
continue;
632 rep = plotter -> getDataRep ( index );
635 rep = plotter -> getTarget ();
660 for (
unsigned int i = 0; i <
m_func_reps.size(); i++ ) {
662 const DataRep * drep = rep -> getTarget ();
663 if ( drep != datarep )
continue;
665 if ( frep -> isComposite () )
break;
674 const std::string & name )
677 assert ( frep != 0 );
680 string nullstring (
"");
682 return controller->
createTextView ( factory, frep, name, nullstring );
705 return fitter -> name ();
724 Fitter * fitter = factory -> create ( name );
734 factory -> setDefault ( name );
741 const std::string & name )
744 Fitter * proto = factory -> prototype ( name );
751 Fitter * new_fitter = factory -> create ( name );
752 new_fitter -> copyFrom ( old_fitter );
764 return rep -> objectiveValue ();
767 const vector < vector < double > > &
773 return rep -> covarianceMatrix ();
784 const DataRep * rep = plotter -> getTarget ();
795 return rep -> degreesOfFreedom ();
804 int old_index = plotter -> activePlotIndex ();
806 if ( index < 0 )
return 0;
808 int saved_index = plotter -> activePlotIndex ();
810 pb -> setActivePlot ( index,
false );
812 pb -> setActivePlot ( saved_index,
false );
814 pb -> setActivePlot ( old_index,
true );
817 const DataRep * target = rep -> getTarget ();
822 vector < double > & x = ntuple -> getColumn ( 0 );
823 vector < double > & y = ntuple -> getColumn ( 1 );
824 unsigned int size = ntuple -> rows ();
825 vector < double > values ( size );
826 vector < double > residuals ( size );
828 for (
unsigned int i = 0; i < x.size(); i++ ) {
829 values [i] =
function -> operator () ( x[i] );
830 residuals [i] = y[i] - values [i];
833 ntuple -> addColumn ( function->name (), values );
834 ntuple -> addColumn (
"Residuals", residuals );
845 ntuple -> setTitle ( plotter -> getTitle () );
846 const string old_label (
"Residuals" );
848 vector < string > bindings ( 4 );
849 bindings[0] = ntuple -> getLabelAt ( 0 );
850 bindings[1] = old_label;
852 bindings[3] = ntuple -> getLabelAt ( 3 );
853 bool axis_scaled = plotter -> isAxisScaled (
Axes::Y );
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;
863 PlotterBase * new_plotter = controller -> createDisplay (
"XY Plot",
867 double factor = plotter -> getScaleFactor (
Axes::Y );
868 new_plotter -> setScaleFactor (
Axes::Y, factor );
870 controller->
setLog( new_plotter,
"x",
871 controller -> getLog( plotter,
"x" ) );
881 double xmin, xmax, ymin, ymax;
885 ellipsoidNTuple ( masterPlot, rep, ntuple, n, xmin, xmax, ymin, ymax );
886 ntuple -> setTitle ( masterPlot -> getTitle () );
888 vector < string > bindings ( 1 );
889 bindings[0] = ntuple -> getLabelAt ( 0 );
895 createDisplay (
"Image", *ntuple, bindings );
900 ellipsePlot -> setBinWidth(
"x", ( xmax - xmin ) / ( n - 1.0 ) );
901 ellipsePlot -> setBinWidth(
"y", ( ymax - ymin ) / ( n - 1.0 ) );
903 ellipsePlot -> setOffsets ( xmin, ymin );
905 ellipsePlot -> setRange (
string(
"X" ), xmin, xmax );
906 ellipsePlot -> setRange (
string(
"Y" ), ymin, ymax );
911 int index = dcontroller -> activeDataRepIndex( masterPlot );
912 ellipsePlot -> setParentDataRepIndex( index );
913 ellipsePlot -> setParentPlotter( masterPlot );
922 double & xmin,
double & xmax,
923 double & ymin,
double & ymax )
925 string text (
"Confidence ellipse " );
927 nt -> setName ( text );
930 controller -> registerNTuple ( nt );
938 vector< vector< double > > Sigma( 2 );
939 Sigma[0].resize( 2, 0.0 );
940 Sigma[1].resize( 2, 0.0 );
942 const vector < vector < double > > & covariance
943 = frep -> covarianceMatrix ();
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];
951 vector< vector< double > > SigmaInv;
957 vector< double > xbar( 2 );
960 vector < double > free_parms;
961 fitter -> fillFreeParameters ( free_parms );
962 xbar[ 0 ] = free_parms [
m_x ];
963 xbar[ 1 ] = free_parms [
m_y ];
968 unsigned int mu = free_parms.size();
971 xmin = boundingTuple -> minElement ( 0 );
972 xmax = boundingTuple -> maxElement ( 0 );
973 ymin = boundingTuple -> minElement ( 1 );
974 ymax = boundingTuple -> maxElement ( 1 );
976 delete boundingTuple;
981 vector< double > p( n * n ), a( 2 );
983 double dx = ( xmax - xmin ) / ( n - 1.0 );
984 double dy = ( ymax - ymin ) / ( n - 1.0 );
987 for(
int i = 0; i < n; i++ )
988 for(
int j = 0; j < n; j++ )
990 a[ 0 ] = xmin + i * dx - xbar[ 0 ] ;
991 a[ 1 ] = ymin + j * dy - xbar[ 1 ];
993 p[ n * i + j ] = 1 -
gammq( mu / 2.0, delta / 2.0 );
997 nt -> addColumn(
"Percent", 100 * p );
1007 double xmin, xmax, ymin, ymax;
1013 ellipsoidNTuple ( masterPlot, frep, ntuple, n, xmin, xmax, ymin, ymax );
1014 ntuple -> setTitle ( masterPlot -> getTitle () );
1019 DataRep * drep = plot -> selectedDataRep();
1022 ntProjector -> setNTuple ( ntuple );
1025 nt -> addObserver ( ntProjector );
1030 plot -> setBinWidth(
"x", ( xmax - xmin ) / ( n - 1.0 ) );
1031 plot -> setBinWidth(
"y", ( ymax - ymin ) / ( n - 1.0 ) );
1033 plot -> setOffsets ( xmin, ymin );
1035 plot -> setRange (
string(
"X" ), xmin, xmax );
1036 plot -> setRange (
string(
"Y" ), ymin, ymax );
1044 std::vector< std::vector < double > > & SigmaInv,
1048 assert( xbar.size() == 2 );
1051 assert( SigmaInv.size() == 2 );
1052 assert( SigmaInv[0].
size() == 2 || SigmaInv[1].
size() == 2 );
1057 double temp = ( SigmaInv[0][1] + SigmaInv[1][0] ) / 2;
1058 SigmaInv[0][1] = temp;
1059 SigmaInv[1][0] = temp;
1062 assert( Csquare > DBL_EPSILON );
1065 double b = -( SigmaInv[0][0] + SigmaInv[1][1] );
1066 double c = SigmaInv[0][0] * SigmaInv[1][1] - SigmaInv[0][1] * SigmaInv[1][0];
1068 double lambda1 = ( -b + sqrt( b * b - 4 * c ) ) / 2;
1069 double lambda2 = ( -b - sqrt( b * b - 4 * c ) ) / 2;
1072 double alpha1 = atan( (SigmaInv[0][0] - lambda1) / SigmaInv[0][1] );
1073 double a = sqrt( Csquare / lambda1 );
1074 b = sqrt( Csquare / lambda2 );
1077 double Rot00 = cos( alpha1 );
1078 double Rot11 = cos( alpha1 );
1079 double Rot01 = -sin( alpha1 );
1080 double Rot10 = sin( alpha1 );
1084 vector< double > x, y;
1089 for(
int i = 0; i < N; i++)
1092 double theta = ( 2 * M_PI * i ) / ( (
double) (N-1) );
1093 double x0 = a * cos( theta );
1094 double x1 = b * sin( theta );
1097 double xrot0 = Rot00 * x0 + Rot01 * x1;
1098 double xrot1 = Rot10 * x0 + Rot11 * x1;
1101 x[i] = xrot0 + xbar[0];
1102 y[i] = xrot1 + xbar[1];
1107 ntuple -> addColumn(
"X", x );
1108 ntuple -> addColumn(
"Y", y );
1124 return EXIT_SUCCESS;
1133 return factory -> exists ( name );
1139 const std::string & fitter )
1144 FunctionBase * fun_proto = fun_fac -> prototype (
function );
1147 Fitter * fit_proto = fit_fac -> prototype ( fitter );
1154 const vector < string > &