00001 /*---------------------------------------------------------------------------*/
00002 /*
00003   This file contains routines for implementing the 3ddelay functions.
00004         Based on fim+.c  by B. Douglas Ward
00005 
00006   This software is copyrighted and owned by the Medical College of Wisconsin.
00007   See the file README.Copyright for details.
00008 
00009 */
00010 
00011 /*---------------------------------------------------------------------------*/
00012 /*
00013   Include software for linear regression analysis and sorting numbers.
00014 */
00015 
00016 #include "RegAna.c"
00017 #include "ranks.c"
00018 
00019 
00020 /*---------------------------------------------------------------------------*/
00021 
00022 /*---------------------------------------------------------------------------*/
00023 /* 
00024    Initialize independent variable X matrix 
00025 */
00026 
00027 int init_indep_var_matrix 
00028 (
00029   int q,                      /* number of parameters in the baseline model */
00030   int p,                      /* number of parameters in the baseline model 
00031                                  plus number of ideals */
00032   int NFirst,                 /* index of first image to use in the analysis */
00033   int N,                      /* total number of images used in the analysis */
00034   int polort,                 /* degree of polynomial ort */
00035   int num_ort_files,          /* number of ort time series files */
00036   int num_ideal_files,        /* number of ideal time series files */
00037   MRI_IMAGE ** ort_array,     /* ort time series arrays */
00038   int ** ort_list,            /* list of ort time series */
00039   MRI_IMAGE ** ideal_array,   /* ideal time series arrays */
00040   int ** ideal_list,          /* list of ideal time series */
00041   float * x_bot,              /* minimum of stimulus time series */
00042   float * x_ave,              /* average of stimulus time series */
00043   float * x_top,              /* maximum of stimulus time series */
00044   int * good_list,            /* list of good time points */
00045   matrix * x                  /* independent variable matrix */
00046 )
00047 
00048 {
00049   const int BIG_NUMBER = 33333;
00050   int i;                    /* time index */
00051   int m;                    /* X matrix column index */
00052   int n;                    /* X matrix row index */
00053   int is;                   /* input ideal index */
00054   float * far = NULL;
00055   int nx, ny, iq, nq;
00056   int Ngood;
00057   matrix xgood;
00058 
00059 
00060   /*----- Initialize X matrix -----*/
00061   matrix_create (N, p, x);
00062   matrix_initialize (&xgood);
00063 
00064 
00065   /*----- Set up columns of X matrix corresponding to polynomial orts -----*/
00066   for (m = 0;  m <= polort;  m++)
00067     for (n = 0;  n < N;  n++)
00068       {
00069         i = n + NFirst;
00070         (*x).elts[n][m] = pow ((double)i, (double)m);
00071       }
00072 
00073 
00074   /*----- Set up columns of X matrix corresponding to ort time series -----*/
00075   for (is = 0;  is < num_ort_files;  is++)
00076     {
00077       far = MRI_FLOAT_PTR (ort_array[is]);
00078       nx = ort_array[is]->nx;
00079       ny = ort_array[is]->ny;
00080 
00081       if (ort_list[is] == NULL)
00082         for (iq = 0;  iq < ny;  iq++)
00083           {
00084             for (n = 0;  n < N;  n++)
00085               {
00086                 i = n + NFirst;
00087                 (*x).elts[n][m] = *(far + iq*nx + i);
00088               }
00089             m++;
00090           }
00091       else
00092         {
00093           nq = ort_list[is][0];
00094           for (iq = 1;  iq <= nq;  iq++)
00095             {
00096               for (n = 0;  n < N;  n++)
00097                 {
00098                   i = n + NFirst;
00099                   (*x).elts[n][m] = *(far + ort_list[is][iq]*nx + i);
00100                 }
00101               m++;
00102             }
00103         }
00104     }
00105 
00106 
00107   /*----- Set up columns of X matrix corresponding to ideal time series -----*/
00108   for (is = 0;  is < num_ideal_files;  is++)
00109     {
00110       far = MRI_FLOAT_PTR (ideal_array[is]);
00111       nx = ideal_array[is]->nx;
00112       ny = ideal_array[is]->ny;
00113 
00114       if (ideal_list[is] == NULL)
00115         for (iq = 0;  iq < ny;  iq++)
00116           {
00117             for (n = 0;  n < N;  n++)
00118               {
00119                 i = n + NFirst;
00120                 (*x).elts[n][m] = *(far + iq*nx + i);
00121               }
00122             
00123             m++;
00124           }
00125       else
00126         {
00127           nq = ideal_list[is][0];
00128           for (iq = 1;  iq <= nq;  iq++)
00129             {
00130               for (n = 0;  n < N;  n++)
00131                 {
00132                   i = n + NFirst;
00133                   (*x).elts[n][m] = *(far + ideal_list[is][iq]*nx + i);
00134                 }
00135 
00136               m++;
00137             }
00138         }
00139     }
00140 
00141 
00142   /*----- Remove row if ideal contains value over 33333 -----*/
00143   Ngood = N;
00144   m = 0;
00145   for (n = 0;  n < N;  n++)
00146     {
00147       for (is = q;  is < p;  is++)
00148         {
00149           if ((*x).elts[n][is] >= BIG_NUMBER)  break;
00150         }
00151       if (is < p)
00152         {
00153           Ngood--;
00154         }
00155       else
00156         {
00157           good_list[m] = n;
00158           m++;
00159         }
00160     }
00161   matrix_extract_rows ((*x), Ngood, good_list, &xgood);
00162   matrix_equate (xgood, x);
00163 
00164 
00165   /*----- Find min, max, and ave for each column of the X matrix -----*/
00166   for (is = 0;  is < p;  is++)
00167     {      
00168       x_bot[is] = x_top[is] = (*x).elts[0][is];
00169       x_ave[is] = 0.0;
00170       for (n = 0;  n < Ngood;  n++)
00171         {
00172           if ((*x).elts[n][is] < x_bot[is])  x_bot[is] = (*x).elts[n][is];  
00173           if ((*x).elts[n][is] > x_top[is])  x_top[is] = (*x).elts[n][is];
00174           x_ave[is] += (*x).elts[n][is] / Ngood;
00175         }
00176     }
00177   
00178   
00179   matrix_destroy (&xgood);
00180 
00181   return (Ngood);
00182 
00183 }
00184 
00185 
00186 /*---------------------------------------------------------------------------*/
00187 /*
00188   Initialization for the delay analysis.
00189 */
00190 
00191 int init_delay 
00192 (
00193   int q,                      /* number of parameters in the baseline model */
00194   int p,                      /* number of parameters in the baseline model 
00195                                  plus number of ideals */
00196   int N,                      /* total number of images used in the analysis */
00197   int num_idealts,            /* number of ideal time series */
00198   matrix xdata,               /* independent variable matrix */
00199   matrix * x_base,            /* extracted X matrix    for baseline model */
00200   matrix * xtxinvxt_base,     /* matrix:  (1/(X'X))X'  for baseline model */
00201   matrix * x_ideal,           /* extracted X matrices  for ideal models */
00202   matrix * xtxinvxt_ideal,    /* matrix:  (1/(X'X))X'  for ideal models */
00203   float ** rarray             /* ranked arrays of ideal time series */
00204 )
00205 
00206 {
00207   int * plist = NULL;         /* list of model parameters */
00208   int ip, it;                 /* parameter indices */
00209   int is, js;                 /* ideal indices */ 
00210   int jm;                     /* lag index */
00211   int ok;                     /* flag for successful matrix calculation */
00212   matrix xtxinv_temp;         /* intermediate results */
00213   vector ideal;               /* ideal vector */
00214   vector coef_temp;           /* intermediate results */
00215   vector xres;                /* vector of residuals */
00216   float sse_base;             /* baseline model error sum of squares */ 
00217         
00218 
00219   /*----- Initialize matrix -----*/
00220   matrix_initialize (&xtxinv_temp);
00221   vector_initialize (&ideal);
00222   vector_initialize (&coef_temp);
00223   vector_initialize (&xres);
00224 
00225 
00226   /*----- Allocate memory -----*/
00227   plist = (int *) malloc (sizeof(int) * p);   MTEST (plist);
00228 
00229 
00230   /*----- Initialize matrices for the baseline model -----*/
00231   for (ip = 0;  ip < q;  ip++)
00232     plist[ip] = ip;
00233   ok = calc_matrices (xdata, q, plist, x_base, &xtxinv_temp, xtxinvxt_base);
00234   if (!ok)  { matrix_destroy (&xtxinv_temp);  return (0); };
00235 
00236 
00237   /*----- Initialize matrices for ideal functions -----*/
00238   for (is = 0;  is < num_idealts;  is++)
00239     {
00240       for (ip = 0;  ip < q;  ip++)
00241         {
00242           plist[ip] = ip;
00243         }
00244 
00245       plist[q] = q + is;
00246 
00247       ok = calc_matrices (xdata, q+1, plist, 
00248                           &(x_ideal[is]), &xtxinv_temp, &(xtxinvxt_ideal[is]));
00249       if (!ok)  { matrix_destroy (&xtxinv_temp);  return (0); };
00250     }
00251 
00252 
00253   /*----- Set up the ranked array for each ideal -----*/
00254   for (is = 0;  is < num_idealts;  is++)
00255     {
00256       /*----- Convert column of matrix to vector -----*/
00257       column_to_vector (xdata, q+is, &ideal);
00258 
00259       /*----- Calculate regression coefficients for baseline model -----*/
00260       calc_coef (*xtxinvxt_base, ideal, &coef_temp);
00261 
00262       /*----- Calculate the error sum of squares for the baseline model -----*/
00263       sse_base = calc_resids (*x_base, coef_temp, ideal, &xres);
00264     
00265       /*----- Form rank array from residual array -----*/
00266       rarray[is] = rank_darray (N, xres.elts);
00267 
00268     }
00269 
00270 
00271   /*----- Destroy matrix -----*/
00272   matrix_destroy (&xtxinv_temp);
00273   vector_destroy (&ideal);
00274   vector_destroy (&coef_temp);
00275   vector_destroy (&xres);
00276 
00277 
00278   /*----- Deallocate memory -----*/
00279   free (plist);   plist = NULL;
00280 
00281 
00282   return (1);
00283 }
00284 
00285 
00286 /*---------------------------------------------------------------------------*/
00287 
00288 
00289 /*---------------------------------------------------------------------------*/
00290 
00291 
00292 /*---------------------------------------------------------------------------*/
00293 /* 
00294    Calculate the sign function.
00295 */
00296 
00297 float sign (float x)
00298 
00299 {
00300   if (x > 0.0)  return (1.0);
00301   if (x < 0.0)  return (-1.0);
00302   return (0.0);
00303 }
00304 
00305 
00306 /*---------------------------------------------------------------------------*/
00307 
00308 int Read_part_file_delay (float *x,
00309                                         char *f_name,
00310                                         int a,
00311                                         int b)
00312    
00313     { 
00314      
00315      int cnt=0,ex,line_num;
00316      float buf;
00317      FILE*file_in;
00318      
00319      file_in = fopen (f_name,"r");
00320      if (file_in == NULL) {
00321             printf ("\aCould not open %s \n",f_name);
00322            printf ("Exiting @ Read_file function\n");
00323             exit (0);
00324          }
00325      
00326      if (a > b || a==0) {
00327                                 printf ("\a\n\33[1mError in (from , to) line numbers\n\33\[0m");
00328                                 printf ("Exiting @Read_part_file function \n");
00329                                 exit (0);
00330                    }
00331      
00332      line_num = 1;      
00333      if (a == 1) {
00334                         ex = fscanf (file_in,"%f",&x[cnt]);     
00335                         ++cnt;
00336                         }                                       
00337       else  ex = fscanf (file_in,"%f",&buf);                                            
00338      ++line_num;
00339      while (ex != EOF && line_num <= b)
00340       {
00341         if (line_num >=a && line_num <=b) 
00342         {
00343          ex = fscanf (file_in,"%f",&x[cnt]);
00344          ++cnt;
00345          if (ex == EOF) --cnt;
00346          }
00347         else 
00348         {
00349          ex = fscanf (file_in,"%f",&buf);
00350          }
00351         ++line_num;
00352         
00353       }
00354       
00355       if (ex == EOF) 
00356         {
00357             --line_num;
00358                 printf ("\n\33[1mEOF reached before line \33[0m%d\n",b);
00359                 printf ("Only %d lines were read, from line %d to %d\n",cnt,a,line_num-1);
00360         }
00361       
00362       fclose (file_in);
00363       return (cnt);                                                          
00364    }
00365 
00366 
00367 
00368 
00369 /*---------------------------------------------------------------------------*/
00370 
00371 /*---------------------------------------------------------------------------*/
00372   
00373 
00374 /*---------------------------------------------------------------------------*/
00375 /*---------------------------------------------------------------------------*/
00376 
00377 
00378 
00379