/***************************************************************************** Major portions of this software are copyrighted by the Medical College of Wisconsin, 1994-2000, and are released under the Gnu General Public License, Version 2. See the file README.Copyright for details. ******************************************************************************/ /* Plugin to calculate a nonlinear regression for an input time series, using a specified nonlinear model. File: plug_nlfit.c Author: B. Douglas Ward Date: 17 June 1997 Mod: Added external time reference capability and added option for absolute noise parameter constraints. Date: 22 August 1997 Mod: Added options for percent signal change above baseline, and calculation of area under the signal above baseline. Date: 26 November 1997 Mod: Print error message if unable to locate any signal models or any noise models. Date: 26 December 1997 Mod: Added novar flag to eliminate unnecessary calculations. Date: 13 July 1999 Mod: Adjust F-statistics if parameter constraints force a parameter to be a constant. Date: 08 February 2000 Mod: Changes for output of R^2 (coefficient of multiple determination), and standard deviation of residuals from full model fit. Added global variable calc_tstats. Also, added screen display of p-values. Date: 10 May 2000 */ /*---------------------------------------------------------------------------*/ #define PROGRAM_NAME "plug_nlfit" /* name of this program */ #define PROGRAM_AUTHOR "B. Douglas Ward" /* program author */ #define PROGRAM_DATE "10 May 2000" /* date of last program mod */ /*---------------------------------------------------------------------------*/ #include "afni.h" #ifndef ALLOW_PLUGINS # error "Plugins not properly set up -- see machdep.h" #endif #include #include #include "mrilib.h" #include "matrix.h" #include "simplex.h" #include "NLfit.h" #include "matrix.c" #include "simplex.c" #include "NLfit.c" /***** 22 July 1998 -- RWCox: Modified to allow DELT to be set from the TR of the input file *****/ static float DELT = 1.0; /* default */ static int inTR = 0 ; /* set to 1 if -inTR option is used */ static float dsTR = 0.0 ; /* TR of the input file */ /*********************************************************************** Plugin to provide nonlinear least squares fitting 1D function for graphs ************************************************************************/ /*------------- string to 'help' the user -------------*/ static char helpstring[] = " Purpose: Control the 'NLfit' and 'NLerr' 1D functions.\n" "\n" " Control: Ignore = Number of points to ignore at start \n" " of each timeseries. \n" " NRandom = Number of random test points. \n" " NBest = Find opt. soln. for these best test points. \n" " \n" " Models: Noise = Label of noise model to use. \n" " Signal = Label of signal model to use. \n" " Noise Constr = Relative or Absolute noise constraints. \n" " \n" " Noise: Parameter = Index of noise model parameter. \n" " Min Constr = Minimum value for noise model parameter. \n" " Max Constr = Maximum value for noise model parameter. \n" " \n" " Signal: Parameter = Index of signal model parameter.\n" " Min Constr = Minimum value for signal model parameter. \n" " Max Constr = Maximum value for signal model parameter. \n" " \n" " Time Scale: Reference = Internal or External time reference. \n" " File = External time reference file name. \n" "Author -- BD Ward" ; /*--------------- prototypes for internal routines ---------------*/ char * NL_main( PLUGIN_interface * ) ; /* the entry point */ void NL_fitter( int nt, double to, double dt, float * vec, char ** label ) ; void NL_error ( int nt, double to, double dt, float * vec, char ** label ) ; void NL_worker( int nt, double dt, float * vec, int dofit, char ** label ) ; /*---------------- global data -------------------*/ static PLUGIN_interface * global_plint = NULL ; static int initialize=1 ; /*------- Global data for models ------*/ char * constr_types[2] = {"Relative", "Absolute"}; /* option labels */ char * time_refs[3] = {"Internal", "External" , "-inTR" }; int plug_ignore = 3; int plug_nrand = 100; int plug_nbest = 5; int plug_nabs = 0; int plug_timeref = 0; char plug_tfilename[MAX_NAME_LENGTH] = ""; /*----- declare reduced (noise) model variables -----*/ int num_noise_models; /* number of noise models */ int plug_noise_index; /* index of noise model */ char * noise_labels[MAX_MODELS]; /* names of noise models */ vfp plug_nmodel[MAX_MODELS]; /* pointer to noise model */ int plug_r[MAX_MODELS]; /* number of parameters in the noise model */ char * noise_plabels[MAX_MODELS][MAX_PARAMETERS]; float plug_min_nconstr[MAX_MODELS][MAX_PARAMETERS]; /* minimum parameter constraints for noise model */ float plug_max_nconstr[MAX_MODELS][MAX_PARAMETERS]; /* maximum parameter constraints for noise model */ /*----- declare full (signal+noise) model variables -----*/ int num_signal_models; /* number of signal models */ int plug_signal_index; /* index of signal model */ char * signal_labels[MAX_MODELS]; /* names of signal models */ vfp plug_smodel[MAX_MODELS]; /* pointer to signal model */ int plug_p[MAX_MODELS]; /* number of parameters in the signal model */ char * signal_plabels[MAX_MODELS][MAX_PARAMETERS]; float plug_min_sconstr[MAX_MODELS][MAX_PARAMETERS]; /* minimum parameter constraints for signal model */ float plug_max_sconstr[MAX_MODELS][MAX_PARAMETERS]; /* maximum parameter constraints for signal model */ /*---------------------------------------------------------------------------*/ /* Routine to initialize the input options. */ void initialize_options ( int * im1, /* index of 1st image in time series for analysis */ char ** nname, /* noise model name */ char ** sname, /* signal model name */ vfp * nmodel, /* pointer to noise model */ vfp * smodel, /* pointer to signal model */ int * r, /* number of parameters in the noise model */ int * p, /* number of parameters in the signal model */ char *** npname, /* noise parameter names */ char *** spname, /* signal parameter names */ float ** min_nconstr, /* minimum parameter constraints for noise model */ float ** max_nconstr, /* maximum parameter constraints for noise model */ float ** min_sconstr, /* minimum parameter constraints for signal model */ float ** max_sconstr, /* maximum parameter constraints for signal model */ int * nabs, /* use absolute constraints for noise parameters */ int * nrand, /* number of random vectors to generate */ int * nbest, /* number of random vectors to keep */ float * rms_min, /* minimum rms error to reject reduced model */ char ** tfilename /* file name for time point series */ ) { int ip; /* parameter index */ int ok; /* boolean for specified model exists */ char message[MAX_NAME_LENGTH]; /* error message */ *im1 = 1; *nrand = plug_nrand; *nbest = plug_nbest; *nabs = plug_nabs; *rms_min = 0.0; *tfilename = plug_tfilename; *nname = noise_labels[plug_noise_index]; *sname = signal_labels[plug_signal_index]; *nmodel = plug_nmodel[plug_noise_index]; *smodel = plug_smodel[plug_signal_index]; *r = plug_r[plug_noise_index]; *p = plug_p[plug_signal_index]; *npname = noise_plabels[plug_noise_index]; *spname = signal_plabels[plug_signal_index]; /*----- allocate memory for parameter constraints -----*/ *min_nconstr = (float *) malloc (sizeof(float) * (*r)); if (*min_nconstr == NULL) NLfit_error ("Unable to allocate memory for min_nconstr"); *max_nconstr = (float *) malloc (sizeof(float) * (*r)); if (*max_nconstr == NULL) NLfit_error ("Unable to allocate memory for max_nconstr"); *min_sconstr = (float *) malloc (sizeof(float) * (*p)); if (*min_sconstr == NULL) NLfit_error ("Unable to allocate memory for min_sconstr"); *max_sconstr = (float *) malloc (sizeof(float) * (*p)); if (*max_sconstr == NULL) NLfit_error ("Unable to allocate memory for max_sconstr"); /*----- initialize constraints -----*/ for (ip = 0; ip < (*r); ip++) { (*min_nconstr)[ip] = plug_min_nconstr[plug_noise_index][ip]; (*max_nconstr)[ip] = plug_max_nconstr[plug_noise_index][ip]; } for (ip = 0; ip < (*p); ip++) { (*min_sconstr)[ip] = plug_min_sconstr[plug_signal_index][ip]; (*max_sconstr)[ip] = plug_max_sconstr[plug_signal_index][ip]; } } /*---------------------------------------------------------------------------*/ /* Routine to check for valid inputs. */ void check_for_valid_inputs () { } /*---------------------------------------------------------------------------*/ void initialize_program ( int * im1, /* index of 1st image in time series for analysis */ char ** nname, /* noise model name */ char ** sname, /* signal model name */ vfp * nmodel, /* pointer to noise model */ vfp * smodel, /* pointer to signal model */ int * r, /* number of parameters in the noise model */ int * p, /* number of parameters in the signal model */ char *** npname, /* noise parameter names */ char *** spname, /* signal parameter names */ float ** min_nconstr, /* minimum parameter constraints for noise model */ float ** max_nconstr, /* maximum parameter constraints for noise model */ float ** min_sconstr, /* minimum parameter constraints for signal model */ float ** max_sconstr, /* maximum parameter constraints for signal model */ int * nabs, /* use absolute constraints for noise parameters */ int * nrand, /* number of random vectors to generate */ int * nbest, /* number of random vectors to keep */ float * rms_min, /* minimum rms error to reject reduced model */ float ** par_rdcd, /* estimated parameters for the reduced model */ float ** par_full, /* estimated parameters for the full model */ float ** tpar_full, /* t-statistic of parameters in the full model */ int ts_length, /* length of time series data */ char ** tfilename, /* file name for time point series */ float *** x_array, /* independent variable matrix */ float ** fit ) { int dimension; /* dimension of full model */ int ip; /* parameter index */ int it; /* time index */ MRI_IMAGE * im, * flim; /* pointers to image structures -- used to read 1D ASCII */ int nt; /* number of points in 1D x data file */ float * tar; /*----- intialize options -----*/ initialize_options (im1, nname, sname, nmodel, smodel, r, p, npname, spname, min_nconstr, max_nconstr, min_sconstr, max_sconstr, nabs, nrand, nbest, rms_min, tfilename); /*----- check for valid inputs -----*/ check_for_valid_inputs (); /*----- allocate space for independent variable matrix -----*/ *x_array = (float **) malloc (sizeof(float *) * ts_length); if (*x_array == NULL) NLfit_error ("Unable to allocate memory for x_array"); for (it = 0; it < ts_length; it++) { (*x_array)[it] = (float *) malloc (sizeof(float) * 3); if ((*x_array)[it] == NULL) NLfit_error ("Unable to allocate memory for x_array[it]"); } /*----- initialize independent variable matrix -----*/ if (!plug_timeref) { static float old_DELT = -1.0 ; DELT = (inTR && dsTR > 0.0) ? dsTR : 1.0 ; /* 22 July 1998 */ if( DELT != old_DELT ){ old_DELT = DELT ; printf("NLfit: switch to TR = %g\n",DELT) ; } for (it = 0; it < ts_length; it++) { (*x_array)[it][0] = 1.0; (*x_array)[it][1] = it * DELT; (*x_array)[it][2] = (it * DELT) * (it * DELT); } } else { flim = mri_read_1D (*tfilename); if (flim == NULL) NLfit_error ("Unable to read time reference file \n"); nt = flim -> nx; if (nt < ts_length) NLfit_error ("Time reference array is too short"); tar = MRI_FLOAT_PTR(flim) ; for (it = 0; it < ts_length; it++) { (*x_array)[it][0] = 1.0; (*x_array)[it][1] = tar[it] ; (*x_array)[it][2] = tar[it] * tar[it]; } mri_free (flim); } /*--- 24 Jul 2006: special change to x_array[][2] for Linear+Ort [RWCox] ---*/ if( strcmp(*nname,"Linear+Ort") == 0 ){ char *fname=NULL; MRI_IMAGE *fim=NULL; int nx; float *far; static int nwarn=0; fname = my_getenv("AFNI_ORTMODEL_REF") ; if( fname == NULL ){ ERROR_message("Linear+Ort model: 'AFNI_ORTMODEL_REF' not set") ; goto PLO_done ; } fim = mri_read_1D(fname) ; if( fim == NULL || fim->nx < 2 ){ ERROR_message( "Linear+Ort model: can't read AFNI_ORTMODEL_REF='%s'",fname) ; goto PLO_done ; } if( fim->ny > 1 && nwarn < 2 ){ WARNING_message( "Linear+Ort model: file AFNI_ORTMODEL_REF='%s' has more than 1 column", fname ) ; nwarn++ ; } nx = fim->nx ; far = MRI_FLOAT_PTR(fim) ; if( nx != ts_length && nwarn ){ WARNING_message("Linear+Ort: length(%s)=%d but length(dataset)=%d", fname , nx , ts_length ) ; nwarn++ ; } for( it=0 ; it < ts_length; it++) (*x_array)[it][2] = (it < nx) ? far[it] : 0.0f ; PLO_done: ; /* nada */ } dimension = (*r) + (*p); /*----- allocate memory space -----*/ *par_rdcd = (float *) malloc (sizeof(float) * dimension); if (*par_rdcd == NULL) NLfit_error ("Unable to allocate memory for par_rdcd"); *par_full = (float *) malloc (sizeof(float) * dimension); if (*par_full == NULL) NLfit_error ("Unable to allocate memory for par_full"); *tpar_full = (float *) malloc (sizeof(float) * dimension); if (*tpar_full == NULL) NLfit_error ("Unable to allocate memory for tpar_full"); *fit = (float *) malloc (sizeof(float) * (ts_length)); if (*fit == NULL) NLfit_error ("Unable to allocate memory for fit"); } /*---------------------------------------------------------------------------*/ /* Release all allocated memory space. */ void terminate_program ( int r, /* number of parameters in the noise model */ int p, /* number of parameters in the signal model */ int ts_length, /* length of time series data */ float *** x_array, /* independent variable matrix */ float ** par_rdcd, /* estimated parameters for the reduced model */ float ** min_nconstr, /* min parameter constraints for noise model */ float ** max_nconstr, /* max parameter constraints for noise model */ float ** par_full, /* estimated parameters for the full model */ float ** tpar_full, /* t-statistic of parameters in full model */ float ** min_sconstr, /* min parameter constraints for signal model */ float ** max_sconstr /* max parameter constraints for signal model */ ) { int ip; /* parameter index */ int it; /* time index */ /*----- deallocate memory for parameter constraints -----*/ if (*min_nconstr != NULL) { free (*min_nconstr); *min_nconstr = NULL; } if (*max_nconstr != NULL) { free (*max_nconstr); *max_nconstr = NULL; } if (*min_sconstr != NULL) { free (*min_sconstr); *min_sconstr = NULL; } if (*max_sconstr != NULL) { free (*max_sconstr); *max_sconstr = NULL; } /*----- deallocate space for independent variable matrix -----*/ for (it = 0; it < ts_length; it++) if ((*x_array)[it] != NULL) { free ((*x_array)[it]); (*x_array)[it] = NULL; } if (*x_array != NULL) { free (*x_array); *x_array = NULL; } /*----- deallocate space for parameters -----*/ if (*par_rdcd != NULL) { free (*par_rdcd); *par_rdcd = NULL; } if (*par_full != NULL) { free (*par_full); *par_full = NULL; } if (*tpar_full != NULL) { free (*tpar_full); *tpar_full = NULL; } } /*---------------------------------------------------------------------------*/ float * nlfit ( int ts_length, /* length of time series data */ float * ts_array, /* input time series array */ char ** label /* label output for this voxel */ ) { float * fit; /* nonlinear fit of time series data */ /*----- declare input option variables -----*/ int nabs; /* use absolute constraints for noise parameters */ int nrand; /* number of random vectors to generate */ int nbest; /* number of random vectors to keep */ float rms_min; /* minimum rms error to reject reduced model */ /*----- declare time series variables -----*/ int im1; /* index of 1st image in time series for analysis */ float ** x_array = NULL; /* independent variable matrix */ char * tfilename = NULL; /* file name of time points */ /*----- declare reduced (noise) model variables -----*/ char * nname = NULL; /* noise model name */ vfp nmodel; /* pointer to noise model */ int r; /* number of parameters in the noise model */ char ** npname = NULL; /* noise parameter labels */ float * par_rdcd = NULL; /* estimated parameters for the reduced model */ float sse_rdcd; /* error sum of squares for the reduced model */ float * min_nconstr = NULL; /* min parameter constraints for noise model */ float * max_nconstr = NULL; /* max parameter constraints for noise model */ /*----- declare full (signal+noise) model variables -----*/ char * sname = NULL; /* signal model name */ vfp smodel; /* pointer to signal model */ int p; /* number of parameters in the signal model */ char ** spname = NULL; /* signal parameter labels */ float * par_full = NULL; /* estimated parameters for the full model */ float sse_full; /* error sum of squares for the full model */ float * tpar_full = NULL; /* t-statistic of parameters in full model */ float freg; /* f-statistic for the full regression model */ float rmsreg; /* rms for the full regression model */ float rsqr; /* R^2 (coef. of multiple determination) */ float smax; /* signed maximum of signal */ float tmax; /* epoch of signed maximum of signal */ float pmax; /* percentage change due to signal */ float area; /* area between signal and baseline */ float parea; /* percent area between signal and baseline */ float * min_sconstr = NULL; /* min parameter constraints for signal model */ float * max_sconstr = NULL; /* max parameter constraints for signal model */ int novar; /* flag for insufficient variation in the data */ /*----- program initialization -----*/ initialize_program (&im1, &nname, &sname, &nmodel, &smodel, &r, &p, &npname, &spname, &min_nconstr, &max_nconstr, &min_sconstr, &max_sconstr, &nabs, &nrand, &nbest, &rms_min, &par_rdcd, &par_full, &tpar_full, ts_length, &tfilename, &x_array, &fit); /*----- calculate the reduced (noise) model -----*/ calc_reduced_model (ts_length, r, x_array, ts_array, par_rdcd, &sse_rdcd); /*----- calculate the full (signal+noise) model -----*/ calc_full_model (nmodel, smodel, r, p, min_nconstr, max_nconstr, min_sconstr, max_sconstr, ts_length, x_array, ts_array, par_rdcd, sse_rdcd, nabs, nrand, nbest, rms_min, par_full, &sse_full, &novar); /*----- create estimated time series using the full model parameters -----*/ full_model (nmodel, smodel, par_full, par_full + r, ts_length, x_array, fit); /*----- calculate statistics for the full model -----*/ analyze_results (nmodel, smodel, r, p, novar, min_nconstr, max_nconstr, min_sconstr, max_sconstr, ts_length, x_array, par_rdcd, sse_rdcd, par_full, sse_full, &rmsreg, &freg, &rsqr, &smax, &tmax, &pmax, &area, &parea, tpar_full); /*----- report results for this voxel -----*/ report_results (nname, sname, r, p, npname, spname, ts_length, par_rdcd, sse_rdcd, par_full, sse_full, tpar_full, rmsreg, freg, rsqr, smax, tmax, pmax, area, parea, label); printf ("\nVoxel Results: \n"); printf ("%s \n", *label); /*----- end of program -----*/ terminate_program (r, p, ts_length, &x_array, &par_rdcd, &min_nconstr, &max_nconstr, &par_full, &tpar_full, &min_sconstr, &max_sconstr); return (fit); } /*********************************************************************** Set up the interface to the user: 1) Create a new interface using "PLUTO_new_interface"; 2) For each line of inputs, create the line with "PLUTO_add_option" (this line of inputs can be optional or mandatory); 3) For each item on the line, create the item with "PLUTO_add_dataset" for a dataset chooser, "PLUTO_add_string" for a string chooser, "PLUTO_add_number" for a number chooser, "PLUTO_add_timeseries" for a timeseries chooser. ************************************************************************/ DEFINE_PLUGIN_PROTOTYPE PLUGIN_interface * PLUGIN_init( int ncall ) { int ii ; PLUGIN_interface * plint ; /* will be the output of this routine */ int ok; NLFIT_MODEL_array * model_array = NULL; /* array of SO models */ int im; /* model index */ int ip; /* parameter index */ char message[MAX_NAME_LENGTH]; /* error message */ if( ncall > 0 ) return NULL ; /* generate interfaces for ncall 0 */ jump_on_NLfit_error = 1 ; /* 01 May 2003: */ if( setjmp(NLfit_error_jmpbuf) != 0 ){ /* NLfit_error() was invoked */ jump_on_NLfit_error = 0 ; /* somewhere below here */ fprintf(stderr,"\n*** Can't load NLfit plugin! ***\n"); return NULL ; } /***** otherwise, do interface # 0 *****/ /*---------------- set titles and call point ----------------*/ plint = PLUTO_new_interface( "NLfit & NLerr" , "Control NLfit and NLerr Functions" , helpstring , PLUGIN_CALL_VIA_MENU , NL_main ) ; { char *eee = getenv("AFNI_NLFIM_METHOD") , str[94] ; if( eee == NULL || strcasecmp(eee,"simplex") == 0 ) N_newuoa = 0 ; else if( strcasecmp(eee,"powell") == 0 ) N_newuoa = 1 ; else if( strcasecmp(eee,"both") == 0 ) N_newuoa = 2 ; else N_newuoa = 0 ; sprintf(str,"Optimizer (AFNI_NLFIM_METHOD) is %s" , (N_newuoa==0) ? "SIMPLEX" :(N_newuoa==1) ? "POWELL" : "BOTH (SIMPLEX+POWELL)" ) ; PLUTO_report(plint,str) ; } PLUTO_add_hint( plint , "Control NLfit and NLerr Functions" ) ; global_plint = plint ; /* make global copy */ PLUTO_set_sequence( plint , "A:funcs:fitting" ) ; PLUTO_set_runlabels( plint , "Set+Keep" , "Set+Close" ) ; /* 04 Nov 2003 */ /*----- initialize the model array -----*/ model_array = NLFIT_get_many_MODELs (); if ((model_array == NULL) || (model_array->num == 0)) #if 1 { PLUTO_report( plint , "Found no models!") ; jump_on_NLfit_error = 0 ; return NULL ; } #else NLfit_error ("Unable to locate any models"); #endif #if 1 else { char str[64] ; sprintf(str,"Found %d models",model_array->num) ; PLUTO_report( plint , str ) ; } #endif /*----- read parameters for noise models -----*/ ii = 0; for (im = 0; im < model_array->num; im++) { if (model_array->modar[im]->interface->model_type == MODEL_NOISE_TYPE) { noise_labels[ii] = (char *) malloc (sizeof(char)*MAX_NAME_LENGTH); strncpy (noise_labels[ii], model_array->modar[im]->interface->label, MAX_NAME_LENGTH); plug_nmodel[ii] = model_array->modar[im]->interface->call_func; if (plug_nmodel[ii] == NULL) { sprintf (message, "Noise model %s improperly defined. \n", noise_labels[ii]); NLfit_error (message); } plug_r[ii] = model_array->modar[im]->interface->params; if ((plug_r[ii] < 0) || (plug_r[ii] > MAX_PARAMETERS)) { sprintf (message, "Illegal number of parameters for noise model %s", noise_labels[ii]); NLfit_error (message); } for (ip = 0; ip < plug_r[ii]; ip++) { noise_plabels[ii][ip] = (char *) malloc (sizeof(char)*MAX_NAME_LENGTH); strncpy (noise_plabels[ii][ip], model_array->modar[im]->interface->plabel[ip], MAX_NAME_LENGTH); plug_min_nconstr[ii][ip] = model_array->modar[im]->interface->min_constr[ip]; plug_max_nconstr[ii][ip] = model_array->modar[im]->interface->max_constr[ip]; if (plug_min_nconstr[ii][ip] > plug_max_nconstr[ii][ip]) NLfit_error ("Must have noise parameter min cnstrnts <= max cnstrnts"); } ii += 1; } } num_noise_models = ii; if (num_noise_models <= 0) NLfit_error ("Unable to locate any noise models"); plug_noise_index = 1; /*----- read parameters for signal models -----*/ ii = 0; for (im = 0; im < model_array->num; im++) { if (model_array->modar[im]->interface->model_type == MODEL_SIGNAL_TYPE) { signal_labels[ii] = (char *) malloc (sizeof(char)*MAX_NAME_LENGTH); strncpy (signal_labels[ii], model_array->modar[im]->interface->label, MAX_NAME_LENGTH); plug_smodel[ii] = model_array->modar[im]->interface->call_func; if (plug_smodel[ii] == NULL) { sprintf (message, "Signal model %s improperly defined. \n", signal_labels[ii]); NLfit_error (message); } plug_p[ii] = model_array->modar[im]->interface->params; if ((plug_p[ii] < 0) || (plug_p[ii] > MAX_PARAMETERS)) { sprintf (message, "Illegal number of parameters for signal model %s", signal_labels[ii]); NLfit_error (message); } for (ip = 0; ip < plug_p[ii]; ip++) { signal_plabels[ii][ip] = (char *) malloc (sizeof(char)*MAX_NAME_LENGTH); strncpy (signal_plabels[ii][ip], model_array->modar[im]->interface->plabel[ip], MAX_NAME_LENGTH); plug_min_sconstr[ii][ip] = model_array->modar[im]->interface->min_constr[ip]; plug_max_sconstr[ii][ip] = model_array->modar[im]->interface->max_constr[ip]; if (plug_min_sconstr[ii][ip] > plug_max_sconstr[ii][ip]) NLfit_error ("Must have signal parameter min cnstrnts <= max cnstrnts"); } ii += 1; } } num_signal_models = ii; if (num_signal_models <= 0) NLfit_error ("Unable to locate any signal models"); plug_signal_index = 0; /*----- Control -----*/ PLUTO_add_option (plint , "Control" , "Control" , TRUE); PLUTO_add_number (plint , "Ignore" , 0,20,0, plug_ignore , FALSE ); PLUTO_add_number (plint , "NRandom" , 10,99999,0, plug_nrand , TRUE ); PLUTO_add_number (plint , "NBest" , 1,10,0, plug_nbest , FALSE); /*----- Models -----*/ PLUTO_add_option (plint , "Models" , "Models" , TRUE); PLUTO_add_string (plint, "Noise Model", num_noise_models, noise_labels, plug_noise_index); PLUTO_add_string (plint, "Signal Model", num_signal_models, signal_labels, plug_signal_index); PLUTO_add_string (plint, "Noise Constr", 2, constr_types, 0); /*----- Noise Model Parameters -----*/ PLUTO_add_option (plint, "Noise", "Noise", FALSE); PLUTO_add_number (plint, "Parameter" , 0, MAX_PARAMETERS, 0, 0, FALSE); PLUTO_add_number (plint, "Min Constr", -99999, 99999, 0, 0, TRUE); PLUTO_add_number (plint, "Max Constr", -99999, 99999, 0, 0, TRUE); /*----- Signal Model Parameters -----*/ PLUTO_add_option (plint, "Signal", "Signal", FALSE); PLUTO_add_number (plint, "Parameter", 0, MAX_PARAMETERS, 0, 0, FALSE); PLUTO_add_number (plint, "Min Constr", -99999, 99999, 0, 0, TRUE); PLUTO_add_number (plint, "Max Constr", -99999, 99999, 0, 0, TRUE); /*----- External Time Reference -----*/ PLUTO_add_option (plint, "Time Scale", "Time Scale", FALSE); PLUTO_add_string (plint, "Reference", 3, time_refs, 0); PLUTO_add_string (plint, "File", 0, NULL, 19); /*--------- done with interface setup ---------*/ PLUTO_register_1D_funcstr ("NLfit" , NL_fitter); PLUTO_register_1D_funcstr ("NLerr" , NL_error); /*----- discard the model array -----*/ #if 0 DESTROY_MODEL_ARRAY (model_array); #endif jump_on_NLfit_error = 0 ; return plint ; } /*************************************************************************** Main routine for this plugin (will be called from AFNI). If the return string is not NULL, some error transpired, and AFNI will popup the return string in a message box. ****************************************************************************/ char * NL_main( PLUGIN_interface * plint ) { char * str ; int ii, ival, ip ; float * tsar ; float min_constr, max_constr; MRI_IMAGE * im; /* pointer to image structures */ /*--------- go to first input line ---------*/ PLUTO_next_option(plint) ; plug_ignore = PLUTO_get_number(plint) ; plug_nrand = PLUTO_get_number(plint) ; plug_nbest = PLUTO_get_number(plint) ; /*------ loop over remaining options, check their tags, process them -----*/ do { str = PLUTO_get_optiontag(plint) ; if( str == NULL ) break ; if( strcmp(str,"Models") == 0 ) { str = PLUTO_get_string(plint) ; for (ii = 0; ii < num_noise_models; ii++) if (strcmp (str, noise_labels[ii]) == 0) plug_noise_index = ii; str = PLUTO_get_string(plint) ; for (ii = 0; ii < num_signal_models; ii++) if (strcmp (str, signal_labels[ii]) == 0) plug_signal_index = ii; str = PLUTO_get_string(plint); if (strcmp (str, "Absolute") == 0) plug_nabs = 1; else plug_nabs = 0; } else if( strcmp(str,"Noise") == 0 ) { ival = PLUTO_get_number(plint); min_constr = PLUTO_get_number(plint); max_constr = PLUTO_get_number(plint); if (min_constr > max_constr) return "**********************************\n" " Require min constr <= max constr \n" "**********************************" ; plug_min_nconstr[plug_noise_index][ival] = min_constr; plug_max_nconstr[plug_noise_index][ival] = max_constr; } else if( strcmp(str,"Signal") == 0 ) { ival = PLUTO_get_number(plint); min_constr = PLUTO_get_number(plint); max_constr = PLUTO_get_number(plint); if (min_constr > max_constr) return "**********************************\n" " Require min constr <= max constr \n" "**********************************" ; plug_min_sconstr[plug_signal_index][ival] = min_constr; plug_max_sconstr[plug_signal_index][ival] = max_constr; } else if( strcmp(str,"Time Scale") == 0 ) { str = PLUTO_get_string(plint); if (strcmp (str, "External") == 0){ plug_timeref = 1; str = PLUTO_get_string(plint); im = mri_read_1D (str); if (im == NULL) return "************************************\n" " Unable to read time reference file \n" "************************************" ; mri_free(im); strcpy (plug_tfilename, str); } else if( strcmp(str,"-inTR") == 0 ){ /* 22 July 1998 */ inTR = 1 ; plug_timeref = 0; } else { plug_timeref = 0; inTR = 0 ; /* 22 July 1998 */ } } else { return "************************\n" "Illegal optiontag found!\n" "************************" ; } } while(1) ; /*----- Identify software -----*/ printf ("\n\n"); printf ("Program: %s \n", PROGRAM_NAME); printf ("Author: %s \n", PROGRAM_AUTHOR); printf ("Date: %s \n", PROGRAM_DATE); printf ("\n"); /*----- show current input options -----*/ printf ("\nControls: \n"); printf ("Ignore = %5d \n", plug_ignore); printf ("Num Random = %5d \n", plug_nrand); printf ("Num Best = %5d \n", plug_nbest); printf ("Noise Constr = %s \n", constr_types[plug_nabs]); printf ("\nNoise Model = %s \n", noise_labels[plug_noise_index]); for (ip = 0; ip < plug_r[plug_noise_index]; ip++) { printf ("gn[%d]: min =%10.3f max =%10.3f %s \n", ip, plug_min_nconstr[plug_noise_index][ip], plug_max_nconstr[plug_noise_index][ip], noise_plabels[plug_noise_index][ip]); } printf ("\nSignal Model = %s \n", signal_labels[plug_signal_index]); for (ip = 0; ip < plug_p[plug_signal_index]; ip++) { printf ("gs[%d]: min =%10.3f max =%10.3f %s \n", ip, plug_min_sconstr[plug_signal_index][ip], plug_max_sconstr[plug_signal_index][ip], signal_plabels[plug_signal_index][ip]); } if (plug_timeref) printf ("\nExternal Time Reference = %s \n", plug_tfilename); else if( inTR ) printf ("\n-inTR Time Reference\n") ; else printf ("\nInternal Time Reference \n"); /*--- nothing left to do until data arrives ---*/ initialize = 1 ; /* force re-initialization */ return NULL ; } /*---------------------------------------------------------------------------*/ void NL_fitter( int nt , double to , double dt , float * vec, char ** label ) { NL_worker( nt , dt , vec , TRUE, label ) ; return ; } /*---------------------------------------------------------------------------*/ void NL_error( int nt , double to , double dt , float * vec, char ** label ) { NL_worker( nt , dt , vec , FALSE, label ) ; return ; } /*---------------------------------------------------------------------------*/ void NL_worker( int nt , double dt , float * vec , int dofit, char ** label ) { float * fit; int ii, nlen; float val; nlen = nt - plug_ignore; dsTR = dt ; /** find least squares fit coefficients **/ fit = nlfit (nlen, vec+plug_ignore, label); for (ii = 0; ii < plug_ignore; ii++) if (dofit) vec[ii] = fit[0]; else vec[ii] = vec[plug_ignore] - fit[0]; for (ii=plug_ignore; ii < nt; ii++) { if (dofit) vec[ii] = fit[ii-plug_ignore]; else vec[ii] = vec[ii] - fit[ii-plug_ignore] ; } free(fit) ; return ; }