/*---------------------------------------------------------------------------*/ /* This file contains routines for implementing the 3ddelay functions. Based on fim+.c by B. Douglas Ward This software is copyrighted and owned by the Medical College of Wisconsin. See the file README.Copyright for details. */ /*---------------------------------------------------------------------------*/ /* Include software for linear regression analysis and sorting numbers. */ #include "RegAna.c" #include "ranks.c" /*---------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------*/ /* Initialize independent variable X matrix */ int init_indep_var_matrix ( int q, /* number of parameters in the baseline model */ int p, /* number of parameters in the baseline model plus number of ideals */ int NFirst, /* index of first image to use in the analysis */ int N, /* total number of images used in the analysis */ int polort, /* degree of polynomial ort */ int num_ort_files, /* number of ort time series files */ int num_ideal_files, /* number of ideal time series files */ MRI_IMAGE ** ort_array, /* ort time series arrays */ int ** ort_list, /* list of ort time series */ MRI_IMAGE ** ideal_array, /* ideal time series arrays */ int ** ideal_list, /* list of ideal time series */ float * x_bot, /* minimum of stimulus time series */ float * x_ave, /* average of stimulus time series */ float * x_top, /* maximum of stimulus time series */ int * good_list, /* list of good time points */ matrix * x /* independent variable matrix */ ) { const int BIG_NUMBER = 33333; int i=-1; /* time index */ int m=-1; /* X matrix column index */ int n=-1; /* X matrix row index */ int is=-1; /* input ideal index */ float * far = NULL; int nx=-1, ny=-1, iq=-1, nq=-1; int Ngood=-1; int maxcols=-1; matrix xgood; if (num_ort_files) { fprintf(stderr,"Not ready to deal with ort files at all!\n"); exit(1); } /*----- Initialize X matrix -----*/ maxcols = 0; for (is = 0; is < num_ort_files; is++) {/* add orts */ if (ort_list[is] == NULL) maxcols += ort_array[is]->ny; else maxcols += ort_list[is][0]+1; } for (is = 0; is < num_ideal_files; is++) {/* add ideals */ if (ideal_list[is] == NULL) maxcols += ideal_array[is]->ny; else maxcols += ideal_list[is][0]+1; } maxcols = maxcols + polort + 1; /* add polorts */ /* fprintf(stderr,"N=%d,p=%d, polort=%d, maxcols = %d\n", N, p, polort, maxcols); */ matrix_create (N+NFirst, maxcols, x); matrix_initialize (&xgood); /*----- Set up columns of X matrix corresponding to polynomial orts -----*/ for (m = 0; m <= polort; m++) for (n = 0; n < N; n++) { i = n + NFirst; (*x).elts[n][m] = pow ((double)i, (double)m); } /*----- Set up columns of X matrix corresponding to ort time series -----*/ for (is = 0; is < num_ort_files; is++) { far = MRI_FLOAT_PTR (ort_array[is]); nx = ort_array[is]->nx; ny = ort_array[is]->ny; if (ort_list[is] == NULL) for (iq = 0; iq < ny; iq++) { for (n = 0; n < N; n++) { i = n + NFirst; (*x).elts[n][m] = *(far + iq*nx + i); } m++; } else { nq = ort_list[is][0]; for (iq = 1; iq <= nq; iq++) { for (n = 0; n < N; n++) { i = n + NFirst; (*x).elts[n][m] = *(far + ort_list[is][iq]*nx + i); } m++; } } } /*----- Set up columns of X matrix corresponding to ideal time series -----*/ for (is = 0; is < num_ideal_files; is++) { far = MRI_FLOAT_PTR (ideal_array[is]); nx = ideal_array[is]->nx; ny = ideal_array[is]->ny; if (ideal_list[is] == NULL) for (iq = 0; iq < ny; iq++) { for (n = 0; n < N; n++) { i = n + NFirst; (*x).elts[n][m] = *(far + iq*nx + i); } m++; } else { nq = ideal_list[is][0]; for (iq = 1; iq <= nq; iq++) { for (n = 0; n < N; n++) { i = n + NFirst; (*x).elts[n][m] = *(far + ideal_list[is][iq]*nx + i); } m++; } } } /*----- Remove row if ideal contains value over 33333 -----*/ Ngood = N; m = 0; for (n = 0; n < N; n++) { for (is = q; is < p; is++) { if ((*x).elts[n][is] >= BIG_NUMBER) break; } if (is < p) { Ngood--; } else { good_list[m] = n; m++; } } matrix_extract_rows ((*x), Ngood, good_list, &xgood); matrix_equate (xgood, x); /*----- Find min, max, and ave for each column of the X matrix -----*/ for (is = 0; is < p; is++) { x_bot[is] = x_top[is] = (*x).elts[0][is]; x_ave[is] = 0.0; for (n = 0; n < Ngood; n++) { if ((*x).elts[n][is] < x_bot[is]) x_bot[is] = (*x).elts[n][is]; if ((*x).elts[n][is] > x_top[is]) x_top[is] = (*x).elts[n][is]; x_ave[is] += (*x).elts[n][is] / Ngood; } } matrix_destroy (&xgood); return (Ngood); } /*---------------------------------------------------------------------------*/ /* Initialization for the delay analysis. */ int init_delay ( int q, /* number of parameters in the baseline model */ int p, /* number of parameters in the baseline model plus number of ideals */ int N, /* total number of images used in the analysis */ int num_idealts, /* number of ideal time series */ matrix xdata, /* independent variable matrix */ matrix * x_base, /* extracted X matrix for baseline model */ matrix * xtxinvxt_base, /* matrix: (1/(X'X))X' for baseline model */ matrix * x_ideal, /* extracted X matrices for ideal models */ matrix * xtxinvxt_ideal, /* matrix: (1/(X'X))X' for ideal models */ float ** rarray /* ranked arrays of ideal time series */ ) { int * plist = NULL; /* list of model parameters */ int ip, it; /* parameter indices */ int is, js; /* ideal indices */ int jm; /* lag index */ int ok; /* flag for successful matrix calculation */ matrix xtxinv_temp; /* intermediate results */ vector ideal; /* ideal vector */ vector coef_temp; /* intermediate results */ vector xres; /* vector of residuals */ float sse_base; /* baseline model error sum of squares */ /*----- Initialize matrix -----*/ matrix_initialize (&xtxinv_temp); vector_initialize (&ideal); vector_initialize (&coef_temp); vector_initialize (&xres); /*----- Allocate memory -----*/ plist = (int *) malloc (sizeof(int) * p); MTEST (plist); /*----- Initialize matrices for the baseline model -----*/ for (ip = 0; ip < q; ip++) plist[ip] = ip; ok = calc_matrices (xdata, q, plist, x_base, &xtxinv_temp, xtxinvxt_base); if (!ok) { matrix_destroy (&xtxinv_temp); return (0); }; /*----- Initialize matrices for ideal functions -----*/ for (is = 0; is < num_idealts; is++) { for (ip = 0; ip < q; ip++) { plist[ip] = ip; } plist[q] = q + is; ok = calc_matrices (xdata, q+1, plist, &(x_ideal[is]), &xtxinv_temp, &(xtxinvxt_ideal[is])); if (!ok) { matrix_destroy (&xtxinv_temp); return (0); }; } /*----- Set up the ranked array for each ideal -----*/ for (is = 0; is < num_idealts; is++) { /*----- Convert column of matrix to vector -----*/ column_to_vector (xdata, q+is, &ideal); /*----- Calculate regression coefficients for baseline model -----*/ calc_coef (*xtxinvxt_base, ideal, &coef_temp); /*----- Calculate the error sum of squares for the baseline model -----*/ sse_base = calc_resids (*x_base, coef_temp, ideal, &xres); /*----- Form rank array from residual array -----*/ rarray[is] = rank_darray (N, xres.elts); } /*----- Destroy matrix -----*/ matrix_destroy (&xtxinv_temp); vector_destroy (&ideal); vector_destroy (&coef_temp); vector_destroy (&xres); /*----- Deallocate memory -----*/ free (plist); plist = NULL; return (1); } /*---------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------*/ /* Calculate the sign function. */ float sign (float x) { if (x > 0.0) return (1.0); if (x < 0.0) return (-1.0); return (0.0); } /*---------------------------------------------------------------------------*/ int Read_part_file_delay(float *x, char *f_name, int a, int b) { float *far=NULL; int i=0, cnt=0; MRI_IMAGE *im_data=NULL; im_data = mri_read_1D(f_name); if (!im_data) { fprintf(stderr,"Failed to read 1D file."); return(-1);} far = MRI_FLOAT_PTR(im_data); if (im_data->ny != 1) { fprintf(stderr,"More than one column in 1D file."); return(-1);} if (b>=im_data->nx) { fprintf(stderr,"Nlast >= number of values in 1D file!"); return(-1);} cnt = 0; for (i=a;i<=b;++i) { x[cnt] = far[i]; #ifdef ZDBG if (cnt < 10) { fprintf(stderr,"%d: %.3f\t", i, x[cnt]);} #endif ++cnt; } mri_free(im_data); im_data = NULL; return(cnt); } int Read_part_file_delay_OLD (float *x, char *f_name, int a, int b) { int cnt=0,ex,line_num; float buf; FILE*file_in; file_in = fopen (f_name,"r"); if (file_in == NULL) { printf ("\aCould not open %s \n",f_name); printf ("Exiting @ Read_file function\n"); exit (0); } if (a > b || a==0) { printf ("\a\n\33[1mError in (from , to) line numbers\n\33\[0m"); printf ("Exiting @Read_part_file function \n"); exit (0); } line_num = 1; if (a == 1) { ex = fscanf (file_in,"%f",&x[cnt]); ++cnt; } else ex = fscanf (file_in,"%f",&buf); ++line_num; while (ex != EOF && line_num <= b) { if (line_num >=a && line_num <=b) { ex = fscanf (file_in,"%f",&x[cnt]); ++cnt; if (ex == EOF) --cnt; } else { ex = fscanf (file_in,"%f",&buf); } ++line_num; } if (ex == EOF) { --line_num; printf ("\n\33[1mEOF reached before line \33[0m%d\n",b); printf ("Only %d lines were read, from line %d to %d\n",cnt,a,line_num-1); } fclose (file_in); return (cnt); } /*---------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------*/