Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
delay.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #include "RegAna.c"
00017 #include "ranks.c"
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 int init_indep_var_matrix
00028 (
00029 int q,
00030 int p,
00031
00032 int NFirst,
00033 int N,
00034 int polort,
00035 int num_ort_files,
00036 int num_ideal_files,
00037 MRI_IMAGE ** ort_array,
00038 int ** ort_list,
00039 MRI_IMAGE ** ideal_array,
00040 int ** ideal_list,
00041 float * x_bot,
00042 float * x_ave,
00043 float * x_top,
00044 int * good_list,
00045 matrix * x
00046 )
00047
00048 {
00049 const int BIG_NUMBER = 33333;
00050 int i;
00051 int m;
00052 int n;
00053 int is;
00054 float * far = NULL;
00055 int nx, ny, iq, nq;
00056 int Ngood;
00057 matrix xgood;
00058
00059
00060
00061 matrix_create (N, p, x);
00062 matrix_initialize (&xgood);
00063
00064
00065
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
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
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
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
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
00189
00190
00191 int init_delay
00192 (
00193 int q,
00194 int p,
00195
00196 int N,
00197 int num_idealts,
00198 matrix xdata,
00199 matrix * x_base,
00200 matrix * xtxinvxt_base,
00201 matrix * x_ideal,
00202 matrix * xtxinvxt_ideal,
00203 float ** rarray
00204 )
00205
00206 {
00207 int * plist = NULL;
00208 int ip, it;
00209 int is, js;
00210 int jm;
00211 int ok;
00212 matrix xtxinv_temp;
00213 vector ideal;
00214 vector coef_temp;
00215 vector xres;
00216 float sse_base;
00217
00218
00219
00220 matrix_initialize (&xtxinv_temp);
00221 vector_initialize (&ideal);
00222 vector_initialize (&coef_temp);
00223 vector_initialize (&xres);
00224
00225
00226
00227 plist = (int *) malloc (sizeof(int) * p); MTEST (plist);
00228
00229
00230
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
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
00254 for (is = 0; is < num_idealts; is++)
00255 {
00256
00257 column_to_vector (xdata, q+is, &ideal);
00258
00259
00260 calc_coef (*xtxinvxt_base, ideal, &coef_temp);
00261
00262
00263 sse_base = calc_resids (*x_base, coef_temp, ideal, &xres);
00264
00265
00266 rarray[is] = rank_darray (N, xres.elts);
00267
00268 }
00269
00270
00271
00272 matrix_destroy (&xtxinv_temp);
00273 vector_destroy (&ideal);
00274 vector_destroy (&coef_temp);
00275 vector_destroy (&xres);
00276
00277
00278
00279 free (plist); plist = NULL;
00280
00281
00282 return (1);
00283 }
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
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