01062 {
01063   const float EPSILON = 1.0e-10;
01064   int dimension;                 /* dimension of full model */
01065   int ip;                        /* parameter index */
01066   int df_rdcd;                   /* degrees of freedom for reduced model */
01067   int df_full;                   /* degrees of freedom for full model */
01068   float mse_full;                /* MSE for full model */
01069   float mse_reg;                 /* MSE for the regression */
01070   int ok;                        /* boolean for successful matrix inversion */
01071   float * y_array = NULL;        /* estimated signal time series */
01072   float * base_array = NULL;     /* baseline time series */
01073   float barea;                   /* area under baseline */
01074   int it;                        /* time series index */
01075   float y1, y2, y3;              /* temporary values for calculating area */
01076 
01077 
01078   /*----- initialization -----*/
01079   dimension = r + p;
01080   *rmsreg = *freg = *rsqr = *smax = *tmax = *pmax = *area = *parea = 0.0;
01081   for (ip = 0;  ip < dimension;  ip++)
01082     tpar_full[ip] = 0.0;
01083 
01084 
01085   /*----- check for insufficient variation in the data -----*/
01086   if (novar)  return;
01087 
01088 
01089   /*----- Adjust dof if constraints force a parameter to be a constant -----*/
01090   df_rdcd = ts_length - r;
01091   df_full = ts_length - dimension;
01092 
01093   for (ip = 0;  ip < r;  ip++)
01094     if (min_nconstr[ip] == max_nconstr[ip])
01095       {
01096         df_rdcd++;
01097         df_full++;
01098       }
01099 
01100   for (ip = 0;  ip < p;  ip++)
01101     if (min_sconstr[ip] == max_sconstr[ip])
01102       df_full++;
01103 
01104 
01105   /*----- MSE for full model -----*/
01106   mse_full = sse_full / df_full;
01107 
01108   /*----- MSE due to regression -----*/
01109   if (df_rdcd == df_full)
01110     mse_reg = 0.0;
01111   else
01112     mse_reg = (sse_rdcd - sse_full) / (df_rdcd - df_full);
01113   if (mse_reg < 0.0) mse_reg = 0.0;
01114 
01115   /*----- f-statistic for significance of the regression -----*/
01116   if (mse_full > EPSILON)
01117     *freg = mse_reg / mse_full;
01118   else
01119     *freg = 0.0;
01120 
01121   /*----- root-mean-square error for the regression -----*/
01122   *rmsreg = sqrt(mse_full);
01123 
01124   /*----- R^2 (coefficient of multiple determination) -----*/
01125   *rsqr = calc_rsqr (sse_full, sse_rdcd);
01126 
01127   /*----- generate time series corresponding to the signal model -----*/
01128   y_array = (float *) malloc (sizeof(float) * (ts_length));
01129   if (y_array == NULL)
01130     NLfit_error ("Unable to allocate memory for y_array");
01131 #if 0
01132   smodel (par_full+r, ts_length, x_array, y_array);
01133 #else
01134   AFNI_CALL_VOID_4ARG(smodel ,
01135                       float *,(par_full+r) , int,ts_length,
01136                       float **,x_array , float *,y_array ) ;
01137 #endif
01138 
01139   /*----- generate time series corresponding to the noise model -----*/
01140   base_array = (float *) malloc (sizeof(float) * (ts_length));
01141   if (base_array == NULL)
01142     NLfit_error ("Unable to allocate memory for base_array");
01143 #if 0
01144   nmodel (par_full, ts_length, x_array, base_array);
01145 #else
01146   AFNI_CALL_VOID_4ARG(nmodel ,
01147                       float *,par_full , int,ts_length,
01148                       float **,x_array , float *,base_array ) ;
01149 #endif
01150 
01151   /*----- initialize signal parameters -----*/
01152   *tmax = x_array[0][1];
01153   *smax = y_array[0];
01154   if (fabs(base_array[0]) > EPSILON)
01155     *pmax = 100.0 * y_array[0] / fabs(base_array[0]);
01156   else
01157     *pmax = 0.0;
01158   *area = 0.0;
01159   barea = 0.0;
01160   *parea = 0.0;
01161 
01162   /*----- calculate signed maximum of the signal, percent change, area -----*/
01163   for (it = 1;  it < ts_length;  it++)
01164     {
01165       /*----- calculate signed maximum of signal, and
01166         calculate max. percentage change of signal wrt baseline -----*/
01167       if (fabs(y_array[it]) > fabs(*smax))
01168         {
01169           *tmax = x_array[it][1];
01170           *smax = y_array[it];
01171           if (fabs(base_array[it]) > EPSILON)
01172             *pmax = 100.0 * y_array[it] / fabs(base_array[it]);
01173         }
01174       
01175       /*----- calculate area and percentage area under the signal -----*/
01176       y1 = y_array[it-1];   y2 = y_array[it];
01177       if ((y1 > 0.0) && (y2 > 0.0))
01178         {
01179           *area += (y1 + y2) / 2.0;
01180           *parea += (y1 + y2) / 2.0;
01181         }
01182       else if ((y1 < 0.0) && (y2 < 0.0))
01183         {
01184           *area -= (y1 + y2) / 2.0;
01185           *parea += (y1 + y2) / 2.0;
01186         }
01187       else
01188         {
01189           y3 = fabs(y1) + fabs(y2);
01190           if (y3 > EPSILON)
01191             {
01192               *area += (y1*y1 + y2*y2) / (2.0 * y3);
01193               if (y1 > y2)
01194                 *parea += (y1*y1 - y2*y2) / (2.0 * y3);
01195               else
01196                 *parea -= (y1*y1 - y2*y2) / (2.0 * y3);
01197             }
01198         }
01199       
01200       y1 = base_array[it-1];   y2 = base_array[it];
01201       if ((y1 > 0.0) && (y2 > 0.0))
01202         {
01203           barea += (y1 + y2) / 2.0;
01204         }
01205       else if ((y1 < 0.0) && (y2 < 0.0))
01206         {
01207           barea -= (y1 + y2) / 2.0;
01208         }
01209       else
01210         {
01211           y3 = fabs(y1) + fabs(y2);
01212           if (y3 > EPSILON)
01213             {
01214               barea += (y1*y1 + y2*y2) / (2.0 * y3);
01215             }
01216         }
01217     }  /* it */
01218 
01219   if (barea > EPSILON)
01220     *parea *= 100.0/barea;
01221   else
01222     *parea = 0.0;
01223 
01224   free (base_array);   base_array = NULL;
01225   free (y_array);      y_array = NULL;
01226 
01227 
01228   /*----- Calculate t-statistics? -----*/
01229   if (calc_tstats)
01230     {
01231       float stddev;                 /* est. std. dev. for a parameter */
01232       matrix d, dt, dtd, dtdinv;    /* matrices used for calc. of covariance */
01233 
01234       /*----- initialize matrices -----*/
01235       matrix_initialize (&d);
01236       matrix_initialize (&dt);
01237       matrix_initialize (&dtd);
01238       matrix_initialize (&dtdinv);
01239       
01240       /*----- calculate the matrix of partial derivatives D -----*/
01241       matrix_create (ts_length, dimension, &d);
01242       calc_partial_derivatives (nmodel, smodel, r, p, 
01243                                 min_nconstr, max_nconstr, 
01244                                 min_sconstr, max_sconstr, 
01245                                 ts_length, x_array, par_full, d);
01246       
01247       /*----- calculate variance-covariance matrix -----*/
01248       matrix_transpose (d, &dt);
01249       matrix_multiply (dt, d, &dtd);
01250       ok = matrix_inverse (dtd, &dtdinv);
01251       if (ok)
01252         for (ip = 0;  ip < dimension;  ip++)
01253           {
01254             stddev = sqrt((sse_full/(df_full)) * dtdinv.elts[ip][ip]);
01255             if (stddev > EPSILON)
01256               tpar_full[ip] = par_full[ip] / stddev;
01257             else
01258               tpar_full[ip] = 0.0;
01259           }
01260       else
01261         for (ip = 0;  ip < dimension;  ip++)
01262           {
01263             tpar_full[ip] = 0.0;
01264           }
01265       
01266       
01267       /*----- dispose of matrices -----*/
01268       matrix_destroy (&dtdinv);
01269       matrix_destroy (&dtd);
01270       matrix_destroy (&dt); 
01271       matrix_destroy (&d);
01272 
01273     }  /* if (calc_tstats) */    
01274 
01275 }

00350 {
00351   int i;                  /* parameter index */
00352 
00353 
00354   /*----- check noise model constraints -----*/
00355   if (nabs)   /*--- absolute noise parameter constraints ---*/
00356     for (i = 0;  i < r;  i++)
00357       {
00358         if (vertex[i] < min_nconstr[i])  return (1);
00359         if (vertex[i] > max_nconstr[i])  return (1);
00360       }
00361   else        /*--- relative noise parameter constraints ---*/
00362     for (i = 0;  i < r;  i++)
00363       {
00364         if (vertex[i] < min_nconstr[i] + b_array[i])  return (1);
00365         if (vertex[i] > max_nconstr[i] + b_array[i])  return (1);
00366       }
00367 
00368   /*----- check signal model constraints -----*/
00369   for (i = 0;  i < p;  i++)
00370     {
00371       if (vertex[i+r] < min_sconstr[i])  return (1);
00372       if (vertex[i+r] > max_sconstr[i])  return (1);
00373     }
00374 
00375   return (0);
00376 }

00843 {
00844   int iv;                        /* vector index */
00845   int ip;                        /* parameter index */ 
00846   float ** parameters = NULL;    /* array of parameter vectors */
00847   float * sse = NULL;            /* array of sse's */
00848 
00849 
00850   /*----- if this is the null signal model, 
00851           or if rms error for reduced model is very small, 
00852           just use the reduced model -----*/
00853   if ( (p < 1) || (sqrt(sse_rdcd/(ts_length - r)) < rms_min) ) 
00854     {
00855       *novar = 1;
00856       for (ip = 0;  ip < r;  ip++)
00857         par_full[ip] = par_rdcd[ip];
00858       for (ip = r;  ip < r+p;  ip++)
00859         par_full[ip] = 0.0;
00860       *sse_full = sse_rdcd;
00861       return;
00862     }
00863   else
00864     *novar = 0;
00865 
00866 
00867   /*----- initialize random number generator -----*/
00868   srand48 (1234567);
00869 
00870   /*----- allocate memory for calculation of full model-----*/
00871   initialize_full_model (r+p, nbest, &parameters, &sse); 
00872 
00873   /*----- evaluate randomly chosen vectors in the parameter space -----*/
00874   random_search (nmodel, smodel, r, p, nabs,
00875                  min_nconstr, max_nconstr, min_sconstr, max_sconstr,
00876                  ts_length, x_array, ts_array, par_rdcd, nrand, nbest, 
00877                  parameters, sse);
00878   
00879   /*----- use the best random vectors as the starting points for 
00880     nonlinear optimization -----*/
00881   for (iv = 0;  iv < nbest;  iv++)
00882     {
00883       simplex_optimization (nmodel, smodel, r, p, 
00884                             min_nconstr, max_nconstr, min_sconstr, max_sconstr,
00885                             nabs, ts_length, x_array, ts_array, par_rdcd,
00886                             parameters[iv], &sse[iv]);
00887     }
00888 
00889   /*----- save the best result from the nonlinear optimization -----*/      
00890   *sse_full = 1.0e+30;
00891   for (iv = 0;  iv < nbest;  iv++)
00892     {
00893       if (sse[iv] < *sse_full)
00894         {
00895           *sse_full = sse[iv];
00896           for (ip = 0;  ip < r+p;  ip++)
00897             par_full[ip] = parameters[iv][ip];
00898         }
00899     }
00900 
00901 
00902   /*----- release memory space -----*/
00903   for (iv = 0;  iv < nbest;  iv++)
00904     {
00905       free (parameters[iv]);
00906       parameters[iv] = NULL;
00907     }
00908   free (parameters);       parameters = NULL;     
00909   free (sse);              sse = NULL;
00910 
00911 }

00223 {
00224   matrix xt, xtx, xtxinv, xtxinvxt;     /* intermediate matrix results */
00225   vector yhat;                          /* estimate of observation vector */
00226   vector e;                             /* error in fit */
00227   int ok;                               /* successful matrix inversion? */
00228 
00229 
00230   /*----- initialize matrices -----*/
00231   matrix_initialize (&xt);
00232   matrix_initialize (&xtx);
00233   matrix_initialize (&xtxinv);
00234   matrix_initialize (&xtxinvxt);
00235   vector_initialize (&yhat);
00236   vector_initialize (&e);
00237 
00238   /*----- calculate least squares regression coefficients -----*/
00239   matrix_transpose (x, &xt);
00240   matrix_multiply (xt, x, &xtx);
00241   ok = matrix_inverse (xtx, &xtxinv);
00242   if (!ok)  NLfit_error ("Unable to invert matrix");
00243   matrix_multiply (xtxinv, xt, &xtxinvxt);
00244   vector_multiply (xtxinvxt, y, b);
00245 
00246   /*----- calculate least squares fit -----*/
00247   vector_multiply (x, *b, &yhat);
00248 
00249   /*----- calculate error sum of squares -----*/
00250   vector_subtract (y, yhat, &e);
00251   *sse = vector_dot (e, e);
00252 
00253   /*----- dispose of matrices -----*/
00254   vector_destroy (&e);
00255   vector_destroy (&yhat);
00256   matrix_destroy (&xtxinvxt);
00257   matrix_destroy (&xtxinv);
00258   matrix_destroy (&xtx);
00259   matrix_destroy (&xt);
00260 
00261 }

00936 {
00937   const float EPSILON = 1.0e-10;
00938   int dimension;          /* dimension of full model */
00939   int ip, jp;             /* parameter indices */ 
00940   int it;                 /* time index */
00941   float delp;             /* delta parameter */
00942   float * par = NULL;            /* perturbed parameter array */
00943   float * ts1_array = NULL;      /* estimated time series */
00944   float * ts2_array = NULL;      /* perturbed estimated time series */
00945 
00946 
00947   /*----- dimension of full model -----*/
00948   dimension = r + p;
00949 
00950   /*----- allocate memory -----*/
00951   ts1_array = (float *) malloc (sizeof(float) * ts_length);
00952   ts2_array = (float *) malloc (sizeof(float) * ts_length);
00953   par = (float *) malloc (sizeof(float) * dimension);
00954 
00955   /*----- fitted time series at estimated parameter vector -----*/
00956   full_model (nmodel, smodel, par_full, par_full + r, 
00957               ts_length, x_array, ts1_array);
00958 
00959   /*----- loop over parameters in model -----*/
00960   for (jp = 0;  jp < dimension;  jp++)
00961     {
00962       /*----- initialize parameters -----*/
00963       for (ip = 0;  ip < dimension;  ip++)
00964         par[ip] = par_full[ip];
00965 
00966       /*----- add small increment to the jpth parameter -----*/
00967       if (jp < r)
00968         delp = (max_nconstr[jp] - min_nconstr[jp]) / 1000.0;
00969       else
00970         delp = (max_sconstr[jp-r] - min_sconstr[jp-r]) / 1000.0;
00971       par[jp] += delp;
00972 
00973       /*----- fit time series for perturbed parameter -----*/
00974       full_model (nmodel, smodel, par, par + r, 
00975                   ts_length, x_array, ts2_array);
00976 
00977       /*----- estimate partial derivative -----*/
00978       if (delp > EPSILON)
00979         for (it = 0;  it < ts_length;  it++)
00980           d.elts[it][jp] = (ts2_array[it] - ts1_array[it]) / delp;
00981       else
00982         for (it = 0;  it < ts_length;  it++)
00983           d.elts[it][jp] = 0.0;
00984       
00985     }
00986 
00987   /*----- free memory -----*/
00988   free (par);        par = NULL;
00989   free (ts2_array);  ts2_array = NULL;
00990   free (ts1_array);  ts1_array = NULL;
00991 
00992 }  

00280 {
00281   matrix x;              /* independent variable matrix */
00282   vector y;              /* observation vector */
00283   vector b;              /* regression coefficient vector */
00284 
00285 
00286   /*----- initialize matrices and vectors -----*/
00287   matrix_initialize (&x);
00288   vector_initialize (&y);
00289   vector_initialize (&b);
00290 
00291   /*----- convert data to matrix and vector types -----*/
00292   array_to_matrix (n, r, x_array, &x);
00293   array_to_vector (n, y_array, &y);
00294 
00295   /*----- calculate linear regression -----*/
00296   calc_linear_regression (x, y, &b, sse);
00297 
00298   /*----- convert vector to array -----*/
00299   vector_to_array (b, b_array);
00300 
00301   /*----- dispose of matrices and vectors -----*/
00302   vector_destroy (&b);
00303   vector_destroy (&y);
00304   matrix_destroy (&x);
00305 }

01006 {
01007   const float EPSILON = 1.0e-5;     /* protection against divide by zero */
01008   float rsqr;                       /* coeff. of multiple determination R^2  */
01009 
01010 
01011   /*----- coefficient of multiple determination R^2 -----*/
01012   if (sser < EPSILON)
01013     rsqr = 0.0;
01014   else
01015     rsqr = (sser - ssef) / sser;
01016 
01017 
01018   /*----- Limit range of values for R^2 -----*/
01019   if (rsqr < 0.0)   rsqr = 0.0;
01020   if (rsqr > 1.0)   rsqr = 1.0;
01021 
01022 
01023   /*----- Return coefficient of multiple determination R^2 -----*/
01024   return (rsqr);
01025 }

00614 {
00615   const float BIG_NUMBER = 1.0e+10;   /* return large number if constraints
00616                                          are violated */
00617   int i;                              /* time point index */
00618   float t;                            /* time */
00619   float diff;                         /* error at individual time point */
00620   float sse;                          /* error sum of squares */
00621   float * y_array = NULL;             /* estimated time series */
00622 
00623 
00624   /*----- allocate memory for estimated time series -----*/
00625   y_array = (float *) malloc (sizeof(float) * ts_length);
00626 
00627   /*----- apply constraints? -----*/
00628   if (calc_constraints (r, p, nabs, b_array, min_nconstr, max_nconstr, 
00629                         min_sconstr, max_sconstr, vertex)) 
00630     { 
00631       free (y_array);   y_array = NULL;
00632       return (BIG_NUMBER);
00633     }    
00634   
00635   /*----- create estimated time series using the full model parameters -----*/
00636   full_model (nmodel, smodel, vertex, vertex + r, ts_length, x_array, y_array);
00637 
00638   /*----- calculate error sum of squares -----*/
00639   sse = 0.0;
00640   for (i = 0;  i < ts_length;  i++)
00641     {
00642       diff = ts_array[i] - y_array[i];
00643       sse += diff * diff;
00644     }
00645 
00646   /*----- release memory -----*/
00647   free (y_array);   y_array = NULL;
00648 
00649   /*----- return error sum of squares -----*/
00650   return (sse);
00651 }

01286 {
01287    int which , status ;
01288    double p , q , f , dfn , dfd , bound ;
01289 
01290    which  = 1 ;
01291    p      = 0.0 ;
01292    q      = 0.0 ;
01293    f      = ff ;
01294    dfn    = dofnum ;
01295    dfd    = dofden ;
01296 
01297    cdff( &which , &p , &q , &f , &dfn , &dfd , &status , &bound ) ;
01298 
01299    if( status == 0 ) return q ;
01300    else              return 1.0 ;
01301 }

00538 {
00539   int it;                     /* time index */
00540   float * y_array = NULL;     /* estimated signal time series */
00541 #ifdef SAVE_RAN
00542   int use_model = ( RAN_sind < 0 || ts_length != OLD_ts_length ) ;
00543 #endif
00544 
00545   /*----- generate time series corresponding to signal model -----*/
00546 
00547 #ifdef SAVE_RAN
00548   if( use_model ){  /* don't use saved time series */
00549 #endif
00550      y_array = (float *) malloc (sizeof(float) * (ts_length));
00551      if (y_array == NULL)
00552        NLfit_error ("Unable to allocate memory for y_array");
00553 #if 0
00554      smodel (gs, ts_length, x_array, y_array);
00555 #else
00556      AFNI_CALL_VOID_4ARG(smodel ,
00557                          float *,gs , int,ts_length,
00558                          float **,x_array , float *,y_array ) ;
00559 #endif
00560 
00561 #ifdef SAVE_RAN
00562   } else            /* recall a saved time series */
00563      y_array = RAN_sts + (ts_length*RAN_sind) ;
00564 #endif
00565 
00566   /*----- generate time series corresponding to the noise model -----*/
00567 #if 0
00568   nmodel (gn, ts_length, x_array, yhat_array);
00569 #else
00570   AFNI_CALL_VOID_4ARG(nmodel ,
00571                       float *,gn , int,ts_length,
00572                       float **,x_array , float *,yhat_array ) ;
00573 #endif
00574 
00575   /*----- add signal and noise model time series -----*/
00576   for (it = 0;  it < ts_length;  it++)
00577     yhat_array[it] += y_array[it];
00578   
00579 
00580   /*----- deallocate memory -----*/
00581 #ifdef SAVE_RAN
00582   if( use_model )
00583 #endif
00584     free (y_array) ;
00585 
00586   y_array = NULL;  /* not really needed since this array is 'auto' */
00587 }

00321 {
00322   int i;                    /* parameter vector index */
00323 
00324   *sse = (float *) malloc (sizeof(float) * nbest);   
00325   *parameters = (float **) malloc (sizeof(float *) * nbest);
00326   for (i = 0;  i < nbest;  i++)
00327     (*parameters)[i] = (float *) malloc (sizeof(float) * dimension);
00328        
00329 }

00159 {
00160   int im, ip, index;
00161   char message[MAX_NAME_LENGTH];    /* error message */
00162 
00163 
00164   index = -1;
00165   for (im = 0;  im < model_array->num;  im++)
00166     if (strncmp(model_array->modar[im]->interface->label, 
00167                 nname, MAX_NAME_LENGTH) == 0)  index = im;
00168   if (index < 0)
00169     {
00170       sprintf (message, "Unable to locate noise model %s", nname);
00171       NLfit_error (message);
00172     }
00173 
00174   if (model_array->modar[index]->interface->model_type != MODEL_NOISE_TYPE)
00175     {
00176       printf ("type = %d \n", 
00177               model_array->modar[index]->interface->model_type);
00178       sprintf (message, "%s has not been declared a noise model", nname);
00179       NLfit_error (message);
00180     }
00181 
00182   *nmodel = model_array->modar[index]->interface->call_func;
00183   if (*nmodel == NULL) 
00184     {
00185       sprintf (message, "Noise model %s not properly implemented", nname);
00186       NLfit_error (message);
00187     }
00188 
00189   *r = model_array->modar[index]->interface->params;
00190   if ((*r < 0) || (*r > MAX_PARAMETERS)) 
00191     {
00192       sprintf (message, "Illegal number of parameters for noise model %s", 
00193                nname);
00194       NLfit_error (message);
00195     }
00196 
00197   for (ip = 0;  ip < *r;  ip++)
00198     {
00199       strncpy (npname[ip], model_array->modar[index]->interface->plabel[ip],
00200                MAX_NAME_LENGTH);
00201       min_nconstr[ip] = model_array->modar[index]->interface->min_constr[ip]; 
00202       max_nconstr[ip] = model_array->modar[index]->interface->max_constr[ip]; 
00203       if (min_nconstr[ip] > max_nconstr[ip])
00204         NLfit_error("Must have noise parameter min cnstrnts <= max cnstrnts");
00205     }
00206 
00207 }

00091 {
00092   int im, ip, index;
00093   char message[MAX_NAME_LENGTH];    /* error message */
00094 
00095 
00096   index = -1;
00097   for (im = 0;  im < model_array->num;  im++)
00098     if (strncmp(model_array->modar[im]->interface->label, 
00099                 sname, MAX_NAME_LENGTH) == 0)  index = im;
00100   if (index < 0)
00101     {
00102       sprintf (message, "Unable to locate signal model %s", sname);
00103       NLfit_error (message);
00104     }
00105 
00106   if (model_array->modar[index]->interface->model_type != MODEL_SIGNAL_TYPE)
00107     {
00108       sprintf (message, "%s has not been declared a signal model", sname);
00109       NLfit_error (message);
00110     }
00111 
00112   *smodel = model_array->modar[index]->interface->call_func;
00113   if (*smodel == NULL) 
00114     {
00115       sprintf (message, "Signal model %s not properly implemented", sname);
00116       NLfit_error (message);
00117     }
00118 
00119   *p = model_array->modar[index]->interface->params;
00120   if ((*p < 0) || (*p > MAX_PARAMETERS)) 
00121     {
00122       sprintf (message, "Illegal number of parameters for signal model %s", 
00123                sname);
00124       NLfit_error (message);
00125     }
00126 
00127   for (ip = 0;  ip < *p;  ip++)
00128     {
00129       strncpy (spname[ip], model_array->modar[index]->interface->plabel[ip],
00130                MAX_NAME_LENGTH);
00131       min_sconstr[ip] = model_array->modar[index]->interface->min_constr[ip]; 
00132       max_sconstr[ip] = model_array->modar[index]->interface->max_constr[ip]; 
00133       if (min_sconstr[ip] > max_sconstr[ip])
00134         NLfit_error("Must have signal parameter min cnstrnts <= max cnstrnts");
00135     }
00136 
00137 }

00065 {
00066    fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00067    if( jump_on_NLfit_error ) longjmp( NLfit_error_jmpbuf , 1 ) ;  /* 01 May 2003 */
00068    exit(1);
00069 }

00413 {
00414    int ii ;
00415    if( n < 1 || a == NULL || b == NULL ) return 1 ;
00416    for( ii=0 ; ii < n ; ii++ ) if( a[ii] != b[ii] ) return 1 ;
00417    return 0 ;
00418 }

00439 {
00440    int ip , ipt ;
00441    float * par ;
00442    float * ts ;
00443 
00444    /*** check if the signal model is dead on arrival ***/
00445 
00446    if( smodel == NULL ){
00447       OLD_smodel = NULL ;
00448       OLD_p      = -1 ;
00449       if( OLD_min_sconstr != NULL ){ free(OLD_min_sconstr); OLD_min_sconstr = NULL; }
00450       if( OLD_max_sconstr != NULL ){ free(OLD_max_sconstr); OLD_max_sconstr = NULL; }
00451       if( RAN_spar        != NULL ){ free(RAN_spar)       ; RAN_spar        = NULL; }
00452       if( RAN_sts         != NULL ){ free(RAN_sts )       ; RAN_sts         = NULL; }
00453       return ;
00454    }
00455 
00456    /*** check if the signal model has changed ***/
00457 
00458    if( smodel        != OLD_smodel                     ||
00459        p             != OLD_p                          ||
00460        ts_length     != OLD_ts_length                  ||
00461        nrand         != OLD_nrand                      ||
00462        x_array[0][1] != OLD_x0                         ||
00463        x_array[1][1] != OLD_x1                         ||
00464        RAN_compare_vect(p,min_sconstr,OLD_min_sconstr) ||
00465        RAN_compare_vect(p,max_sconstr,OLD_max_sconstr)   ){
00466 
00467       /* save parameters of new signal model */
00468 
00469 #if 0
00470       fprintf(stderr,"++ NLfit: initializing random signal models") ;
00471 #endif
00472 
00473       OLD_smodel    = smodel ;
00474       OLD_p         = p ;
00475       OLD_ts_length = ts_length ;
00476       OLD_nrand     = nrand ;
00477       OLD_x0        = x_array[0][1] ;
00478       OLD_x1        = x_array[1][1] ;
00479       if( OLD_min_sconstr != NULL ) free(OLD_min_sconstr) ;
00480       if( OLD_max_sconstr != NULL ) free(OLD_max_sconstr) ;
00481       OLD_min_sconstr = (float *) malloc(sizeof(float)*p) ;
00482       OLD_max_sconstr = (float *) malloc(sizeof(float)*p) ;
00483       memcpy( OLD_min_sconstr , min_sconstr , sizeof(float)*p ) ;
00484       memcpy( OLD_max_sconstr , max_sconstr , sizeof(float)*p ) ;
00485 
00486       /* create space for new random signal model vectors */
00487 
00488       if( RAN_spar != NULL ) free(RAN_spar) ;
00489       if( RAN_sts  != NULL ) free(RAN_sts ) ;
00490       RAN_spar = (float *) malloc(sizeof(float)*nrand*p) ;
00491       RAN_sts  = (float *) malloc(sizeof(float)*nrand*ts_length) ;
00492 
00493       /* compute new random signal vectors */
00494 
00495       for( ipt=0 ; ipt < nrand ; ipt++ ){
00496 
00497          par = RAN_spar + (ipt*p) ;         /* ipt-th parameter vector */
00498          ts  = RAN_sts  + (ipt*ts_length) ; /* ipt-th signal model time series */
00499 
00500          for( ip=0 ; ip < p ; ip++ )                 /* parameter vector */
00501             par[ip] = get_random_value(min_sconstr[ip], max_sconstr[ip]) ;
00502 
00503 #if 0
00504          smodel( par , ts_length , x_array , ts ) ;  /* time series vector */
00505 #else
00506          AFNI_CALL_VOID_4ARG(smodel ,
00507                              float *,par , int,ts_length,
00508                              float **,x_array , float *,ts ) ;
00509 #endif
00510       }
00511 
00512 #if 0
00513       fprintf(stderr," - done\n") ;
00514 #endif
00515    } /* end of signal model stowage */
00516 
00517    return ;
00518 }

00680 {
00681   int ip;                 /* parameter index */
00682   int iv, ipt;            /* vector indices */
00683   float sse;              /* error sum of squares */
00684   float * par = NULL;     /* test parameter vector */
00685 
00686 
00687 #ifdef SAVE_RAN
00688   /*** 28 July 1998: set up to choose the best noise model for each case,
00689                      rather than a completely random one -- RWCox.       ***/
00690 #undef  BEST_NOISE
00691 #ifdef  BEST_NOISE
00692    float * qts = (float *) malloc(sizeof(float)*ts_length) ;
00693    float junk ;
00694 #endif
00695 #endif /* SAVE_RAN */
00696 
00697   /*** 27 July 1998: generate some random signal parameters,
00698                      but only if the signal model is new this time. ***/
00699 #ifdef SAVE_RAN
00700    RAN_setup( nmodel , smodel , r , p , nabs ,
00701               min_nconstr, max_nconstr,
00702               min_sconstr, max_sconstr, 
00703               ts_length, x_array, nrand ) ;
00704 #endif
00705 
00706   /*----- allocate memory for test parameter vector -----*/
00707   par = (float *) malloc (sizeof(float) * (r+p));
00708 
00709   /*----- initialize response values -----*/
00710   for (iv = 0;  iv < nbest;  iv++)
00711     response[iv] = 1.0e+30;
00712 
00713 
00714   /*----- loop over random vectors -----*/
00715   for (ipt = 0;  ipt < nrand;  ipt++)
00716     {
00717 
00718       /*** 28 July 1998: get the best noise parameter model ***/
00719 #ifdef BEST_NOISE
00720       for( ip=0 ; ip < ts_length ; ip++ )  /* subtract signal model from data */
00721          qts[ip] = ts_array[ip] - RAN_sts[ts_length*ipt+ip] ;
00722       calc_reduced_model( ts_length , r , x_array , qts , par , &junk ) ;
00723 #else
00724       /*----- select random parameters -----*/
00725       if (nabs)   /*--- absolute noise parameter constraints ---*/
00726         for (ip = 0;  ip < r;  ip++)
00727           par[ip] = get_random_value (min_nconstr[ip], max_nconstr[ip]);
00728       else        /*--- relative noise parameter constraints ---*/
00729         for (ip = 0;  ip < r;  ip++)
00730           par[ip] = get_random_value (min_nconstr[ip] + par_rdcd[ip], 
00731                                       max_nconstr[ip] + par_rdcd[ip]);
00732 #endif /* BEST_NOISE */
00733 
00734       /*** 27 July 1998: get the signal parameters from the saved set ***/
00735 #ifdef SAVE_RAN
00736       for( ip=0 ; ip < p ; ip++ )
00737          par[ip+r] = RAN_spar[ipt*p+ip] ;
00738 #else
00739       for (ip = 0;  ip < p;  ip++)
00740           par[ip+r] = get_random_value (min_sconstr[ip], max_sconstr[ip]);
00741 #endif
00742 
00743 
00744       /*----- evaluate this random vector -----*/
00745 #ifdef SAVE_RAN
00746       RAN_sind = ipt ;
00747 #endif
00748 
00749       sse = calc_sse (nmodel, smodel, r, p, nabs, min_nconstr, max_nconstr, 
00750                       min_sconstr, max_sconstr, par_rdcd, par, 
00751                       ts_length, x_array, ts_array);
00752 
00753       /*----- save best random vectors -----*/
00754 
00755       /*** The previous code for this is erroneous, since it does not
00756            allow for the case where an in-between best and worst vector
00757            is found.  For example, suppose the response array is now
00758            [ 1 3 5 ] (with nbest==3).  Then we get sse=2.  This would
00759            replace the 3 value, leaving us with [1 2 5], instead of
00760            [ 1 2 3 ], which is a better set.  The solution is to push
00761            the displaced vectors up in the array.  In this way, the
00762            response array is always sorted, and keeps the truly best. ***/
00763 
00764 #define PUSH_BEST  /* 27 July 1998 -- RWCox */
00765 #ifdef  PUSH_BEST
00766       for (iv = 0;  iv < nbest;  iv++)
00767         if (sse < response[iv]){
00768            int jv ;
00769            for( jv=nbest-1 ; jv > iv ; jv-- ){  /* the push-up code */
00770               response[jv] = response[jv-1] ;
00771               for( ip=0 ; ip < r+p ; ip++ )
00772                  parameters[jv][ip] = parameters[jv-1][ip] ;
00773            }
00774 
00775            response[iv] = sse ;                 /* now can copy new */
00776            for( ip=0 ; ip < r+p ;  ip++ )       /* values into place */
00777               parameters[iv][ip] = par[ip] ;
00778            break ;
00779         }
00780 #else
00781       for (iv = 0;  iv < nbest;  iv++)          /* this is the old code */
00782         if (sse < response[iv])
00783           {
00784             for (ip = 0;  ip < r+p;  ip++)
00785               parameters[iv][ip] = par[ip];
00786             response[iv] = sse;
00787             break;
00788           }
00789 #endif
00790   } /* end of loop over random vectors */
00791 
00792   /*----- release memory space -----*/
00793   free (par);       par = NULL;
00794 
00795 #ifdef SAVE_RAN
00796       RAN_sind = -1 ;
00797 #endif
00798 
00799 #if 0
00800    /*** Some debugging stuff -- 28 July 1998 ***/
00801    { static count=-1 ;
00802      count++ ;
00803      if( count%10 == 0 ){
00804         printf("Random search sse:\n") ;
00805         for( iv=0 ; iv < nbest ; iv++ ) printf(" %g",response[iv]) ;
00806         printf("\n") ;
00807      }
00808    }
00809 #endif
00810  
00811 }

01333 {
01334   int ip;                  /* parameter index */
01335   double pvalue;
01336 
01337 
01338   /*----- create label if desired by calling program -----*/
01339   if (label == NULL)  return;
01340 
01341   /*----- make this a 0 length string to start -----*/
01342   lbuf[0] = '\0';
01343       
01344   /*----- write reduced model parameters -----*/
01345   sprintf (sbuf, "Reduced (%s) Model: \n", nname);
01346   strcat (lbuf, sbuf);
01347   for (ip = 0;  ip < r;  ip++)
01348     {
01349       sprintf (sbuf, "b[%d]= %12.6f  %s \n", ip, par_rdcd[ip], npname[ip]); 
01350       strcat (lbuf, sbuf);
01351     }
01352     
01353   /*----- write full model parameters -----*/  
01354   sprintf (sbuf, "\nFull (%s + %s) Model: \n", nname, sname);
01355   strcat (lbuf, sbuf);
01356   for (ip = 0;  ip < r;  ip++)
01357     {
01358       sprintf (sbuf, "gn[%d]=%12.6f  %s \n", ip, par_full[ip], npname[ip]);
01359       strcat (lbuf, sbuf);
01360                     
01361     }
01362   for (ip = 0;  ip < p;  ip++)
01363     {
01364       sprintf (sbuf, "gs[%d]=%12.6f  %s \n", ip, par_full[ip+r], spname[ip]);
01365       strcat (lbuf, sbuf);            
01366     }
01367   
01368   sprintf (sbuf, "\nSignal Tmax  = %12.3f \n", tmax);
01369   strcat (lbuf, sbuf);
01370   sprintf (sbuf,   "Signal Smax  = %12.3f \n", smax);
01371   strcat (lbuf, sbuf);
01372   sprintf (sbuf,   "Signal PSmax = %12.3f \n", pmax);
01373   strcat (lbuf, sbuf);
01374   sprintf (sbuf,   "Signal Area  = %12.3f \n", area);
01375   strcat (lbuf, sbuf);
01376   sprintf (sbuf,   "Signal PArea = %12.3f \n", parea);
01377   strcat (lbuf, sbuf);
01378 
01379   sprintf (sbuf, "\nRMSE Rdcd = %8.3f \n", sqrt(sse_rdcd/(ts_length-r)));
01380   strcat (lbuf, sbuf);
01381   sprintf (sbuf, "RMSE Full = %8.3f \n", sqrt(sse_full/(ts_length-r-p)));
01382   strcat (lbuf, sbuf);
01383             
01384   sprintf (sbuf, "\nR^2       = %7.3f \n", rsqr);
01385   strcat (lbuf, sbuf);
01386   sprintf (sbuf, "F[%2d,%3d] = %7.3f \n", p, ts_length-r-p, freg);
01387   strcat (lbuf, sbuf);
01388   pvalue = fstat_t2p ( (double) freg, (double) p, (double) ts_length-r-p);
01389   sprintf (sbuf, "p-value   = %e  \n", pvalue);
01390   strcat (lbuf, sbuf);
01391  
01392 
01393   /*----- send address of lbuf back in what label points to -----*/
01394   *label = lbuf;
01395   
01396 }