00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*
00008    Plugin to calculate a nonlinear regression for an input time series, 
00009    using a specified nonlinear model.
00010 
00011    File:     plug_nlfit.c
00012    Author:   B. Douglas Ward
00013    Date:     17 June 1997
00014 
00015    Mod:      Added external time reference capability and added option for 
00016              absolute noise parameter constraints.
00017    Date:     22 August 1997
00018 
00019    Mod:      Added options for percent signal change above baseline, and
00020              calculation of area under the signal above baseline.
00021    Date:     26 November 1997
00022 
00023    Mod:      Print error message if unable to locate any signal models or 
00024              any noise models.
00025    Date:     26 December 1997
00026 
00027    Mod:      Added novar flag to eliminate unnecessary calculations.
00028    Date:     13 July 1999
00029 
00030    Mod:      Adjust F-statistics if parameter constraints force a parameter
00031              to be a constant.
00032    Date:     08 February 2000
00033 
00034    Mod:      Changes for output of R^2 (coefficient of multiple determination),
00035              and standard deviation of residuals from full model fit.
00036              Added global variable calc_tstats.
00037              Also, added screen display of p-values.
00038    Date:     10 May 2000
00039 
00040 */
00041 
00042 /*---------------------------------------------------------------------------*/
00043 
00044 #define PROGRAM_NAME "plug_nlfit"                    /* name of this program */
00045 #define PROGRAM_AUTHOR "B. Douglas Ward"                   /* program author */
00046 #define PROGRAM_DATE "10 May 2000"               /* date of last program mod */
00047 
00048 /*---------------------------------------------------------------------------*/
00049 
00050 
00051 #include "afni.h"
00052 
00053 #ifndef ALLOW_PLUGINS
00054 #  error "Plugins not properly set up -- see machdep.h"
00055 #endif
00056 
00057 
00058 #include <math.h>
00059 #include <stdlib.h>
00060 
00061 #include "mrilib.h"
00062 
00063 #include "matrix.h"
00064 #include "simplex.h"
00065 #include "NLfit.h"
00066 
00067 #include "matrix.c"
00068 #include "simplex.c"
00069 #include "NLfit.c"
00070 
00071 /***** 22 July 1998 -- RWCox:
00072        Modified to allow DELT to be set from the TR of the input file *****/
00073 
00074 static float DELT = 1.0;   /* default */
00075 static int   inTR = 0 ;    /* set to 1 if -inTR option is used */
00076 static float dsTR = 0.0 ;  /* TR of the input file */
00077 
00078 
00079 /***********************************************************************
00080   Plugin to provide nonlinear least squares fitting 1D function for graphs
00081 ************************************************************************/
00082 
00083 /*------------- string to 'help' the user -------------*/
00084 
00085 static char helpstring[] =
00086    " Purpose: Control the 'NLfit' and 'NLerr' 1D functions.\n"
00087    "\n"
00088    " Control:    Ignore       = Number of points to ignore at start \n"
00089    "                            of each timeseries. \n"
00090    "             NRandom      = Number of random test points. \n"
00091    "             NBest        = Find opt. soln. for these best test points. \n"
00092    " \n"
00093    " Models:     Noise        = Label of noise model to use. \n"
00094    "             Signal       = Label of signal model to use. \n"
00095    "             Noise Constr = Relative or Absolute noise constraints. \n"
00096    " \n"
00097    " Noise:      Parameter    = Index of noise model parameter. \n"
00098    "             Min Constr   = Minimum value for noise model parameter. \n"
00099    "             Max Constr   = Maximum value for noise model parameter. \n"
00100    " \n"
00101    " Signal:     Parameter    = Index of signal model parameter.\n"
00102    "             Min Constr   = Minimum value for signal model parameter. \n"
00103    "             Max Constr   = Maximum value for signal model parameter. \n"
00104    " \n"
00105    " Time Scale: Reference    = Internal or External time reference. \n"
00106    "             File         = External time reference file name. \n" 
00107    "Author -- BD Ward"
00108 ;
00109 
00110 
00111 
00112 /*--------------- prototypes for internal routines ---------------*/
00113 
00114 char * NL_main( PLUGIN_interface * ) ;  /* the entry point */
00115 
00116 void NL_fitter( int nt, double to, double dt, float * vec, char ** label ) ;
00117 void NL_error ( int nt, double to, double dt, float * vec, char ** label ) ;
00118 void NL_worker( int nt, double dt, float * vec, int dofit, char ** label ) ;
00119 
00120 
00121 /*---------------- global data -------------------*/
00122 static PLUGIN_interface * global_plint = NULL ;
00123 static int initialize=1 ;
00124 
00125 
00126 /*------- Global data for models ------*/
00127 char * constr_types[2] = {"Relative", "Absolute"};   /* option labels */
00128 char * time_refs[3] = {"Internal", "External" , "-inTR" };
00129 int plug_ignore = 3;
00130 int plug_nrand = 100;
00131 int plug_nbest = 5;
00132 int plug_nabs = 0;
00133 int plug_timeref = 0;
00134 char plug_tfilename[MAX_NAME_LENGTH] = "";
00135 
00136 /*----- declare reduced (noise) model variables -----*/
00137 int num_noise_models;                /* number of noise models */
00138 int plug_noise_index;                /* index of noise model */
00139 char * noise_labels[MAX_MODELS];     /* names of noise models */
00140 vfp plug_nmodel[MAX_MODELS];     /* pointer to noise model */
00141 int plug_r[MAX_MODELS];          /* number of parameters in the noise model */
00142 char * noise_plabels[MAX_MODELS][MAX_PARAMETERS];
00143 float plug_min_nconstr[MAX_MODELS][MAX_PARAMETERS]; 
00144                            /* minimum parameter constraints for noise model */
00145 float plug_max_nconstr[MAX_MODELS][MAX_PARAMETERS]; 
00146                            /* maximum parameter constraints for noise model */
00147 
00148 /*----- declare full (signal+noise) model variables -----*/
00149 int num_signal_models;               /* number of signal models */
00150 int plug_signal_index;               /* index of signal model */
00151 char * signal_labels[MAX_MODELS];    /* names of signal models */
00152 vfp plug_smodel[MAX_MODELS];     /* pointer to signal model */
00153 int plug_p[MAX_MODELS];          /* number of parameters in the signal model */
00154 char * signal_plabels[MAX_MODELS][MAX_PARAMETERS];
00155 float plug_min_sconstr[MAX_MODELS][MAX_PARAMETERS];  
00156                            /* minimum parameter constraints for signal model */
00157 float plug_max_sconstr[MAX_MODELS][MAX_PARAMETERS];  
00158                            /* maximum parameter constraints for signal model */
00159 
00160 
00161 /*---------------------------------------------------------------------------*/
00162 /*
00163   Routine to initialize the input options.
00164 */
00165  
00166 void initialize_options 
00167 (
00168   int * im1,               /* index of 1st image in time series for analysis */
00169   char ** nname,           /* noise model name */
00170   char ** sname,           /* signal model name */
00171   vfp * nmodel,            /* pointer to noise model */
00172   vfp * smodel,            /* pointer to signal model */  
00173   int * r,                 /* number of parameters in the noise model */
00174   int * p,                 /* number of parameters in the signal model */
00175   char *** npname,         /* noise parameter names */
00176   char *** spname,         /* signal parameter names */
00177   float ** min_nconstr,    /* minimum parameter constraints for noise model */
00178   float ** max_nconstr,    /* maximum parameter constraints for noise model */
00179   float ** min_sconstr,    /* minimum parameter constraints for signal model */
00180   float ** max_sconstr,    /* maximum parameter constraints for signal model */
00181   int * nabs,              /* use absolute constraints for noise parameters */
00182   int * nrand,             /* number of random vectors to generate */
00183   int * nbest,             /* number of random vectors to keep */
00184   float * rms_min,         /* minimum rms error to reject reduced model */
00185   char ** tfilename        /* file name for time point series */  
00186 )
00187  
00188 {
00189   int ip;                           /* parameter index */
00190   int ok;                           /* boolean for specified model exists */
00191   char message[MAX_NAME_LENGTH];    /* error message */
00192 
00193 
00194   *im1 = 1;
00195   *nrand = plug_nrand;
00196   *nbest = plug_nbest; 
00197   *nabs = plug_nabs;
00198   *rms_min = 0.0;
00199   *tfilename = plug_tfilename;
00200 
00201   *nname = noise_labels[plug_noise_index];
00202   *sname = signal_labels[plug_signal_index];
00203 
00204   *nmodel = plug_nmodel[plug_noise_index];
00205   *smodel = plug_smodel[plug_signal_index];
00206 
00207   *r = plug_r[plug_noise_index];
00208   *p = plug_p[plug_signal_index];
00209 
00210   *npname = noise_plabels[plug_noise_index];
00211   *spname = signal_plabels[plug_signal_index];
00212 
00213   /*----- allocate memory for parameter constraints -----*/
00214   *min_nconstr = (float *) malloc (sizeof(float) * (*r));
00215   if (*min_nconstr == NULL)  
00216     NLfit_error ("Unable to allocate memory for min_nconstr");
00217   *max_nconstr = (float *) malloc (sizeof(float) * (*r));
00218   if (*max_nconstr == NULL)
00219     NLfit_error ("Unable to allocate memory for max_nconstr");
00220   *min_sconstr = (float *) malloc (sizeof(float) * (*p));
00221   if (*min_sconstr == NULL)  
00222     NLfit_error ("Unable to allocate memory for min_sconstr");
00223   *max_sconstr = (float *) malloc (sizeof(float) * (*p));
00224   if (*max_sconstr == NULL)
00225     NLfit_error ("Unable to allocate memory for max_sconstr");
00226   
00227   /*----- initialize constraints -----*/
00228   for (ip = 0;  ip < (*r);  ip++)
00229     {
00230       (*min_nconstr)[ip] = plug_min_nconstr[plug_noise_index][ip];
00231       (*max_nconstr)[ip] = plug_max_nconstr[plug_noise_index][ip];      
00232     }
00233   
00234   for (ip = 0;  ip < (*p);  ip++)
00235     {
00236       (*min_sconstr)[ip] = plug_min_sconstr[plug_signal_index][ip];
00237       (*max_sconstr)[ip] = plug_max_sconstr[plug_signal_index][ip];      
00238     }
00239  
00240 }
00241 
00242 
00243 /*---------------------------------------------------------------------------*/
00244 /*
00245   Routine to check for valid inputs.
00246 */
00247   
00248 void check_for_valid_inputs ()
00249 {
00250 }
00251 
00252 
00253 /*---------------------------------------------------------------------------*/
00254 
00255 void initialize_program 
00256 (
00257   int * im1,               /* index of 1st image in time series for analysis */
00258   char ** nname,           /* noise model name */
00259   char ** sname,           /* signal model name */
00260   vfp * nmodel,            /* pointer to noise model */
00261   vfp * smodel,            /* pointer to signal model */  
00262   int * r,                 /* number of parameters in the noise model */
00263   int * p,                 /* number of parameters in the signal model */
00264   char *** npname,         /* noise parameter names */
00265   char *** spname,         /* signal parameter names */
00266   float ** min_nconstr,    /* minimum parameter constraints for noise model */
00267   float ** max_nconstr,    /* maximum parameter constraints for noise model */
00268   float ** min_sconstr,    /* minimum parameter constraints for signal model */
00269   float ** max_sconstr,    /* maximum parameter constraints for signal model */
00270   int * nabs,              /* use absolute constraints for noise parameters */
00271   int * nrand,             /* number of random vectors to generate */
00272   int * nbest,             /* number of random vectors to keep */
00273   float * rms_min,         /* minimum rms error to reject reduced model */
00274 
00275   float ** par_rdcd,       /* estimated parameters for the reduced model */
00276   float ** par_full,       /* estimated parameters for the full model */
00277   float ** tpar_full,      /* t-statistic of parameters in the full model */
00278 
00279   int ts_length,           /* length of time series data */  
00280   char ** tfilename,       /* file name for time point series */  
00281   float *** x_array,       /* independent variable matrix */
00282 
00283   float ** fit
00284 )
00285 
00286 {
00287   int dimension;           /* dimension of full model */
00288   int ip;                  /* parameter index */
00289   int it;                  /* time index */
00290   MRI_IMAGE * im, * flim;  /* pointers to image structures 
00291                               -- used to read 1D ASCII */
00292   int nt;                  /* number of points in 1D x data file */
00293   float * tar;
00294 
00295 
00296   /*----- intialize options -----*/
00297   initialize_options (im1, nname, sname, nmodel, smodel, r, p, npname, spname, 
00298                       min_nconstr, max_nconstr, min_sconstr, max_sconstr, 
00299                       nabs, nrand, nbest, rms_min, tfilename);
00300 
00301   /*----- check for valid inputs -----*/
00302   check_for_valid_inputs ();
00303 
00304 
00305   /*----- allocate space for independent variable matrix -----*/
00306   *x_array = (float **) malloc (sizeof(float *) * ts_length);
00307   if (*x_array == NULL)
00308     NLfit_error ("Unable to allocate memory for x_array");
00309   for (it = 0;  it < ts_length;  it++)
00310     {
00311       (*x_array)[it] = (float *) malloc (sizeof(float) * 3);
00312       if ((*x_array)[it] == NULL)
00313         NLfit_error ("Unable to allocate memory for x_array[it]");
00314     }
00315     
00316   /*----- initialize independent variable matrix -----*/
00317   if (!plug_timeref)
00318     {
00319       static float old_DELT = -1.0 ;
00320       DELT = (inTR && dsTR > 0.0) ? dsTR : 1.0 ;  /* 22 July 1998 */
00321       if( DELT != old_DELT ){
00322          old_DELT = DELT ;
00323          printf("NLfit: switch to TR = %g\n",DELT) ;
00324       }
00325 
00326       for (it = 0;  it < ts_length;  it++)  
00327         {
00328           (*x_array)[it][0] = 1.0;
00329           (*x_array)[it][1] = it * DELT;
00330           (*x_array)[it][2] = (it * DELT) * (it * DELT);
00331         }
00332     }
00333   else 
00334     {
00335         flim = mri_read_1D (*tfilename); 
00336         if (flim == NULL)
00337           NLfit_error ("Unable to read time reference file \n");
00338         nt = flim -> nx;
00339         if (nt < ts_length)
00340             NLfit_error ("Time reference array is too short");  
00341         tar = MRI_FLOAT_PTR(flim) ;
00342         for (it = 0;  it < ts_length;  it++)  
00343          {
00344            (*x_array)[it][0] = 1.0;
00345            (*x_array)[it][1] = tar[it] ;
00346            (*x_array)[it][2] = tar[it] * tar[it];
00347          }
00348         mri_free (flim);
00349      }
00350   
00351   dimension = (*r) + (*p);
00352 
00353   /*----- allocate memory space -----*/
00354   *par_rdcd = (float *) malloc (sizeof(float) * dimension);
00355   if (*par_rdcd == NULL)
00356     NLfit_error ("Unable to allocate memory for par_rdcd");
00357   *par_full = (float *) malloc (sizeof(float) * dimension);
00358   if (*par_full == NULL)
00359     NLfit_error ("Unable to allocate memory for par_full");
00360   *tpar_full = (float *) malloc (sizeof(float) * dimension);
00361   if (*tpar_full == NULL)
00362     NLfit_error ("Unable to allocate memory for tpar_full");
00363   *fit = (float *) malloc (sizeof(float) * (ts_length));
00364   if (*fit == NULL)
00365     NLfit_error ("Unable to allocate memory for fit");
00366 
00367 }
00368 
00369 
00370 /*---------------------------------------------------------------------------*/
00371 /*
00372   Release all allocated memory space.
00373 */
00374 
00375 void terminate_program 
00376 (
00377   int r,                       /* number of parameters in the noise model */
00378   int p,                       /* number of parameters in the signal model */
00379   int ts_length,               /* length of time series data */
00380   float *** x_array,           /* independent variable matrix */
00381   float ** par_rdcd,      /* estimated parameters for the reduced model */
00382   float ** min_nconstr,   /* min parameter constraints for noise model */
00383   float ** max_nconstr,   /* max parameter constraints for noise model */
00384   float ** par_full,      /* estimated parameters for the full model */
00385   float ** tpar_full,     /* t-statistic of parameters in full model */
00386   float ** min_sconstr,   /* min parameter constraints for signal model */
00387   float ** max_sconstr    /* max parameter constraints for signal model */
00388 )
00389  
00390 {
00391   int ip;                        /* parameter index */
00392   int it;                        /* time index */
00393 
00394 
00395   /*----- deallocate memory for parameter constraints -----*/
00396   if (*min_nconstr != NULL)  { free (*min_nconstr);  *min_nconstr = NULL; }
00397   if (*max_nconstr != NULL)  { free (*max_nconstr);  *max_nconstr = NULL; }
00398   if (*min_sconstr != NULL)  { free (*min_sconstr);  *min_sconstr = NULL; }
00399   if (*max_sconstr != NULL)  { free (*max_sconstr);  *max_sconstr = NULL; }
00400 
00401 
00402   /*----- deallocate space for independent variable matrix -----*/
00403   for (it = 0;  it < ts_length;  it++)
00404     if ((*x_array)[it] != NULL)
00405       { free ((*x_array)[it]);  (*x_array)[it] = NULL; }
00406   if (*x_array != NULL)     { free (*x_array);  *x_array = NULL; }
00407 
00408 
00409   /*----- deallocate space for parameters -----*/
00410   if (*par_rdcd != NULL)    { free (*par_rdcd);    *par_rdcd = NULL; }
00411   if (*par_full != NULL)    { free (*par_full);    *par_full = NULL; }
00412   if (*tpar_full != NULL)   { free (*tpar_full);   *tpar_full = NULL; }
00413 
00414 }
00415 
00416 
00417 /*---------------------------------------------------------------------------*/
00418 
00419 float *  nlfit
00420 (
00421   int ts_length,                     /* length of time series data */
00422   float * ts_array,                  /* input time series array */
00423   char ** label                      /* label output for this voxel */
00424 )
00425 
00426 {
00427   float * fit;                       /* nonlinear fit of time series data */
00428 
00429   /*----- declare input option variables -----*/
00430   int nabs;                /* use absolute constraints for noise parameters */
00431   int  nrand;              /* number of random vectors to generate */
00432   int  nbest;              /* number of random vectors to keep */
00433   float rms_min;           /* minimum rms error to reject reduced model */
00434 
00435   /*----- declare time series variables -----*/
00436   int im1;                 /* index of 1st image in time series for analysis */
00437   float ** x_array = NULL;     /* independent variable matrix */
00438   char * tfilename = NULL;     /* file name of time points */
00439 
00440   /*----- declare reduced (noise) model variables -----*/
00441   char * nname = NULL;         /* noise model name */
00442   vfp nmodel;                  /* pointer to noise model */
00443   int r;                       /* number of parameters in the noise model */
00444   char ** npname = NULL;       /* noise parameter labels */
00445   float * par_rdcd = NULL;     /* estimated parameters for the reduced model */
00446   float sse_rdcd;              /* error sum of squares for the reduced model */
00447   float * min_nconstr = NULL;  /* min parameter constraints for noise model */
00448   float * max_nconstr = NULL;  /* max parameter constraints for noise model */
00449 
00450   /*----- declare full (signal+noise) model variables -----*/
00451   char * sname = NULL;         /* signal model name */
00452   vfp smodel;                  /* pointer to signal model */
00453   int p;                       /* number of parameters in the signal model */
00454   char ** spname = NULL;       /* signal parameter labels */
00455   float * par_full = NULL;     /* estimated parameters for the full model */
00456   float sse_full;              /* error sum of squares for the full model */
00457   float * tpar_full = NULL;    /* t-statistic of parameters in full model */
00458   float freg;                  /* f-statistic for the full regression model */
00459   float rmsreg;                /* rms for the full regression model */
00460   float rsqr;                  /* R^2 (coef. of multiple determination) */
00461   float smax;                  /* signed maximum of signal */
00462   float tmax;                  /* epoch of signed maximum of signal */
00463   float pmax;                  /* percentage change due to signal */
00464   float area;                  /* area between signal and baseline */
00465   float parea;                 /* percent area between signal and baseline */
00466   float * min_sconstr = NULL;  /* min parameter constraints for signal model */
00467   float * max_sconstr = NULL;  /* max parameter constraints for signal model */
00468 
00469   int novar;               /* flag for insufficient variation in the data */
00470 
00471    
00472   /*----- program initialization -----*/
00473   initialize_program (&im1, &nname, &sname, &nmodel, &smodel, 
00474                       &r, &p, &npname, &spname,
00475                       &min_nconstr, &max_nconstr, &min_sconstr, &max_sconstr,
00476                       &nabs, &nrand, &nbest, &rms_min, 
00477                       &par_rdcd, &par_full, &tpar_full, 
00478                       ts_length, &tfilename, &x_array, &fit);
00479   
00480 
00481   /*----- calculate the reduced (noise) model -----*/  
00482   calc_reduced_model (ts_length, r, x_array, ts_array, 
00483                       par_rdcd, &sse_rdcd);
00484       
00485 
00486   /*----- calculate the full (signal+noise) model -----*/
00487   calc_full_model (nmodel, smodel, r, p, 
00488                    min_nconstr, max_nconstr, min_sconstr, max_sconstr,
00489                    ts_length, x_array, ts_array, par_rdcd, sse_rdcd, 
00490                    nabs, nrand, nbest, rms_min, par_full, &sse_full, &novar);
00491 
00492 
00493   /*----- create estimated time series using the full model parameters -----*/
00494   full_model (nmodel, smodel, par_full, par_full + r, 
00495               ts_length, x_array, fit);
00496       
00497 
00498   /*----- calculate statistics for the full model -----*/
00499   analyze_results (nmodel, smodel, r, p, novar,
00500                    min_nconstr, max_nconstr, min_sconstr, max_sconstr, 
00501                    ts_length, x_array,
00502                    par_rdcd, sse_rdcd, par_full, sse_full,
00503                    &rmsreg, &freg, &rsqr, &smax, &tmax, &pmax, &area, &parea,
00504                    tpar_full);
00505 
00506 
00507   /*----- report results for this voxel -----*/
00508   report_results (nname, sname, r, p, npname, spname, ts_length,
00509                   par_rdcd, sse_rdcd, par_full, sse_full, tpar_full,
00510                   rmsreg, freg, rsqr, smax, tmax, pmax, area, parea, label);
00511   printf ("\nVoxel Results: \n");
00512   printf ("%s \n", *label);
00513   
00514 
00515   /*----- end of program -----*/
00516   terminate_program (r, p, ts_length, &x_array,
00517                      &par_rdcd, &min_nconstr, &max_nconstr, 
00518                      &par_full, &tpar_full, 
00519                      &min_sconstr, &max_sconstr); 
00520 
00521   return (fit);
00522   
00523 }
00524 
00525 
00526 /***********************************************************************
00527    Set up the interface to the user:
00528     1) Create a new interface using "PLUTO_new_interface";
00529 
00530     2) For each line of inputs, create the line with "PLUTO_add_option"
00531          (this line of inputs can be optional or mandatory);
00532 
00533     3) For each item on the line, create the item with
00534         "PLUTO_add_dataset"    for a dataset chooser,
00535         "PLUTO_add_string"     for a string chooser,
00536         "PLUTO_add_number"     for a number chooser,
00537         "PLUTO_add_timeseries" for a timeseries chooser.
00538 ************************************************************************/
00539 
00540 
00541 DEFINE_PLUGIN_PROTOTYPE
00542 
00543 PLUGIN_interface * PLUGIN_init( int ncall )
00544 {
00545    int ii ;
00546    PLUGIN_interface * plint ;     /* will be the output of this routine */
00547    int ok;
00548 
00549    NLFIT_MODEL_array * model_array = NULL;   /* array of SO models */
00550    int im;                                   /* model index */
00551    int ip;                                   /* parameter index */
00552 
00553    char message[MAX_NAME_LENGTH];    /* error message */
00554 
00555 
00556    if( ncall > 0 ) return NULL ;  /* generate interfaces for ncall 0 */
00557 
00558    jump_on_NLfit_error = 1 ;                 /* 01 May 2003: */
00559    if( setjmp(NLfit_error_jmpbuf) != 0 ){    /* NLfit_error() was invoked */
00560      jump_on_NLfit_error = 0 ;               /* somewhere below here */
00561      fprintf(stderr,"\n*** Can't load NLfit plugin! ***\n");
00562      return NULL ;
00563    }
00564 
00565    /***** otherwise, do interface # 0 *****/
00566 
00567    /*---------------- set titles and call point ----------------*/
00568 
00569    plint = PLUTO_new_interface( "NLfit & NLerr" ,
00570                                 "Control NLfit and NLerr Functions" ,
00571                                 helpstring ,
00572                                 PLUGIN_CALL_VIA_MENU , NL_main ) ;
00573 
00574    PLUTO_add_hint( plint , "Control NLfit and NLerr Functions" ) ;
00575 
00576    global_plint = plint ;  /* make global copy */
00577 
00578    PLUTO_set_sequence( plint , "A:funcs:fitting" ) ;
00579 
00580    PLUTO_set_runlabels( plint , "Set+Keep" , "Set+Close" ) ;  /* 04 Nov 2003 */
00581 
00582    /*----- initialize the model array -----*/
00583    model_array = NLFIT_get_many_MODELs ();
00584    if ((model_array == NULL) || (model_array->num == 0))
00585 #if 1
00586      { PLUTO_report( plint , "Found no models!") ;
00587        jump_on_NLfit_error = 0 ; return NULL ; }
00588 #else
00589      NLfit_error ("Unable to locate any models");
00590 #endif
00591 
00592 #if 1
00593    else
00594    { char str[64] ;
00595      sprintf(str,"Found %d models",model_array->num) ;
00596      PLUTO_report( plint , str ) ;
00597    }
00598 #endif
00599 
00600    /*----- read parameters for noise models -----*/
00601    ii = 0;
00602    for (im = 0;  im < model_array->num;  im++)
00603      {
00604        if (model_array->modar[im]->interface->model_type == MODEL_NOISE_TYPE)
00605          {
00606            noise_labels[ii] = (char *) malloc (sizeof(char)*MAX_NAME_LENGTH);
00607            strncpy (noise_labels[ii], model_array->modar[im]->interface->label,
00608                     MAX_NAME_LENGTH);
00609 
00610            plug_nmodel[ii] = model_array->modar[im]->interface->call_func;
00611            if (plug_nmodel[ii] == NULL)
00612             {
00613               sprintf (message, "Noise model %s improperly defined. \n", 
00614                        noise_labels[ii]);
00615               NLfit_error (message);
00616             }
00617                   
00618            plug_r[ii] = model_array->modar[im]->interface->params;
00619            if ((plug_r[ii] < 0) || (plug_r[ii] > MAX_PARAMETERS))
00620             {
00621               sprintf (message, 
00622                        "Illegal number of parameters for noise model %s", 
00623                        noise_labels[ii]);
00624               NLfit_error (message);
00625             }
00626           
00627            for (ip = 0;  ip < plug_r[ii];  ip++)
00628              {
00629                noise_plabels[ii][ip] = 
00630                  (char *) malloc (sizeof(char)*MAX_NAME_LENGTH);
00631                strncpy (noise_plabels[ii][ip], 
00632                        model_array->modar[im]->interface->plabel[ip],
00633                         MAX_NAME_LENGTH);
00634                plug_min_nconstr[ii][ip] =
00635                  model_array->modar[im]->interface->min_constr[ip]; 
00636                plug_max_nconstr[ii][ip] = 
00637                  model_array->modar[im]->interface->max_constr[ip]; 
00638                if (plug_min_nconstr[ii][ip] > plug_max_nconstr[ii][ip])
00639                  NLfit_error
00640                    ("Must have noise parameter min cnstrnts <= max cnstrnts");
00641              }
00642            ii += 1;
00643          }
00644      }
00645    num_noise_models = ii;
00646    if (num_noise_models <= 0)  
00647      NLfit_error ("Unable to locate any noise models");
00648    plug_noise_index = 1;
00649 
00650 
00651    /*----- read parameters for signal models -----*/
00652    ii = 0;
00653    for (im = 0;  im < model_array->num;  im++)
00654      {
00655        if (model_array->modar[im]->interface->model_type == MODEL_SIGNAL_TYPE)
00656          {
00657            signal_labels[ii] = (char *) malloc (sizeof(char)*MAX_NAME_LENGTH);
00658            strncpy (signal_labels[ii],
00659                     model_array->modar[im]->interface->label,
00660                     MAX_NAME_LENGTH);
00661 
00662            plug_smodel[ii] = model_array->modar[im]->interface->call_func;
00663            if (plug_smodel[ii] == NULL)
00664             {
00665               sprintf (message, "Signal model %s improperly defined. \n", 
00666                        signal_labels[ii]);
00667               NLfit_error (message);
00668             }
00669                   
00670            plug_p[ii] = model_array->modar[im]->interface->params;
00671            if ((plug_p[ii] < 0) || (plug_p[ii] > MAX_PARAMETERS))
00672             {
00673               sprintf (message, 
00674                        "Illegal number of parameters for signal model %s", 
00675                        signal_labels[ii]);
00676               NLfit_error (message);
00677             }
00678           
00679 
00680            for (ip = 0;  ip < plug_p[ii];  ip++)
00681              {
00682                signal_plabels[ii][ip] = 
00683                  (char *) malloc (sizeof(char)*MAX_NAME_LENGTH);
00684                strncpy (signal_plabels[ii][ip], 
00685                         model_array->modar[im]->interface->plabel[ip],
00686                         MAX_NAME_LENGTH);
00687                plug_min_sconstr[ii][ip] =
00688                  model_array->modar[im]->interface->min_constr[ip]; 
00689                plug_max_sconstr[ii][ip] = 
00690                  model_array->modar[im]->interface->max_constr[ip]; 
00691                if (plug_min_sconstr[ii][ip] > plug_max_sconstr[ii][ip])
00692                  NLfit_error
00693                    ("Must have signal parameter min cnstrnts <= max cnstrnts");
00694              }
00695            ii += 1;
00696          }
00697      }
00698    num_signal_models = ii;
00699    if (num_signal_models <= 0)  
00700      NLfit_error ("Unable to locate any signal models");
00701    plug_signal_index = 0;
00702 
00703 
00704    /*----- Control -----*/
00705    PLUTO_add_option (plint , "Control" , "Control" , TRUE);
00706    PLUTO_add_number (plint , "Ignore" , 0,20,0,
00707                      plug_ignore , FALSE );
00708    PLUTO_add_number (plint , "NRandom" , 10,99999,0,
00709                      plug_nrand , TRUE );
00710    PLUTO_add_number (plint , "NBest" , 1,10,0,
00711                      plug_nbest , FALSE);
00712 
00713 
00714    /*----- Models -----*/
00715    PLUTO_add_option (plint , "Models" , "Models" , TRUE);
00716    PLUTO_add_string (plint, "Noise Model", num_noise_models, noise_labels, 
00717                      plug_noise_index);
00718    PLUTO_add_string (plint, "Signal Model", num_signal_models, signal_labels,
00719                      plug_signal_index);
00720    PLUTO_add_string (plint, "Noise Constr", 2, constr_types, 0);
00721 
00722 
00723    /*----- Noise Model Parameters -----*/
00724    PLUTO_add_option (plint, "Noise", "Noise", FALSE);
00725    PLUTO_add_number (plint, "Parameter" , 0, MAX_PARAMETERS, 0, 0, FALSE);
00726    PLUTO_add_number (plint, "Min Constr", -99999, 99999, 0, 0, TRUE);
00727    PLUTO_add_number (plint, "Max Constr", -99999, 99999, 0, 0, TRUE);
00728 
00729 
00730    /*----- Signal Model Parameters -----*/
00731    PLUTO_add_option (plint, "Signal", "Signal", FALSE);
00732    PLUTO_add_number (plint, "Parameter", 0, MAX_PARAMETERS, 0, 0, FALSE);
00733    PLUTO_add_number (plint, "Min Constr", -99999, 99999, 0, 0, TRUE);
00734    PLUTO_add_number (plint, "Max Constr", -99999, 99999, 0, 0, TRUE);
00735 
00736 
00737    /*----- External Time Reference -----*/
00738    PLUTO_add_option (plint, "Time Scale", "Time Scale", FALSE);
00739    PLUTO_add_string (plint, "Reference", 3, time_refs, 0);
00740    PLUTO_add_string (plint, "File", 0, NULL, 19);
00741 
00742 
00743    /*--------- done with interface setup ---------*/
00744 
00745    PLUTO_register_1D_funcstr ("NLfit" , NL_fitter);
00746    PLUTO_register_1D_funcstr ("NLerr" , NL_error);
00747 
00748 
00749    /*----- discard the model array -----*/
00750 #if 0
00751    DESTROY_MODEL_ARRAY (model_array);
00752 #endif
00753 
00754    jump_on_NLfit_error = 0 ; return plint ;
00755 }
00756 
00757 /***************************************************************************
00758   Main routine for this plugin (will be called from AFNI).
00759   If the return string is not NULL, some error transpired, and
00760   AFNI will popup the return string in a message box.
00761 ****************************************************************************/
00762 
00763 char * NL_main( PLUGIN_interface * plint )
00764 {
00765    char * str ;
00766    int  ii, ival, ip ;
00767    float * tsar ;
00768    float min_constr, max_constr;
00769    MRI_IMAGE * im;  /* pointer to image structures */
00770 
00771 
00772    /*--------- go to first input line ---------*/
00773 
00774    PLUTO_next_option(plint) ;
00775 
00776    plug_ignore = PLUTO_get_number(plint) ;
00777    plug_nrand = PLUTO_get_number(plint) ;
00778    plug_nbest = PLUTO_get_number(plint) ;
00779 
00780 
00781    /*------ loop over remaining options, check their tags, process them -----*/
00782 
00783    do 
00784      {
00785        str = PLUTO_get_optiontag(plint) ; 
00786        if( str == NULL ) break ;
00787 
00788        if( strcmp(str,"Models") == 0 )
00789          {
00790            str = PLUTO_get_string(plint) ;
00791            for (ii = 0;  ii < num_noise_models;  ii++)
00792              if (strcmp (str, noise_labels[ii]) == 0)
00793                plug_noise_index = ii;
00794          
00795            str = PLUTO_get_string(plint) ;
00796            for (ii = 0;  ii < num_signal_models;  ii++)
00797              if (strcmp (str, signal_labels[ii]) == 0)
00798                plug_signal_index = ii;
00799 
00800            str = PLUTO_get_string(plint);
00801            if (strcmp (str, "Absolute") == 0)
00802              plug_nabs = 1;
00803            else
00804              plug_nabs = 0;
00805          } 
00806 
00807        else if( strcmp(str,"Noise") == 0 )
00808          {
00809            ival = PLUTO_get_number(plint);
00810            min_constr = PLUTO_get_number(plint);
00811            max_constr = PLUTO_get_number(plint);
00812            if (min_constr > max_constr)
00813              return "**********************************\n"
00814                     " Require min constr <= max constr \n"
00815                     "**********************************"  ;
00816            plug_min_nconstr[plug_noise_index][ival] = min_constr;
00817            plug_max_nconstr[plug_noise_index][ival] = max_constr;
00818          } 
00819 
00820        else if( strcmp(str,"Signal") == 0 )
00821          {
00822            ival = PLUTO_get_number(plint);
00823            min_constr = PLUTO_get_number(plint);
00824            max_constr = PLUTO_get_number(plint);
00825            if (min_constr > max_constr)
00826              return "**********************************\n"
00827                     " Require min constr <= max constr \n"
00828                     "**********************************"  ;
00829            plug_min_sconstr[plug_signal_index][ival] = min_constr;
00830            plug_max_sconstr[plug_signal_index][ival] = max_constr;
00831          } 
00832 
00833        else if( strcmp(str,"Time Scale") == 0 )
00834          {
00835            str = PLUTO_get_string(plint);
00836            if (strcmp (str, "External") == 0){
00837                plug_timeref = 1;
00838                str = PLUTO_get_string(plint);
00839                im = mri_read_1D (str); 
00840                if (im == NULL)
00841                  return "************************************\n"
00842                         " Unable to read time reference file \n"
00843                         "************************************"  ;
00844                mri_free(im);
00845                strcpy (plug_tfilename, str);
00846 
00847            } else if( strcmp(str,"-inTR") == 0 ){  /* 22 July 1998 */
00848               inTR = 1 ;
00849               plug_timeref = 0;
00850 
00851            } else {
00852              plug_timeref = 0;
00853              inTR = 0 ;                       /* 22 July 1998 */
00854            }
00855 
00856          } 
00857 
00858        else 
00859          {
00860            return "************************\n"
00861                   "Illegal optiontag found!\n"
00862                   "************************"  ;
00863          }
00864      } while(1) ;
00865 
00866 
00867   
00868   /*----- Identify software -----*/
00869   printf ("\n\n");
00870   printf ("Program: %s \n", PROGRAM_NAME);
00871   printf ("Author:  %s \n", PROGRAM_AUTHOR); 
00872   printf ("Date:    %s \n", PROGRAM_DATE);
00873   printf ("\n");
00874 
00875    
00876    /*----- show current input options -----*/
00877    printf ("\nControls: \n");
00878    printf ("Ignore       = %5d \n", plug_ignore);
00879    printf ("Num Random   = %5d \n", plug_nrand);
00880    printf ("Num Best     = %5d \n", plug_nbest);
00881    printf ("Noise Constr = %s  \n", constr_types[plug_nabs]);
00882    printf ("\nNoise  Model = %s \n", noise_labels[plug_noise_index]);
00883    for (ip = 0;  ip < plug_r[plug_noise_index];  ip++)
00884      {
00885        printf ("gn[%d]:   min =%10.3f   max =%10.3f   %s \n", 
00886                ip, plug_min_nconstr[plug_noise_index][ip], 
00887                plug_max_nconstr[plug_noise_index][ip],
00888                noise_plabels[plug_noise_index][ip]);
00889      }
00890    printf ("\nSignal Model = %s \n", signal_labels[plug_signal_index]);
00891    for (ip = 0;  ip < plug_p[plug_signal_index];  ip++)
00892      {
00893        printf ("gs[%d]:   min =%10.3f   max =%10.3f   %s \n", 
00894                ip, plug_min_sconstr[plug_signal_index][ip], 
00895                plug_max_sconstr[plug_signal_index][ip],
00896                signal_plabels[plug_signal_index][ip]);
00897      }
00898 
00899    if (plug_timeref)
00900      printf ("\nExternal Time Reference = %s \n", plug_tfilename);
00901    else if( inTR )
00902      printf ("\n-inTR Time Reference\n") ;
00903    else
00904      printf ("\nInternal Time Reference \n");
00905    
00906    
00907    /*--- nothing left to do until data arrives ---*/
00908    
00909    initialize = 1 ;  /* force re-initialization */
00910    
00911    
00912    return NULL ;
00913 }
00914 
00915 
00916 /*---------------------------------------------------------------------------*/
00917 
00918 void NL_fitter( int nt , double to , double dt , float * vec, char ** label )
00919 {
00920    NL_worker( nt , dt , vec , TRUE, label ) ;
00921    return ;
00922 }
00923 
00924 
00925 /*---------------------------------------------------------------------------*/
00926 
00927 void NL_error( int nt , double to , double dt , float * vec, char ** label )
00928 {
00929    NL_worker( nt , dt , vec , FALSE, label ) ;
00930    return ;
00931 }
00932  
00933 
00934 /*---------------------------------------------------------------------------*/
00935 
00936 void NL_worker( int nt , double dt , float * vec , int dofit, char ** label )
00937 {
00938    float * fit;
00939    int ii, nlen;
00940    float val;
00941 
00942 
00943    nlen = nt - plug_ignore;
00944 
00945    dsTR = dt ;
00946 
00947    /** find least squares fit coefficients **/
00948 
00949    fit = nlfit (nlen, vec+plug_ignore, label);
00950 
00951    for (ii = 0;  ii < plug_ignore;  ii++)
00952      if (dofit)
00953        vec[ii] = fit[0];
00954      else
00955        vec[ii] = vec[plug_ignore] - fit[0];
00956 
00957    for (ii=plug_ignore; ii < nt; ii++)
00958      {
00959        if (dofit)
00960          vec[ii] = fit[ii-plug_ignore];
00961        else
00962          vec[ii] = vec[ii] - fit[ii-plug_ignore] ;
00963      }
00964 
00965    free(fit) ;
00966    return ;
00967 }
00968 
00969 
00970 
00971 
00972 
00973 
00974 
00975 
00976 
00977 
00978