1 #include "mrilib.h"
2 #include "remla.c"
3
main(int argc,char * argv[])4 int main( int argc , char *argv[] )
5 {
6 MRI_IMAGE *inim , *outim ; float *iv ;
7 int iarg , ii,jj , nreg , ntime , *tau=NULL , rnum ;
8 NI_element *nelmat=NULL ; char *matname=NULL ;
9 MTYPE rhomax=0.7 , bmax=0.7 ; int rhonum=7 , bnum=14 ;
10 char *cgl , *rst ;
11 matrix X ; vector y ;
12 float cput ;
13 reml_collection *rrcol ;
14
15 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
16 printf(
17 "Usage: 1dREMLfit [option] file.1D\n"
18 "Least squares fit with REML estimation of the ARMA(1,1) noise.\n"
19 "\n"
20 "Options (the first one is mandatory)\n"
21 "------------------------------------\n"
22 " -matrix mmm = Read the matrix 'mmm', which should have been\n"
23 " output from 3dDeconvolve via the '-x1D' option.\n"
24 " -MAXrho rm = Set the max allowed rho parameter to 'rm' (default=0.7).\n"
25 " -Nrho nr = Use 'nr' values for the rho parameter (default=7).\n"
26 " -MAXb bm = Set max allow MA b parameter to 'bm' (default=0.7).\n"
27 " -Nb nb = Use 'nb' values for the b parameter (default=7).\n"
28 ) ;
29 PRINT_COMPILE_DATE ; exit(0) ;
30 }
31
32 iarg = 1 ;
33 while( iarg < argc && argv[iarg][0] == '-' ){
34
35 /** rho and del params **/
36
37 if( strcmp(argv[iarg],"-MAXrho") == 0 ){
38 rhomax = (MTYPE)strtod(argv[++iarg],NULL) ;
39 if( rhomax < 0.3 ) rhomax = 0.3 ;
40 else if( rhomax > 0.9 ) rhomax = 0.9 ;
41 iarg++ ; continue ;
42 }
43 if( strcmp(argv[iarg],"-Nrho") == 0 ){
44 rhonum = (int)strtod(argv[++iarg],NULL) ;
45 if( rhonum < 2 ) rhonum = 2 ;
46 else if( rhonum > 20 ) rhonum = 20 ;
47 iarg++ ; continue ;
48 }
49 if( strcmp(argv[iarg],"-MAXb") == 0 ){
50 bmax = (MTYPE)strtod(argv[++iarg],NULL) ;
51 if( bmax < 0.3 ) bmax = 0.3 ;
52 else if( bmax > 0.9 ) bmax = 0.9 ;
53 iarg++ ; continue ;
54 }
55 if( strcmp(argv[iarg],"-Nb") == 0 ){
56 bnum = (int)strtod(argv[++iarg],NULL) ;
57 if( bnum < 2 ) bnum = 2 ;
58 else if( bnum > 20 ) bnum = 20 ;
59 iarg++ ; continue ;
60 }
61
62 /** -matrix **/
63
64 if( strcmp(argv[iarg],"-matrix") == 0 ){
65 if( nelmat != NULL ) ERROR_exit("More than 1 -matrix option!");
66 nelmat = NI_read_element_fromfile( argv[++iarg] ) ; /* read NIML file */
67 matname = argv[iarg];
68 if( nelmat == NULL ){ /* try to read as a 1D file */
69 MRI_IMAGE *nim ; float *nar ;
70 nim = mri_read_1D(argv[iarg]) ;
71 if( nim != NULL ){ /* construct a minimal NIML element */
72 nelmat = NI_new_data_element( "matrix" , nim->nx ) ;
73 nar = MRI_FLOAT_PTR(nim) ;
74 for( jj=0 ; jj < nim->ny ; jj++ )
75 NI_add_column( nelmat , NI_FLOAT , nar + nim->nx*jj ) ;
76 mri_free(nim) ;
77 }
78 }
79 if( nelmat == NULL || nelmat->type != NI_ELEMENT_TYPE )
80 ERROR_exit("Can't process -matrix file!");
81 iarg++ ; continue ;
82 }
83
84 ERROR_exit("Unknown option '%s'",argv[iarg]) ;
85 }
86
87 if( iarg >= argc ) ERROR_exit("No 1D file on command line?!") ;
88
89 inim = mri_read_1D( argv[iarg] ) ;
90 if( inim == NULL ) ERROR_exit("Can't read 1D file %s",argv[iarg]) ;
91
92 nreg = nelmat->vec_num ;
93 ntime = nelmat->vec_len ;
94 if( ntime != inim->nx )
95 ERROR_exit("matrix vectors are %d long but input 1D file is %d long",
96 ntime,inim->nx) ;
97
98 cgl = NI_get_attribute( nelmat , "GoodList" ) ;
99 if( cgl != NULL ){
100 int Ngoodlist,*goodlist , Nruns,*runs ;
101 NI_int_array *giar ;
102 giar = NI_decode_int_list( cgl , ";," ) ;
103 if( giar == NULL || giar->num < ntime )
104 ERROR_exit("-matrix 'GoodList' badly formatted?") ;
105 Ngoodlist = giar->num ; goodlist = giar->ar ;
106 rst = NI_get_attribute( nelmat , "RunStart" ) ;
107 if( rst != NULL ){
108 NI_int_array *riar = NI_decode_int_list( rst , ";,") ;
109 if( riar == NULL ) ERROR_exit("-matrix 'RunStart' badly formatted?") ;
110 Nruns = riar->num ; runs = riar->ar ;
111 } else {
112 Nruns = 1 ; runs = calloc(sizeof(int),1) ;
113 }
114 rnum = 0 ; tau = (int *)malloc(sizeof(int)*ntime) ;
115 for( ii=0 ; ii < ntime ; ii++ ){
116 jj = goodlist[ii] ;
117 for( ; rnum+1 < Nruns && jj >= runs[rnum+1] ; rnum++ ) ; /*nada*/
118 tau[ii] = jj + 10000*rnum ;
119 }
120 }
121
122 matrix_initialize( &X ) ;
123 matrix_create( ntime , nreg , &X ) ;
124 if( nelmat->vec_typ[0] == NI_FLOAT ){
125 float *cd ;
126 for( jj=0 ; jj < nreg ; jj++ ){
127 cd = (float *)nelmat->vec[jj] ;
128 for( ii=0 ; ii < ntime ; ii++ ) X.elts[ii][jj] = (MTYPE)cd[ii] ;
129 }
130 } else if( nelmat->vec_typ[0] == NI_DOUBLE ){
131 double *cd ;
132 for( jj=0 ; jj < nreg ; jj++ ){
133 cd = (double *)nelmat->vec[jj] ;
134 for( ii=0 ; ii < ntime ; ii++ ) X.elts[ii][jj] = (MTYPE)cd[ii] ;
135 }
136 } else {
137 ERROR_exit("-matrix file stored will illegal data type!?") ;
138 }
139
140 cput = COX_cpu_time() ;
141 rrcol = REML_setup( &X , tau , rhonum,rhomax,bnum,bmax ) ;
142 if( rrcol == NULL ) ERROR_exit("REML setup fails?" ) ;
143 cput = COX_cpu_time() - cput ;
144 INFO_message("REML setup: rows=%d cols=%d %d cases CPU=%.2f",
145 ntime,nreg,rrcol->nset,cput) ;
146
147 cput = COX_cpu_time() ;
148 vector_initialize( &y ) ; vector_create_noinit( ntime , &y ) ;
149 for( jj=0 ; jj < inim->ny ; jj++ ){
150 iv = MRI_FLOAT_PTR(inim) + ntime*jj ;
151 for( ii=0 ; ii < ntime ; ii++ ) y.elts[ii] = (MTYPE)iv[ii] ;
152 (void)REML_find_best_case( &y , rrcol ) ;
153 INFO_message(
154 "Vector #%d: best_rho=%.2f best_b=%.2f best_lam=%.2f best_ssq=%g olsq_ssq=%g",
155 jj, REML_best_rho, REML_best_bb, REML_best_lam, REML_best_ssq, REML_olsq_ssq ) ;
156 }
157 cput = COX_cpu_time() - cput ;
158 INFO_message("REML fitting: CPU=%.2f",cput) ;
159 exit(0) ;
160 }
161