/*****************************************************************************
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 <math.h>
#include <stdlib.h>
#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 ;
}
syntax highlighted by Code2HTML, v. 0.9.1