Skip to content

AFNI/NIfTI Server

Sections
Personal tools
You are here: Home » AFNI » Documentation

Doxygen Source Code Documentation


Main Page   Alphabetical List   Data Structures   File List   Data Fields   Globals   Search  

fim+.c File Reference

#include "RegAna.c"
#include "ranks.c"

Go to the source code of this file.


Defines

#define MAX_OUTPUT_TYPE   12
#define FIM_FitCoef   (0)
#define FIM_BestIndex   (1)
#define FIM_PrcntChange   (2)
#define FIM_Baseline   (4)
#define FIM_Correlation   (6)
#define FIM_PrcntFromAve   (3)
#define FIM_Average   (5)
#define FIM_PrcntFromTop   (7)
#define FIM_Topline   (8)
#define FIM_SigmaResid   (9)
#define FIM_SpearmanCC   (10)
#define FIM_QuadrantCC   (11)
#define USE_LEGENDRE

Functions

double legendre (double x, int m)
int init_indep_var_matrix (int q, int p, int NFirst, int N, int polort, int num_ort_files, int num_ideal_files, MRI_IMAGE **ort_array, int **ort_list, MRI_IMAGE **ideal_array, int **ideal_list, float *x_bot, float *x_ave, float *x_top, int *good_list, matrix *x)
int init_regression_analysis (int q, int p, int N, int num_idealts, matrix xdata, matrix *x_base, matrix *xtxinvxt_base, matrix *x_ideal, matrix *xtxinvxt_ideal, float **rarray)
float calc_CC (int N, float *x, float *y)
float calc_SpearmanCC (int N, float *r, float *s)
float sign (float x)
float calc_QuadrantCC (int N, float *r, float *s)
float percent_change (float num, float den)
void regression_analysis (int N, int q, int num_idealts, matrix x_base, matrix xtxinvxt_base, matrix *x_ideal, matrix *xtxinvxt_ideal, vector y, float *x_bot, float *x_ave, float *x_top, float **rarray, int *output_type, float *FimParams)
void report_results (int *output_type, float *FimParams, char **label)

Variables

char * OUTPUT_TYPE_name [MAX_OUTPUT_TYPE]
char lbuf [4096]
char sbuf [256]

Define Documentation

#define FIM_Average   (5)
 

Definition at line 48 of file fim+.c.

Referenced by display_help_menu(), and regression_analysis().

#define FIM_Baseline   (4)
 

Definition at line 45 of file fim+.c.

Referenced by display_help_menu(), and regression_analysis().

#define FIM_BestIndex   (1)
 

Definition at line 43 of file fim+.c.

Referenced by check_for_valid_inputs(), display_help_menu(), and regression_analysis().

#define FIM_Correlation   (6)
 

Definition at line 46 of file fim+.c.

Referenced by calculate_results(), display_help_menu(), regression_analysis(), and write_bucket_data().

#define FIM_FitCoef   (0)
 

Definition at line 42 of file fim+.c.

Referenced by display_help_menu(), and regression_analysis().

#define FIM_PrcntChange   (2)
 

Definition at line 44 of file fim+.c.

Referenced by display_help_menu(), and regression_analysis().

#define FIM_PrcntFromAve   (3)
 

Definition at line 47 of file fim+.c.

Referenced by display_help_menu(), and regression_analysis().

#define FIM_PrcntFromTop   (7)
 

Definition at line 49 of file fim+.c.

Referenced by display_help_menu(), and regression_analysis().

#define FIM_QuadrantCC   (11)
 

Definition at line 53 of file fim+.c.

Referenced by display_help_menu(), regression_analysis(), and write_bucket_data().

#define FIM_SigmaResid   (9)
 

Definition at line 51 of file fim+.c.

Referenced by display_help_menu(), and regression_analysis().

#define FIM_SpearmanCC   (10)
 

Definition at line 52 of file fim+.c.

Referenced by display_help_menu(), regression_analysis(), and write_bucket_data().

#define FIM_Topline   (8)
 

Definition at line 50 of file fim+.c.

Referenced by display_help_menu(), and regression_analysis().

#define MAX_OUTPUT_TYPE   12
 

Definition at line 34 of file fim+.c.

Referenced by allocate_memory(), calculate_results(), get_options(), initialize_options(), report_results(), save_voxel(), terminate_program(), and write_bucket_data().

#define USE_LEGENDRE
 

Definition at line 57 of file fim+.c.


Function Documentation

float calc_CC int    N,
float *    x,
float *    y
 

Definition at line 462 of file fim+.c.

References i.

Referenced by calc_QuadrantCC(), and calc_SpearmanCC().

00468 {
00469   const float EPSILON = 1.0e-10;   /* protect against divide by zero */
00470   float cc;
00471   int i;
00472   float csum, xsum, ysum, prod;
00473 
00474 
00475   csum = xsum = ysum = 0.0;
00476 
00477   for (i = 0;  i < N;  i++)
00478     {
00479       csum += x[i] * y[i];
00480       xsum += x[i] * x[i];
00481       ysum += y[i] * y[i];
00482     }
00483 
00484 
00485   prod = xsum * ysum;
00486   if (prod < EPSILON)
00487     cc = 0.0;
00488   else
00489     cc = csum / sqrt(prod);
00490 
00491 
00492   return (cc);
00493  
00494 }

float calc_QuadrantCC int    N,
float *    r,
float *    s
 

Definition at line 558 of file fim+.c.

References calc_CC(), free, i, malloc, r, and sign.

Referenced by regression_analysis().

00564 {
00565   float cc;
00566   float rbar;
00567   float * dr = NULL;
00568   float * ds = NULL;
00569   int i;
00570 
00571   
00572   dr = (float *) malloc (sizeof(float) * N);
00573   ds = (float *) malloc (sizeof(float) * N);
00574 
00575 
00576   rbar = (N+1.0) / 2.0;
00577   for (i = 0;  i < N;  i++)
00578     {
00579       dr[i] = sign(r[i] - rbar);
00580       ds[i] = sign(s[i] - rbar);
00581     }
00582 
00583 
00584   cc = calc_CC(N, dr, ds);
00585 
00586   free (dr);   dr = NULL;
00587   free (ds);   ds = NULL;
00588  
00589   return (cc);
00590 }

float calc_SpearmanCC int    N,
float *    r,
float *    s
 

Definition at line 503 of file fim+.c.

References calc_CC(), free, i, malloc, and r.

Referenced by regression_analysis().

00509 {
00510   float cc;
00511   float rbar;
00512   float * dr = NULL;
00513   float * ds = NULL;
00514   int i;
00515 
00516   
00517   dr = (float *) malloc (sizeof(float) * N);
00518   ds = (float *) malloc (sizeof(float) * N);
00519 
00520 
00521   rbar = (N+1.0) / 2.0;
00522   for (i = 0;  i < N;  i++)
00523     {
00524       dr[i] = r[i] - rbar;
00525       ds[i] = s[i] - rbar;
00526     }
00527 
00528 
00529   cc = calc_CC(N, dr, ds);
00530 
00531   free (dr);   dr = NULL;
00532   free (ds);   ds = NULL;
00533  
00534   return (cc);
00535 }

int init_indep_var_matrix int    q,
int    p,
int    NFirst,
int    N,
int    polort,
int    num_ort_files,
int    num_ideal_files,
MRI_IMAGE **    ort_array,
int **    ort_list,
MRI_IMAGE **    ideal_array,
int **    ideal_list,
float *    x_bot,
float *    x_ave,
float *    x_top,
int *    good_list,
matrix   x
 

if here, m > 20 ==> use recurrence relation *

Definition at line 185 of file fim+.c.

References BIG_NUMBER, far, i, legendre(), matrix_create(), matrix_destroy(), matrix_equate(), matrix_extract_rows(), matrix_initialize(), MRI_FLOAT_PTR, p, and q.

Referenced by calculate_results().

00205 {
00206   const int BIG_NUMBER = 33333;
00207   int i;                    /* time index */
00208   int m;                    /* X matrix column index */
00209   int n;                    /* X matrix row index */
00210   int is;                   /* input ideal index */
00211   float * far = NULL;
00212   int nx, ny, iq, nq;
00213   int Ngood;
00214   matrix xgood;
00215 #ifdef USE_LEGENDRE
00216   double nfac,nsub ;
00217 #endif
00218 
00219 
00220   /*----- Initialize X matrix -----*/
00221   matrix_create (N, p, x);
00222   matrix_initialize (&xgood);
00223 
00224 #ifdef USE_LEGENDRE
00225   nsub = 0.5*(N-1) ;
00226   nfac = 1.0/nsub  ;
00227 #endif
00228 
00229   /*----- Set up columns of X matrix corresponding to polynomial orts -----*/
00230   for (m = 0;  m <= polort;  m++)
00231     for (n = 0;  n < N;  n++)
00232       {
00233 #ifndef USE_LEGENDRE   /** the old polort way: t^m **/
00234              i = n + NFirst;
00235              (*x).elts[n][m] = pow ((double)i, (double)m);
00236 
00237 #else                  /** the new polort way: Legendre P_m(t) - 29 Mar 2005 **/
00238 
00239              (*x).elts[n][m] = legendre( nfac*(n-nsub) , m ) ;
00240 #endif
00241       }
00242 
00243 
00244   /*----- Set up columns of X matrix corresponding to ort time series -----*/
00245   for (is = 0;  is < num_ort_files;  is++)
00246     {
00247       far = MRI_FLOAT_PTR (ort_array[is]);
00248       nx = ort_array[is]->nx;
00249       ny = ort_array[is]->ny;
00250 
00251       if (ort_list[is] == NULL)
00252         for (iq = 0;  iq < ny;  iq++)
00253           {
00254             for (n = 0;  n < N;  n++)
00255               {
00256                 i = n + NFirst;
00257                 (*x).elts[n][m] = *(far + iq*nx + i);
00258               }
00259             m++;
00260           }
00261       else
00262         {
00263           nq = ort_list[is][0];
00264           for (iq = 1;  iq <= nq;  iq++)
00265             {
00266               for (n = 0;  n < N;  n++)
00267                 {
00268                   i = n + NFirst;
00269                   (*x).elts[n][m] = *(far + ort_list[is][iq]*nx + i);
00270                 }
00271               m++;
00272             }
00273         }
00274     }
00275 
00276 
00277   /*----- Set up columns of X matrix corresponding to ideal time series -----*/
00278   for (is = 0;  is < num_ideal_files;  is++)
00279     {
00280       far = MRI_FLOAT_PTR (ideal_array[is]);
00281       nx = ideal_array[is]->nx;
00282       ny = ideal_array[is]->ny;
00283 
00284       if (ideal_list[is] == NULL)
00285         for (iq = 0;  iq < ny;  iq++)
00286           {
00287             for (n = 0;  n < N;  n++)
00288               {
00289                 i = n + NFirst;
00290                 (*x).elts[n][m] = *(far + iq*nx + i);
00291               }
00292             
00293             m++;
00294           }
00295       else
00296         {
00297           nq = ideal_list[is][0];
00298           for (iq = 1;  iq <= nq;  iq++)
00299             {
00300               for (n = 0;  n < N;  n++)
00301                 {
00302                   i = n + NFirst;
00303                   (*x).elts[n][m] = *(far + ideal_list[is][iq]*nx + i);
00304                 }
00305 
00306               m++;
00307             }
00308         }
00309     }
00310 
00311 
00312   /*----- Remove row if ideal contains value over 33333 -----*/
00313   Ngood = N;
00314   m = 0;
00315   for (n = 0;  n < N;  n++)
00316     {
00317       for (is = q;  is < p;  is++)
00318         {
00319           if ((*x).elts[n][is] >= BIG_NUMBER)  break;
00320         }
00321       if (is < p)
00322         {
00323           Ngood--;
00324         }
00325       else
00326         {
00327           good_list[m] = n;
00328           m++;
00329         }
00330     }
00331   matrix_extract_rows ((*x), Ngood, good_list, &xgood);
00332   matrix_equate (xgood, x);
00333 
00334 
00335   /*----- Find min, max, and ave for each column of the X matrix -----*/
00336   for (is = 0;  is < p;  is++)
00337     {      
00338       x_bot[is] = x_top[is] = (*x).elts[0][is];
00339       x_ave[is] = 0.0;
00340       for (n = 0;  n < Ngood;  n++)
00341         {
00342           if ((*x).elts[n][is] < x_bot[is])  x_bot[is] = (*x).elts[n][is];  
00343           if ((*x).elts[n][is] > x_top[is])  x_top[is] = (*x).elts[n][is];
00344           x_ave[is] += (*x).elts[n][is] / Ngood;
00345         }
00346     }
00347   
00348   
00349   matrix_destroy (&xgood);
00350 
00351   return (Ngood);
00352 
00353 }

int init_regression_analysis int    q,
int    p,
int    N,
int    num_idealts,
matrix    xdata,
matrix   x_base,
matrix   xtxinvxt_base,
matrix   x_ideal,
matrix   xtxinvxt_ideal,
float **    rarray
 

the new polort way: Legendre P_m(t) - 29 Mar 2005 *

Definition at line 362 of file fim+.c.

References calc_coef(), calc_matrices(), calc_resids(), column_to_vector(), vector::elts, free, malloc, matrix_destroy(), matrix_initialize(), MTEST, p, q, rank_darray(), vector_destroy(), and vector_initialize().

Referenced by calculate_results().

00370                                        :  (1/(X'X))X'  for baseline model */
00371   matrix * x_ideal,           /* extracted X matrices  for ideal models */
00372   matrix * xtxinvxt_ideal,    /* matrix:  (1/(X'X))X'  for ideal models */
00373   float ** rarray             /* ranked arrays of ideal time series */
00374 )
00375 
00376 {
00377   int * plist = NULL;         /* list of model parameters */
00378   int ip, it;                 /* parameter indices */
00379   int is, js;                 /* ideal indices */ 
00380   int jm;                     /* lag index */
00381   int ok;                     /* flag for successful matrix calculation */
00382   matrix xtxinv_temp;         /* intermediate results */
00383   vector ideal;               /* ideal vector */
00384   vector coef_temp;           /* intermediate results */
00385   vector xres;                /* vector of residuals */
00386   float sse_base;             /* baseline model error sum of squares */ 
00387         
00388 
00389   /*----- Initialize matrix -----*/
00390   matrix_initialize (&xtxinv_temp);
00391   vector_initialize (&ideal);
00392   vector_initialize (&coef_temp);
00393   vector_initialize (&xres);
00394 
00395 
00396   /*----- Allocate memory -----*/
00397   plist = (int *) malloc (sizeof(int) * p);   MTEST (plist);
00398 
00399 
00400   /*----- Initialize matrices for the baseline model -----*/
00401   for (ip = 0;  ip < q;  ip++)
00402     plist[ip] = ip;
00403   ok = calc_matrices (xdata, q, plist, x_base, &xtxinv_temp, xtxinvxt_base);
00404   if (!ok)  { matrix_destroy (&xtxinv_temp);  return (0); };
00405 
00406 
00407   /*----- Initialize matrices for ideal functions -----*/
00408   for (is = 0;  is < num_idealts;  is++)
00409     {
00410       for (ip = 0;  ip < q;  ip++)
00411         {
00412           plist[ip] = ip;
00413         }
00414 
00415       plist[q] = q + is;
00416 
00417       ok = calc_matrices (xdata, q+1, plist, 
00418                           &(x_ideal[is]), &xtxinv_temp, &(xtxinvxt_ideal[is]));
00419       if (!ok)  { matrix_destroy (&xtxinv_temp);  return (0); };
00420     }
00421 
00422 
00423   /*----- Set up the ranked array for each ideal -----*/
00424   for (is = 0;  is < num_idealts;  is++)
00425     {
00426       /*----- Convert column of matrix to vector -----*/
00427       column_to_vector (xdata, q+is, &ideal);
00428 
00429       /*----- Calculate regression coefficients for baseline model -----*/
00430       calc_coef (*xtxinvxt_base, ideal, &coef_temp);
00431 
00432       /*----- Calculate the error sum of squares for the baseline model -----*/
00433       sse_base = calc_resids (*x_base, coef_temp, ideal, &xres);
00434     
00435       /*----- Form rank array from residual array -----*/
00436       rarray[is] = rank_darray (N, xres.elts);
00437 
00438     }
00439 
00440 
00441   /*----- Destroy matrix -----*/
00442   matrix_destroy (&xtxinv_temp);
00443   vector_destroy (&ideal);
00444   vector_destroy (&coef_temp);
00445   vector_destroy (&xres);
00446 
00447 
00448   /*----- Deallocate memory -----*/
00449   free (plist);   plist = NULL;
00450 
00451 
00452   return (1);
00453 }

double legendre double    x,
int    m
 

Definition at line 60 of file fim+.c.

Referenced by init_indep_var_matrix(), and legendre().

00061 {
00062    if( m < 0 ) return 1.0 ;    /* bad input */
00063 
00064    switch( m ){                /*** P_m(x) for m=0..20 ***/
00065     case 0: return 1.0 ;
00066     case 1: return x ;
00067     case 2: return (3.0*x*x-1.0)/2.0 ;
00068     case 3: return (5.0*x*x-3.0)*x/2.0 ;
00069     case 4: return ((35.0*x*x-30.0)*x*x+3.0)/8.0 ;
00070     case 5: return ((63.0*x*x-70.0)*x*x+15.0)*x/8.0 ;
00071     case 6: return (((231.0*x*x-315.0)*x*x+105.0)*x*x-5.0)/16.0 ;
00072     case 7: return (((429.0*x*x-693.0)*x*x+315.0)*x*x-35.0)*x/16.0 ;
00073     case 8: return ((((6435.0*x*x-12012.0)*x*x+6930.0)*x*x-1260.0)*x*x+35.0)/128.0;
00074 
00075            /** 07 Feb 2005: this part generated by Maple, then hand massaged **/
00076 
00077     case 9:
00078       return (0.24609375e1 + (-0.3609375e2 + (0.140765625e3 + (-0.20109375e3
00079               + 0.949609375e2 * x * x) * x * x) * x * x) * x * x) * x;
00080 
00081     case 10:
00082       return -0.24609375e0 + (0.1353515625e2 + (-0.1173046875e3 +
00083               (0.3519140625e3 + (-0.42732421875e3 + 0.18042578125e3 * x * x)
00084              * x * x) * x * x) * x * x) * x * x;
00085 
00086     case 11:
00087       return (-0.270703125e1 + (0.5865234375e2 + (-0.3519140625e3 +
00088              (0.8546484375e3 + (-0.90212890625e3 + 0.34444921875e3 * x * x)
00089              * x * x) * x * x) * x * x) * x * x) * x;
00090 
00091     case 12:
00092       return 0.2255859375e0 + (-0.17595703125e2 + (0.2199462890625e3 +
00093              (-0.99708984375e3 + (0.20297900390625e4 + (-0.1894470703125e4
00094              + 0.6601943359375e3 * x * x) * x * x) * x * x) * x * x) * x * x)
00095              * x * x;
00096 
00097     case 13:
00098       return (0.29326171875e1 + (-0.87978515625e2 + (0.7478173828125e3 +
00099              (-0.270638671875e4 + (0.47361767578125e4 + (-0.3961166015625e4
00100              + 0.12696044921875e4 * x * x) * x * x) * x * x) * x * x) * x * x)
00101             * x * x) * x;
00102 
00103     case 14:
00104       return -0.20947265625e0 + (0.2199462890625e2 + (-0.37390869140625e3 +
00105              (0.236808837890625e4 + (-0.710426513671875e4 +
00106              (0.1089320654296875e5 + (-0.825242919921875e4 +
00107             0.244852294921875e4 * x * x) * x * x) * x * x) * x * x) * x * x)
00108            * x * x) * x * x;
00109 
00110     case 15:
00111       return (-0.314208984375e1 + (0.12463623046875e3 + (-0.142085302734375e4
00112             + (0.710426513671875e4 + (-0.1815534423828125e5 +
00113               (0.2475728759765625e5 + (-0.1713966064453125e5 +
00114                0.473381103515625e4 * x * x) * x * x) * x * x) * x * x)
00115              * x * x) * x * x) * x * x) * x;
00116 
00117     case 16:
00118       return 0.196380615234375e0 + (-0.26707763671875e2 + (0.5920220947265625e3
00119             + (-0.4972985595703125e4 + (0.2042476226806641e5 +
00120               (-0.4538836059570312e5 + (0.5570389709472656e5 +
00121                (-0.3550358276367188e5 + 0.9171758880615234e4 * x * x) * x * x)
00122             * x * x) * x * x) * x * x) * x * x) * x * x) * x * x;
00123 
00124     case 17:
00125       return (0.3338470458984375e1 + (-0.169149169921875e3 +
00126              (0.2486492797851562e4 + (-0.1633980981445312e5 +
00127              (0.5673545074462891e5 + (-0.1114077941894531e6 +
00128              (0.1242625396728516e6 + (-0.7337407104492188e5 +
00129               0.1780400253295898e5 * x * x) * x * x) * x * x) * x * x)
00130            * x * x) * x * x) * x * x) * x * x) * x;
00131 
00132     case 18:
00133       return -0.1854705810546875e0 + (0.3171546936035156e2 +
00134             (-0.8880331420898438e3 + (0.9531555725097656e4 +
00135             (-0.5106190567016602e5 + (0.153185717010498e6 +
00136             (-0.2692355026245117e6 + (0.275152766418457e6 +
00137             (-0.1513340215301514e6 + 0.3461889381408691e5 * x * x) * x * x)
00138            * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x * x;
00139 
00140     case 19:
00141       return (-0.3523941040039062e1 + (0.2220082855224609e3 +
00142              (-0.4084952453613281e4 + (0.3404127044677734e5 +
00143              (-0.153185717010498e6 + (0.4038532539367676e6 +
00144              (-0.6420231216430664e6 + (0.6053360861206055e6 +
00145              (-0.3115700443267822e6 + 0.6741574058532715e5 * x * x) * x * x)
00146           * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x;
00147 
00148     case 20:
00149       return 0.1761970520019531e0 + (-0.3700138092041016e2 +
00150             (0.127654764175415e4 + (-0.1702063522338867e5 +
00151             (0.1148892877578735e6 + (-0.4442385793304443e6 +
00152             (0.1043287572669983e7 + (-0.1513340215301514e7 +
00153             (0.1324172688388824e7 + (-0.6404495355606079e6 +
00154              0.1314606941413879e6 * x * x) * x * x) * x * x) * x * x) * x * x)
00155             * x * x) * x * x) * x * x) * x * x) * x * x;
00156    }
00157 
00158 #if 0
00159    /* order out of range: return Chebyshev instead (it's easy) */
00160 
00161         if(  x >=  1.0 ) x = 0.0 ;
00162    else if ( x <= -1.0 ) x = 3.14159265358979323846 ;
00163    else                  x = acos(x) ;
00164    return cos(m*x) ;
00165 #else
00166    /** if here, m > 20 ==> use recurrence relation **/
00167 
00168    { double pk, pkm1, pkm2 ; int k ;
00169      pkm2 = legendre( x , 19 ) ;
00170      pkm1 = legendre( x , 20 ) ;
00171      for( k=21 ; k <= m ; k++ , pkm2=pkm1 , pkm1=pk )
00172        pk = ((2.0*k-1.0)*x*pkm1 - (k-1.0)*pkm2)/k ;
00173      return pk ;
00174    }
00175 #endif
00176 }

float percent_change float    num,
float    den
 

Definition at line 599 of file fim+.c.

References sign.

Referenced by regression_analysis().

00600 {
00601   const float EPSILON = 1.0e-10;      /* guard against divide by zero */
00602   const float MAX_PERCENT = 1000.0;   /* limit maximum percent change */ 
00603   float PrcntChange;
00604 
00605 
00606   if (fabs(den) < EPSILON)
00607     PrcntChange = sign(num) * MAX_PERCENT;
00608   else
00609     PrcntChange  = 100.0 * num / den;
00610  
00611   if (PrcntChange >  MAX_PERCENT)   PrcntChange =  MAX_PERCENT;
00612   if (PrcntChange < -MAX_PERCENT)   PrcntChange = -MAX_PERCENT;
00613 
00614   return (PrcntChange);
00615 }

void regression_analysis int    N,
int    q,
int    num_idealts,
matrix    x_base,
matrix    xtxinvxt_base,
matrix   x_ideal,
matrix   xtxinvxt_ideal,
vector    y,
float *    x_bot,
float *    x_ave,
float *    x_top,
float **    rarray,
int *    output_type,
float *    FimParams
 

Definition at line 624 of file fim+.c.

References calc_coef(), calc_QuadrantCC(), calc_resids(), calc_rsqr(), calc_SpearmanCC(), calc_sse(), vector::elts, FIM_Average, FIM_Baseline, FIM_BestIndex, FIM_Correlation, FIM_FitCoef, FIM_PrcntChange, FIM_PrcntFromAve, FIM_PrcntFromTop, FIM_QuadrantCC, FIM_SigmaResid, FIM_SpearmanCC, FIM_Topline, free, percent_change(), q, rank_darray(), vector_destroy(), vector_equate(), and vector_initialize().

Referenced by calculate_results().

00629                                      :  (1/(X'X))X'  for baseline model */
00630   matrix * x_ideal,         /* extracted X matrices  for ideal models */
00631   matrix * xtxinvxt_ideal,  /* matrix:  (1/(X'X))X'  for ideal models */
00632   vector y,                 /* vector of measured data */
00633   float * x_bot,            /* minimum of stimulus time series */
00634   float * x_ave,            /* average of stimulus time series */
00635   float * x_top,            /* maximum of stimulus time series */
00636   float ** rarray,          /* ranked arrays of ideal time series */
00637   int * output_type,        /* list of operator requested outputs */
00638   float * FimParams         /* output fim parameters */
00639 )
00640 
00641 {
00642   const float EPSILON = 1.0e-05;   /* protection against divide by zero */
00643   int is;                   /* input ideal index */
00644   float sse_base;           /* error sum of squares, baseline model */
00645   float sse_ideal;          /* error sum of squares, ideal model */
00646   vector coef_temp;         /* intermediate results */
00647   vector coef_best;         /* best results */
00648   vector yres;              /* vector of residuals */
00649   float rtemp, rbest;       /* best correlation coefficient */  
00650   float stemp, sbest;       /* best Spearman correlation coefficient */
00651   float qtemp, qbest;       /* best quadrant correlation coefficient */
00652   float mse;                /* mean square error (sample variance) */
00653 
00654   float * sarray = NULL;    /* ranked array of measurement data */
00655 
00656                             /* fim output parameters */
00657   float FitCoef      = 0.0; 
00658   int   BestIndex    = 0;
00659   float PrcntChange  = 0.0; 
00660   float Baseline     = 0.0; 
00661   float Correlation  = 0.0; 
00662   float PrcntFromAve = 0.0;
00663   float Average      = 0.0; 
00664   float PrcntFromTop = 0.0; 
00665   float Topline      = 0.0; 
00666   float SigmaResid   = 0.0; 
00667   float SpearmanCC   = 0.0; 
00668   float QuadrantCC   = 0.0;
00669 
00670 
00671   /*----- Initialization -----*/
00672   vector_initialize (&coef_temp);
00673   vector_initialize (&coef_best);
00674   vector_initialize (&yres);
00675 
00676 
00677   /*----- Calculate regression coefficients for baseline model -----*/
00678   calc_coef (xtxinvxt_base, y, &coef_temp);
00679 
00680 
00681   /*----- Calculate the error sum of squares for the baseline model -----*/ 
00682   sse_base = calc_resids (x_base, coef_temp, y, &yres);
00683 
00684     
00685   /*----- Form rank array from y array -----*/
00686   if (output_type[FIM_SpearmanCC] || output_type[FIM_QuadrantCC]) 
00687     sarray = rank_darray (N, yres.elts);
00688 
00689 
00690   /*----- Determine the best ideal reference for this voxel -----*/
00691   rbest = 0.0;  sbest = 0.0;  qbest = 0.0;  BestIndex = -1;
00692   for (is = 0;  is < num_idealts;  is++)
00693     {
00694 
00695       /*----- Calculate regression coefficients for ideal model -----*/
00696       calc_coef (xtxinvxt_ideal[is], y, &coef_temp);
00697 
00698 
00699       /*----- Calculate the error sum of squares for the ideal model -----*/ 
00700       sse_ideal = calc_sse (x_ideal[is], coef_temp, y);
00701 
00702 
00703       /*----- Calculate partial R^2 for this ideal -----*/
00704       rtemp = calc_rsqr (sse_ideal, sse_base);
00705 
00706 
00707       if (rtemp >= rbest)
00708         {
00709           rbest = rtemp;
00710           BestIndex = is;
00711           vector_equate (coef_temp, &coef_best);
00712           if (num_idealts == 1)
00713             mse = sse_ideal / (N-q-1);
00714           else
00715             mse = sse_ideal / (N-q-2);
00716         }
00717 
00718 
00719       /*----- Calculate the Spearman rank correlation coefficient -----*/
00720       if (output_type[FIM_SpearmanCC])
00721         { 
00722           stemp = calc_SpearmanCC (N, rarray[is], sarray);
00723           if (fabs(stemp) > fabs(sbest))  sbest = stemp;
00724         }
00725 
00726 
00727       /*----- Calculate the Quadrant correlation coefficient -----*/
00728       if (output_type[FIM_QuadrantCC])
00729         { 
00730           qtemp = calc_QuadrantCC (N, rarray[is], sarray);
00731           if (fabs(qtemp) > fabs(qbest))  qbest = qtemp;
00732         }
00733     }
00734   
00735   if ((0 <= BestIndex) && (BestIndex < num_idealts))
00736     {
00737       float Top, Ave, Base, Center;
00738       int ip;
00739       
00740       Top  = x_top[q+BestIndex];  
00741       Ave  = x_ave[q+BestIndex];  
00742       Base = x_bot[q+BestIndex];
00743       
00744       
00745       FitCoef = coef_best.elts[q];
00746       Correlation = sqrt(rbest);
00747       if (FitCoef < 0.0)  Correlation = -Correlation;
00748       
00749       Center = 0.0;
00750       for (ip = 0;  ip < q;  ip++)
00751         Center += coef_best.elts[ip] * x_ave[ip];
00752       
00753       Baseline = Center + FitCoef*Base;
00754       Average  = Center + FitCoef*Ave;
00755       Topline  = Center + FitCoef*Top;
00756 
00757       
00758       
00759       PrcntChange  = percent_change (FitCoef * (Top-Base), Baseline);
00760       PrcntFromAve = percent_change (FitCoef * (Top-Base), Average);
00761       PrcntFromTop = percent_change (FitCoef * (Top-Base), Topline);
00762       
00763       SigmaResid = sqrt(mse);
00764       
00765       SpearmanCC = sbest;
00766       
00767       QuadrantCC = qbest;
00768     }
00769   
00770 
00771   /*----- Save output parameters -----*/
00772   FimParams[FIM_FitCoef]      = FitCoef;
00773   FimParams[FIM_BestIndex]    = BestIndex + 1.0;
00774   FimParams[FIM_PrcntChange]  = PrcntChange;
00775   FimParams[FIM_Baseline]     = Baseline;
00776   FimParams[FIM_Correlation]  = Correlation;
00777   FimParams[FIM_PrcntFromAve] = PrcntFromAve;
00778   FimParams[FIM_Average]      = Average;
00779   FimParams[FIM_PrcntFromTop] = PrcntFromTop;
00780   FimParams[FIM_Topline]      = Topline;
00781   FimParams[FIM_SigmaResid]   = SigmaResid;
00782   FimParams[FIM_SpearmanCC]   = SpearmanCC;
00783   FimParams[FIM_QuadrantCC]   = QuadrantCC;
00784 
00785 
00786   /*----- Dispose of vectors -----*/
00787   vector_destroy (&coef_temp);
00788   vector_destroy (&coef_best);
00789   vector_destroy (&yres);
00790 
00791 
00792   /*----- Deallocate memory -----*/
00793   if (sarray != NULL)
00794     { free (sarray);  sarray = NULL; }
00795 
00796 }

void report_results int *    output_type,
float *    FimParams,
char **    label
 

Definition at line 806 of file fim+.c.

References lbuf, MAX_OUTPUT_TYPE, OUTPUT_TYPE_name, and sbuf.

00812 {
00813   int ip;                     /* parameter index */
00814 
00815   
00816   if( label != NULL ){  /* assemble this 1 line at a time from sbuf */
00817     
00818     lbuf[0] = '\0' ;   /* make this a 0 length string to start */
00819     
00820     /** for each reference, make a string into sbuf **/
00821 
00822     for (ip = 0;  ip < MAX_OUTPUT_TYPE;  ip++)
00823       if (output_type[ip])
00824           {
00825             sprintf (sbuf, "%12s = %10.4f \n", 
00826                      OUTPUT_TYPE_name[ip], FimParams[ip]);
00827             strcat (lbuf, sbuf);
00828           }    
00829 
00830     
00831     *label = lbuf ;  /* send address of lbuf back in what label points to */
00832   }
00833   
00834 }

float sign float    x
 

Definition at line 543 of file fim+.c.

00545 {
00546   if (x > 0.0)  return (1.0);
00547   if (x < 0.0)  return (-1.0);
00548   return (0.0);
00549 }

Variable Documentation

char lbuf[4096] [static]
 

Definition at line 801 of file fim+.c.

Referenced by report_results().

char* OUTPUT_TYPE_name[MAX_OUTPUT_TYPE] [static]
 

Initial value:

  { "Fit Coef", "Best Index", "% Change", "% From Ave", "Baseline", "Average", 
    "Correlation", "% From Top", "Topline", "Sigma Resid",
    "Spearman CC", "Quadrant CC" }

Definition at line 36 of file fim+.c.

Referenced by report_results().

char sbuf[256] [static]
 

Definition at line 802 of file fim+.c.

Referenced by report_results().

 

Powered by Plone

This site conforms to the following standards: