1 #include "mrilib.h"
2 #include "zgaussian.c"
3 #include <time.h>
4 #include <sys/types.h>
5 #include <unistd.h>
6
7 #undef MTYPE
8 #undef MPAIR
9 #define MTYPE double
10 #define MPAIR double_pair
11
12 #undef BIGVAL
13 #define BIGVAL 1.e+38
14
15 #undef LAMBDA
16 #define LAMBDA(a,b) ((b+a)*(1.0+a*b)/(1.0+2.0*a*b+b*b))
17
18 static MTYPE corcut = 0.0001 ;
19
20 #undef TAU
21 #define TAU(i) ((tau==NULL) ? (i) : tau[i])
22
23 #define ALLOW_ARMA51
24
25 /*--------------------------------------------------------------------------*/
26 /*! Setup sparse banded correlation matrix (as an rcmat struct):
27 [ 1 lam lam*rho lam*rho^2 lam*rho^3 ... ]
28 which is the ARMA(1,1) model with the AR parameter a = rho,
29 and the MA parameter b such that (b+a)*(1+a*b)/(1+2*a*b+b*b) = lam.
30 * For reasonable models of FMRI noise, 0 < lam < rho < 0.9.
31 * The maximum bandwidth of the matrix is chosen so that the last
32 correlation element is about 0.01.
33 * tau[i] is the 'true' time index of the i-th data point. This
34 lets you allow for censoring and for inter-run gaps.
35 * If tau==NULL, tau[i] is taken to be i -- that is, no censoring/gaps.
36 *//*------------------------------------------------------------------------*/
37
rcmat_arma11(int nt,int * tau,MTYPE rho,MTYPE lam)38 rcmat * rcmat_arma11( int nt, int *tau, MTYPE rho, MTYPE lam )
39 {
40 rcmat *rcm ;
41 LENTYP *len ;
42 MTYPE **rc , *rii , alam ;
43 int ii , jj , bmax , jbot , itt,jtt ;
44
45 if( nt < 2 ) return NULL ;
46
47 rcm = rcmat_init( nt ) ; /* create sparse matrix struct */
48 len = rcm->len ;
49 rc = rcm->rc ;
50
51 if( rho > 0.99 ) rho = 0.99 ; /* max allowed NN correlation */
52 else if( rho < -0.99 ) rho = -0.99 ;
53
54 /* set maximum bandwidth */
55
56 alam = fabs(lam) ;
57 if( alam >= corcut ){
58 if( rho != 0.0 ) /* bmax is such that last element is about 'corcut' */
59 bmax = 1 + (int)ceil( log(corcut/alam) / log(fabs(rho)) ) ;
60 else
61 bmax = 1 ; /* pure MA(1) case */
62 } else {
63 bmax = 0 ; /* identity matrix case */
64 }
65
66 /* special and trivial case: identity matrix */
67
68 if( bmax == 0 ){
69 for( ii=0 ; ii < nt ; ii++ ){
70 len[ii] = 1 ; rc[ii] = malloc(sizeof(MTYPE)) ; rc[ii][0] = 1.0 ;
71 }
72 return rcm ;
73 }
74
75 /* First row/column has only 1 entry = diagonal value = 1 */
76
77 len[0] = 1 ; rc[0] = malloc(sizeof(MTYPE)) ; rc[0][0] = 1.0 ;
78
79 /* Subsequent rows/columns: */
80
81 for( ii=1 ; ii < nt ; ii++ ){
82 itt = TAU(ii) ; /* 'time' of the i'th index */
83 jbot = ii-bmax ; if( jbot < 0 ) jbot = 0 ; /* earliest allow index */
84 for( jj=jbot ; jj < ii ; jj++ ){ /* scan to find bandwith */
85 jtt = itt - TAU(jj) ; /* 'time' difference i-to-j */
86 if( jtt <= bmax ) break ; /* if in OK region, stop now */
87 }
88 jbot = jj ; /* this is the earliest index to be correlated with #i */
89 if( jbot == ii ){ /* a purely diagonal row/colum (inter-run gap?) */
90 len[ii] = 1 ; rc[ii] = malloc(sizeof(MTYPE)) ; rc[ii][0] = 1.0 ;
91 continue ;
92 }
93 len[ii] = ii + 1 - jbot ; /* number of entries in row/column */
94 rc[ii] = calloc(sizeof(MTYPE),len[ii]) ; /* space for the entries */
95 rii = rc[ii] - jbot ; /* shifted pointer to this row/column */
96 rii[ii] = 1.0 ; /* diagonal entry */
97 for( jj=jbot ; jj < ii ; jj++ ){ /* compute off diagonal entries */
98 jtt = itt - TAU(jj) ; /* 'time' difference again */
99 if( jtt == 1 ) rii[jj] = lam ; /* lag==1 means lam */
100 else if( jtt > 1 ) rii[jj] = lam * pow( rho , jtt-1.0 ) ;
101 }
102 }
103
104 return rcm ;
105 }
106
107 /*--------------------------------------------------------------------------*/
108 #include "armacor.c" /* ARMA(p,1) models for p=3 and 5 [01 Jul 2020] */
109 /*--------------------------------------------------------------------------*/
110
111 /*--------------------------------------------------------------------------*/
112 /*! Setup sparse banded correlation matrix (as an rcmat struct) for a
113 general shift-invariate structure, whose lagged correlations are given
114 in corvec.
115 * tau[i] is the 'true' time index of the i-th data point. This
116 lets you allow for censoring and for inter-run gaps.
117 * If tau==NULL, tau[i] is taken to be i -- that is, no censoring/gaps.
118 *//*------------------------------------------------------------------------*/
119
rcmat_arma_gen(int nt,int * tau,doublevec * corvec)120 rcmat * rcmat_arma_gen( int nt, int *tau, doublevec *corvec )
121 {
122 rcmat *rcm ;
123 LENTYP *len ;
124 MTYPE **rc , *rii , alam ;
125 int ii , jj , bmax , jbot , itt,jtt ;
126
127 if( nt < 2 || corvec == NULL ) return NULL ;
128
129 rcm = rcmat_init( nt ) ; /* create sparse matrix struct */
130 len = rcm->len ;
131 rc = rcm->rc ;
132
133 /* set maximum bandwidth */
134
135 bmax = corvec->nar - 1 ; /* we have lags 0 .. bmax */
136
137 /* special and trivial case: identity matrix */
138
139 if( bmax == 0 ){
140 for( ii=0 ; ii < nt ; ii++ ){
141 len[ii] = 1 ; rc[ii] = malloc(sizeof(MTYPE)) ; rc[ii][0] = 1.0 ;
142 }
143 return rcm ;
144 }
145
146 /* First row/column has only 1 entry = diagonal value = 1 */
147
148 len[0] = 1 ; rc[0] = malloc(sizeof(MTYPE)) ; rc[0][0] = 1.0 ;
149
150 /* Subsequent rows/columns: */
151
152 for( ii=1 ; ii < nt ; ii++ ){
153 itt = TAU(ii) ; /* 'time' of the i'th index */
154 jbot = ii-bmax ; if( jbot < 0 ) jbot = 0 ; /* earliest allow index */
155 for( jj=jbot ; jj < ii ; jj++ ){ /* scan to find bandwith */
156 jtt = itt - TAU(jj) ; /* 'time' difference i-to-j */
157 if( jtt <= bmax ) break ; /* if in OK region, stop now */
158 }
159 jbot = jj ; /* this is the earliest index to be correlated with #i */
160 if( jbot == ii ){ /* a purely diagonal row/colum (inter-run gap?) */
161 len[ii] = 1 ; rc[ii] = malloc(sizeof(MTYPE)) ; rc[ii][0] = 1.0 ;
162 continue ;
163 }
164 len[ii] = ii + 1 - jbot ; /* number of entries in row/column */
165 rc[ii] = calloc(sizeof(MTYPE),len[ii]) ; /* space for the entries */
166 rii = rc[ii] - jbot ; /* shifted pointer to this row/column */
167 rii[ii] = 1.0 ; /* diagonal entry */
168 for( jj=jbot ; jj < ii ; jj++ ){ /* compute off diagonal entries */
169 jtt = itt - TAU(jj) ; /* 'time' difference again */
170 rii[jj] = corvec->ar[jtt] ; /* extract correlation from corvec */
171 }
172 }
173
174 return rcm ;
175 }
176
177 /*---------------------------------------------------------------------------*/
178
rcmat_arma31(int nt,int * tau,double a,double r1,double t1,double vrt)179 rcmat * rcmat_arma31( int nt , int *tau ,
180 double a, double r1, double t1, double vrt )
181 {
182 rcmat *rcm ;
183 doublevec *corvec ;
184
185 corvec = arma31_correlations( a , r1 , t1 , vrt , corcut , nt ) ;
186 { int ii ; double ssum=0.0 ;
187 INFO_message("Correlation count: %d",corvec->nar) ;
188 for( ii=0 ; ii < corvec->nar ; ii++ ){
189 fprintf(stderr," %g",corvec->ar[ii]) ; if( ii > 0 ) ssum += fabs(corvec->ar[ii]) ;
190 }
191 fprintf(stderr," : Gsum = %g\n",ssum) ;
192 }
193 rcm = rcmat_arma_gen( nt, tau, corvec ) ;
194 KILL_doublevec( corvec ) ;
195 return rcm ;
196 }
197
198 #ifdef ALLOW_ARMA51
199 /*---------------------------------------------------------------------------*/
200
rcmat_arma51(int nt,int * tau,double a,double r1,double t1,double r2,double t2,double vrt)201 rcmat * rcmat_arma51( int nt , int *tau ,
202 double a, double r1, double t1,
203 double r2, double t2, double vrt )
204 {
205 rcmat *rcm ;
206 doublevec *corvec ;
207
208 corvec = arma51_correlations( a , r1 , t1 , r2 , t2 , vrt , corcut , nt ) ;
209 { int ii ; double ssum=0.0 ;
210 INFO_message("Correlation count: %d",corvec->nar) ;
211 for( ii=0 ; ii < corvec->nar ; ii++ ){
212 fprintf(stderr," %g",corvec->ar[ii]) ; if( ii > 0 ) ssum += fabs(corvec->ar[ii]) ;
213 }
214 fprintf(stderr," : Gsum = %g\n",ssum) ;
215 }
216 rcm = rcmat_arma_gen( nt, tau, corvec ) ;
217 KILL_doublevec( corvec ) ;
218 return rcm ;
219 }
220 #endif
221
222 /*---------------------------------------------------------------------------*/
223
main(int argc,char * argv[])224 int main( int argc , char *argv[] )
225 {
226 int nlen=0 , nvec=1 , iarg=1 , nbad=0 , kk,ii , do_norm=0 ;
227 double aa=-666.0, bb=-666.0 , lam=-666.0 , sig=1.0 ;
228 double *rvec ;
229 rcmat *rcm=NULL ;
230 MRI_IMAGE *outim ;
231 float *outar , *vv ;
232 long seed=0 ;
233
234 int do_arma31 = 0 , do_arma51 = 0 , do_arma11 = 1 ;
235 double r1=0.0,t1=0.0 , r2=0.0,t2=0.0 , vrt=0.0 ;
236
237 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
238 printf(
239 "Program to generate an ARMA(1,1) time series, for simulation studies.\n"
240 "Results are written to stdout.\n"
241 "\n"
242 "Usage: 1dgenARMA11 [options]\n"
243 "\n"
244 "Options:\n"
245 "========\n"
246 " -num N } These equivalent options specify the length of the time\n"
247 " -len N } series vector to generate.\n"
248 "\n"
249 " -nvec M = The number of time series vectors to generate;\n"
250 " if this option is not given, defaults to 1.\n"
251 "\n"
252 " -a a = Specify ARMA(1,1) parameters 'a'.\n"
253 " -b b = Specify ARMA(1,1) parameter 'b' directly.\n"
254 " -lam lam = Specify ARMA(1,1) parameter 'b' indirectly.\n"
255 "\n"
256 " -sig ss = Set standard deviation of results [default=1].\n"
257 " -norm = Normalize time series so sum of squares is 1.\n"
258 " -seed dd = Set random number seed.\n"
259 "\n"
260 " * The correlation coefficient r(k) of noise samples k units apart in time,\n"
261 " for k >= 1, is given by r(k) = lam * a^(k-1)\n"
262 " where lam = (b+a)(1+a*b)/(1+2*a*b+b*b)\n"
263 " (N.B.: lam=a when b=0 -- AR(1) noise has r(k)=a^k for k >= 0)\n"
264 " (N.B.: lam=b when a=0 -- MA(1) noise has r(k)=b for k=1, r(k)=0 for k>1)\n"
265 " * lam can be bigger or smaller than a, depending on the sign of b:\n"
266 " b > 0 means lam > a; b < 0 means lam < a.\n"
267 " * What I call (a,b) here is sometimes called (p,q) in the ARMA literature.\n"
268 " * For a noise model which is the sum of AR(1) and white noise, 0 < lam < a\n"
269 " (i.e., a > 0 and -a < b < 0 ).\n"
270 "\n"
271 " -CORcut cc = The exact ARMA(1,1) correlation matrix (for a != 0)\n"
272 " has no non-zero entries. The calculations in this\n"
273 " program set correlations below a cutoff to zero.\n"
274 " The default cutoff is %.5f, but can be altered with\n"
275 " this option. The usual reason to use this option is\n"
276 " to test the sensitivity of the results to the cutoff.\n"
277 "\n"
278 "-----------------------------\n"
279 "A restricted ARMA(3,1) model:\n"
280 "-----------------------------\n"
281 "Skip the '-a', '-b', and '-lam' options, and use a model with 3 roots\n"
282 "\n"
283 " -arma31 a r theta vrat\n"
284 "\n"
285 " where the roots are z = a, z = r*exp(I*theta), z = r*exp(-I*theta)\n"
286 " and vrat = s^2/(s^2+w^2) [so 0 < vrat < 1], where s = variance\n"
287 " of the pure AR(3) component and w = variance of extra white noise\n"
288 " added to the AR(3) process -- this is the 'restricted' ARMA(3,1).\n"
289 "\n"
290 " If the data has given TR, and you want a frequency of f Hz, in\n"
291 " the noise model, then theta = 2 * PI * TR * f. If theta > PI,\n"
292 " then you are modeling noise beyond the Nyquist frequency and\n"
293 " the gods (and this program) won't be happy.\n"
294 "\n"
295 " # csh syntax for 'set' variable assignment commands\n"
296 " set nt = 500\n"
297 " set tr = 1\n"
298 " set df = `ccalc \"1/($nt*$tr)\"`\n"
299 " set f1 = 0.10\n"
300 " set t1 = `ccalc \"2*PI*$tr*$f1\"`\n"
301 " 1dgenARMA11 -nvec 500 -len $nt -arma31 0.8 0.9 $t1 0.9 -CORcut 0.0001 \\\n"
302 " | 1dfft -nodetrend stdin: > qqq.1D\n"
303 " 3dTstat -mean -prefix stdout: qqq.1D \\\n"
304 " | 1dplot -stdin -num 201 -dt $df -xlabel 'frequency' -ylabel '|FFT|'\n"
305 "---------------------------------------------------------------------------\n"
306 #ifdef ALLOW_ARMA51
307 "A similar option is now available for a restricted ARMA(5,1) model:\n"
308 " -arma51 a r1 theta1 r2 theta2 vrat\n"
309 "where now the roots are\n"
310 " z = a z = r1*exp(I*theta1) z = r1*exp(-I*theta1)\n"
311 " z = r2*exp(I*theta2) z = r2*exp(-I*theta2)\n"
312 "This model allows the simulation of two separate frequencies in the 'noise'.\n"
313 "---------------------------------------------------------------------------\n"
314 #endif
315 "\n"
316 "Author: RWCox [for his own demented and deranged purposes]\n"
317 "\n"
318 "Examples:\n"
319 " 1dgenARMA11 -num 200 -a .8 -lam 0.7 | 1dplot -stdin\n"
320 " 1dgenARMA11 -num 2000 -a .8 -lam 0.7 | 1dfft -nodetrend stdin: stdout: | 1dplot -stdin\n"
321 "\n" ,
322 corcut
323 ) ;
324 exit(0) ;
325 }
326
327 while( iarg < argc ){
328
329 if( strcasecmp(argv[iarg],"-ARMA31") == 0 ){ /* 01 Jul 2020 */
330 if( iarg+4 >= argc ) ERROR_exit("Need 4 arguments after option '%s'",argv[iarg]) ;
331 do_arma31 = 1 ; do_arma51 = do_arma11 = 0 ;
332 aa = strtod( argv[++iarg] , NULL ) ;
333 r1 = strtod( argv[++iarg] , NULL ) ;
334 t1 = strtod( argv[++iarg] , NULL ) ;
335 vrt = strtod( argv[++iarg] , NULL ) ;
336 iarg++ ; continue ;
337 }
338
339 #ifdef ALLOW_ARMA51
340 if( strcasecmp(argv[iarg],"-ARMA51") == 0 ){ /* 01 Jul 2020 */
341 if( iarg+6 >= argc ) ERROR_exit("Need 6 arguments after option '%s'",argv[iarg]) ;
342 do_arma11 = do_arma31 = 0 ; do_arma51 = 1 ;
343 aa = strtod( argv[++iarg] , NULL ) ;
344 r1 = strtod( argv[++iarg] , NULL ) ;
345 t1 = strtod( argv[++iarg] , NULL ) ;
346 r2 = strtod( argv[++iarg] , NULL ) ;
347 t2 = strtod( argv[++iarg] , NULL ) ;
348 vrt = strtod( argv[++iarg] , NULL ) ;
349 iarg++ ; continue ;
350 }
351 #endif
352
353 if( strcmp(argv[iarg],"-norm") == 0 ){
354 do_norm = 1 ; iarg++ ; continue ;
355 }
356
357 if( strcmp(argv[iarg],"-num") == 0 || strcmp(argv[iarg],"-len") == 0 ){
358 iarg++ ;
359 if( iarg >= argc ) ERROR_exit("Need argument after option '%s'",argv[iarg-1]) ;
360 nlen = (int)strtod(argv[iarg],NULL) ;
361 if( nlen < 3 ) ERROR_exit("Time series length must be at least 3!") ;
362 iarg++ ; continue ;
363 }
364
365 if( strcmp(argv[iarg],"-nvec") == 0 ){
366 iarg++ ;
367 if( iarg >= argc ) ERROR_exit("Need argument after option '%s'",argv[iarg-1]) ;
368 nvec = (int)strtod(argv[iarg],NULL) ;
369 if( nvec < 1 ) ERROR_exit("Number of vectors must be at least 1!") ;
370 iarg++ ; continue ;
371 }
372
373 if( strcmp(argv[iarg],"-seed") == 0 ){
374 iarg++ ;
375 if( iarg >= argc ) ERROR_exit("Need argument after option '%s'",argv[iarg-1]) ;
376 seed = (long)strtod(argv[iarg],NULL) ;
377 iarg++ ; continue ;
378 }
379
380 if( strcmp(argv[iarg],"-sig") == 0 ){
381 iarg++ ;
382 if( iarg >= argc ) ERROR_exit("Need argument after option '%s'",argv[iarg-1]) ;
383 sig = strtod(argv[iarg],NULL) ;
384 if( sig <= 0.0 ) ERROR_exit("Value after '-sig' should be positive!") ;
385 iarg++ ; continue ;
386 }
387
388 if( strcmp(argv[iarg],"-a") == 0 ){
389 iarg++ ;
390 if( iarg >= argc ) ERROR_exit("Need argument after option '%s'",argv[iarg-1]) ;
391 aa = strtod(argv[iarg],NULL) ;
392 if( fabs(aa) >= 0.95 ) ERROR_exit("Largest allowed value of a = 0.95") ;
393 iarg++ ; continue ;
394 }
395
396 if( strcmp(argv[iarg],"-b") == 0 ){
397 iarg++ ;
398 if( iarg >= argc ) ERROR_exit("Need argument after option '%s'",argv[iarg-1]) ;
399 bb = strtod(argv[iarg],NULL) ;
400 if( fabs(bb) >= 0.95 ) ERROR_exit("Largest allowed value of b = 0.95") ;
401 iarg++ ; continue ;
402 }
403
404 if( strcmp(argv[iarg],"-lam") == 0 ){
405 iarg++ ;
406 if( iarg >= argc ) ERROR_exit("Need argument after option '%s'",argv[iarg-1]) ;
407 lam = strtod(argv[iarg],NULL) ;
408 if( fabs(lam) >= 0.95 ) ERROR_exit("Largest allowed value of lam = 0.95") ;
409 iarg++ ; continue ;
410 }
411
412 if( strcmp(argv[iarg],"-CORcut") == 0 ){
413 iarg++ ;
414 if( iarg >= argc ) ERROR_exit("Need argument after option '%s'",argv[iarg-1]) ;
415 corcut = strtod(argv[iarg],NULL) ;
416 if( corcut <= 0.0 || corcut > 0.1 ) ERROR_exit("Illegal value after '-CORcut'") ;
417 iarg++ ; continue ;
418 }
419
420 ERROR_exit("Unrecognized option: '%s'",argv[iarg]) ;
421 }
422
423 /*-- check errors --*/
424
425 if( nlen <= 3 ){ ERROR_message("Didn't you give the -num/-len option?"); nbad++; }
426
427 if( do_arma11 ){
428 if( aa <= -1.0 ){ ERROR_message("Didn't you give the -a option?"); nbad++; }
429 if( bb <= -1.0 && lam <= -1.0 ){
430 ERROR_message("Didn't you give the -b or -lam option?"); nbad++;
431 }
432 if( nbad > 0 ) ERROR_exit("Can't continue past the above problem%s" ,
433 (nbad==1)? "\0" : "s" ) ;
434
435 if( bb > -1.0 ){
436 lam = LAMBDA(aa,bb) ; INFO_message("lam(a=%.3f,b=%.3f) = %.3f",aa,bb,lam) ;
437 }
438
439 /* setup */
440
441 rcm = rcmat_arma11( nlen , NULL , aa , lam ) ;
442
443 } else if( do_arma31 ){
444
445 rcm = rcmat_arma31( nlen , NULL , aa,r1,t1,vrt ) ;
446
447 #ifdef ALLOW_ARMA51
448 } else if( do_arma51 ){
449
450 rcm = rcmat_arma51( nlen , NULL , aa,r1,t1,r2,t2,vrt ) ;
451 #endif
452
453 }
454
455 if( rcm == NULL ) ERROR_exit("Can't setup matrix?!") ; /* should be impossible */
456
457 kk = rcmat_choleski( rcm ) ;
458 if( kk > 0 ) ERROR_exit("Choleski fails at row %d",kk) ;
459
460 /* simulate */
461
462 outim = mri_new( nlen , nvec , MRI_float ) ; outar = MRI_FLOAT_PTR(outim) ;
463 rvec = (double *)malloc(sizeof(double)*nlen) ;
464
465 if( seed == 0 ) seed = (long)time(NULL)+(long)getpid() ;
466 srand48(seed) ;
467
468 for( kk=0 ; kk < nvec ; kk++ ){
469 for( ii=0 ; ii < nlen ; ii++ ) rvec[ii] = zgaussian() ;
470 rcmat_lowert_vecmul( rcm , rvec ) ;
471 vv = outar + kk*nlen ;
472 if( do_norm ){
473 sig = 0.0 ;
474 for( ii=0 ; ii < nlen ; ii++ ) sig += rvec[ii]*rvec[ii] ;
475 sig = 1.0 / sqrt(sig) ;
476 }
477 for( ii=0 ; ii < nlen ; ii++ ) vv[ii] = sig * rvec[ii] ;
478 }
479
480 mri_write_1D(NULL,outim) ;
481 exit(0) ;
482 }
483