1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "afni.h"
8
9 #ifndef ALLOW_PLUGINS
10 # error "Plugins not properly set up -- see machdep.h"
11 #endif
12
13 /***********************************************************************
14 Plugin to provide a least squares fitting 1D function for graphs
15 ************************************************************************/
16
17 /*------------- string to 'help' the user -------------*/
18
19 static char helpstring[] =
20 " Purpose: Control the 'LSqFit' and 'LSqDtr 1D functions.\n"
21 "\n"
22 " Parameters: Baseline = 'Constant' or 'Linear'\n"
23 " Is the baseline 'a' or 'a+b*t'?\n"
24 " Ignore = Number of points to ignore at\n"
25 " start of each timeseries.\n"
26 " \n"
27 " Sinusoids: Period = Fundamental period to use.\n"
28 " Harmonics = Number of overtones to use.\n"
29 " \n"
30 " Timeseries: File = Input timeseries file to use.\n"
31 ;
32
33 static char plehstring[] =
34 " Purpose: Generate a timeseries and store in AFNI list\n"
35 "\n"
36 " Parameters: Delta = time step between points\n"
37 " Length = number of points\n"
38 "\n"
39 " Output: Label = String to label timeseries by\n"
40 " in AFNI choosers\n"
41 "\n"
42 " Polynomial: Order = Maximum power to include\n"
43 " (Chebyshev polynomials are used)\n"
44 "\n"
45 " Sinusoid: Period = Fundamental period to use.\n"
46 " Harmonics = Number of overtones to use.\n"
47 ;
48
49 /*------- Strings for baseline control ------*/
50
51 #define NBASE 3
52 static char * baseline_strings[NBASE] = { "Constant" , "Linear" , "Quadratic" } ;
53
54 /*--------------- prototypes for internal routines ---------------*/
55
56 char * LSQ_main( PLUGIN_interface * ) ; /* the entry point */
57
58 void LSQ_fitter( int nt, double to, double dt, float * vec, char ** label ) ;
59 void LSQ_detrend( int nt, double to, double dt, float * vec, char ** label ) ;
60 void LSQ_worker( int nt, double dt, float * vec, int dofit, char ** label ) ;
61
62 PLUGIN_interface * TSGEN_init(void) ;
63 char * TSGEN_main( PLUGIN_interface * ) ;
64
65 PLUGIN_interface * EXP0D_init(void) ;
66 char * EXP0D_main( PLUGIN_interface * ) ;
67 void EXP0D_worker( int num , float * vec ) ;
68
69 #ifdef ALLOW_LOMO
70 PLUGIN_interface * LOMOR_init(void) ;
71 char * LOMOR_main( PLUGIN_interface * ) ;
72 void LOMOR_fitter() ;
73 #endif
74
75 /*---------------- global data -------------------*/
76
77 static PLUGIN_interface * global_plint = NULL ;
78
79 #define NRMAX_SIN 2
80 #define NRMAX_TS 2
81 #define HARM_MAX 22
82
83 static int polort=1 , ignore=3 , nrsin=0 , nrts=0,ntsim=0 , initialize=1 ;
84 static float sinper[NRMAX_SIN] ;
85 static int sinharm[NRMAX_SIN] ;
86 static MRI_IMAGE * tsim[NRMAX_TS] ;
87
88 /***********************************************************************
89 Set up the interface to the user:
90 1) Create a new interface using "PLUTO_new_interface";
91
92 2) For each line of inputs, create the line with "PLUTO_add_option"
93 (this line of inputs can be optional or mandatory);
94
95 3) For each item on the line, create the item with
96 "PLUTO_add_dataset" for a dataset chooser,
97 "PLUTO_add_string" for a string chooser,
98 "PLUTO_add_number" for a number chooser,
99 "PLUTO_add_timeseries" for a timeseries chooser.
100 ************************************************************************/
101
102
103 DEFINE_PLUGIN_PROTOTYPE
104
PLUGIN_init(int ncall)105 PLUGIN_interface * PLUGIN_init( int ncall )
106 {
107 int ii ;
108 PLUGIN_interface * plint ; /* will be the output of this routine */
109
110 if( ncall > 3 ) return NULL ; /* generate interfaces for ncall 0-3 */
111
112 if( ncall == 1 ) return TSGEN_init() ; /* interface # 1 */
113 if( ncall == 2 ) return EXP0D_init() ; /* interface # 2 */
114 #ifdef ALLOW_LOMO
115 if( ncall == 3 ) return LOMOR_init() ; /* interface # 3 */
116 #else
117 if( ncall == 3 ) return NULL ;
118 #endif
119
120 /***** otherwise, do interface # 0 *****/
121
122 /*---------------- set titles and call point ----------------*/
123
124 CHECK_IF_ALLOWED("L2FIT","LSqFit & Dtr") ; /* 30 Sep 2016 */
125
126 plint = PLUTO_new_interface( "LSqFit & Dtr" ,
127 "Control LSqFit and LSqDtr Functions" ,
128 helpstring ,
129 PLUGIN_CALL_VIA_MENU , LSQ_main ) ;
130
131 global_plint = plint ; /* make global copy */
132
133 PLUTO_set_sequence( plint , "A:funcs:fitting" ) ;
134
135 PLUTO_add_hint( plint , "Control LSqFit and LSqDtr Functions" ) ;
136
137 PLUTO_set_runlabels( plint , "Set+Keep" , "Set+Close" ) ; /* 04 Nov 2003 */
138
139 /*----- Parameters -----*/
140
141 PLUTO_add_option( plint , "Parameters" , "Parameters" , TRUE ) ;
142
143 PLUTO_add_string( plint , "Baseline" , NBASE , baseline_strings , 1 ) ;
144
145 PLUTO_add_number( plint , "Ignore" , 0,20,0,3 , FALSE ) ;
146
147 /*----- Sinusoid -----*/
148
149 for( ii=0 ; ii < NRMAX_SIN ; ii++ ){
150 PLUTO_add_option( plint , "Sinusoid" , "Sinusoid" , FALSE ) ;
151 PLUTO_add_number( plint , "Period" , 0,99999,0,20, TRUE ) ;
152 PLUTO_add_number( plint , "Harmonics" , 1,HARM_MAX,0,1 , FALSE ) ;
153 }
154
155 /*----- Timeseries -----*/
156
157 for( ii=0 ; ii < NRMAX_TS ; ii++ ){
158 PLUTO_add_option( plint , "Timeseries" , "Timeseries" , FALSE ) ;
159 PLUTO_add_timeseries( plint , "File" ) ;
160 }
161
162 /*--------- done with interface setup ---------*/
163
164 PLUTO_register_1D_funcstr( "LSqFit" , LSQ_fitter ) ;
165 PLUTO_register_1D_funcstr( "LSqDtr" , LSQ_detrend ) ;
166
167 return plint ;
168 }
169
170 /***************************************************************************
171 Main routine for this plugin (will be called from AFNI).
172 If the return string is not NULL, some error transpired, and
173 AFNI will popup the return string in a message box.
174 ****************************************************************************/
175
LSQ_main(PLUGIN_interface * plint)176 char * LSQ_main( PLUGIN_interface * plint )
177 {
178 char * str ;
179 int ii ;
180 float * tsar ;
181
182 /*--------- go to first input line ---------*/
183
184 PLUTO_next_option(plint) ;
185
186 str = PLUTO_get_string(plint) ;
187 polort = PLUTO_string_index( str , NBASE , baseline_strings ) ;
188
189 ignore = PLUTO_get_number(plint) ;
190
191 /*------ loop over remaining options, check their tags, process them -----*/
192
193 nrsin = nrts = ntsim = 0 ;
194 do {
195 str = PLUTO_get_optiontag(plint) ; if( str == NULL ) break ;
196
197 if( strcmp(str,"Sinusoid") == 0 ){
198
199 sinper[nrsin] = PLUTO_get_number(plint) ;
200 sinharm[nrsin] = PLUTO_get_number(plint) - 1.0 ;
201 if( sinper[nrsin] <= 0.0 )
202 return "************************\n"
203 "Illegal Sinusoid Period!\n"
204 "************************" ;
205
206 nrsin++ ;
207
208 } else if( strcmp(str,"Timeseries") == 0 ){
209
210 tsim[ntsim] = PLUTO_get_timeseries(plint) ;
211
212 if( tsim[ntsim] == NULL || tsim[ntsim]->nx < 3 || tsim[ntsim]->kind != MRI_float )
213 return "*************************\n"
214 "Illegal Timeseries Input!\n"
215 "*************************" ;
216
217 tsar = MRI_FLOAT_PTR(tsim[ntsim]) ;
218 for( ii=ignore ; ii < tsim[ntsim]->nx && tsar[ii] >= WAY_BIG ; ii++ ) ; /* nada */
219 ignore = ii ;
220 nrts += tsim[ntsim]->ny ; ntsim++ ;
221
222 } else {
223 return "************************\n"
224 "Illegal optiontag found!\n"
225 "************************" ;
226 }
227 } while(1) ;
228
229 /*--- nothing left to do until data arrives ---*/
230
231 initialize = 1 ; /* force re-initialization */
232
233 /*** compute how many ref functions are ordered ***/
234
235 { int nref , ks ;
236 char str[64] ;
237
238 nref = (polort+1) + nrts ;
239 for( ks=0 ; ks < nrsin ; ks++ ) nref += 2*(sinharm[ks]+1) ;
240 sprintf(str," \nNumber of fit parameters = %d\n",nref) ;
241 PLUTO_popup_transient( plint , str ) ;
242 }
243
244 return NULL ;
245 }
246
247 /*---------------------------------------------------------------*/
248
249 /** 22 Apr 1997: added label that will go to graphs **/
250
LSQ_fitter(int nt,double to,double dt,float * vec,char ** label)251 void LSQ_fitter( int nt , double to , double dt , float * vec , char ** label )
252 {
253 LSQ_worker( nt , dt , vec , TRUE , label ) ;
254 return ;
255 }
256
LSQ_detrend(int nt,double to,double dt,float * vec,char ** label)257 void LSQ_detrend( int nt , double to , double dt , float * vec , char ** label )
258 {
259 LSQ_worker( nt , dt , vec , FALSE , label ) ;
260 return ;
261 }
262
263 static char lbuf[8192] ; /* 22 Apr 1997: will hold label for graphs */
264 static char sbuf[256] ;
265
LSQ_worker(int nt,double dt,float * vec,int dofit,char ** label)266 void LSQ_worker( int nt , double dt , float * vec , int dofit , char ** label )
267 {
268 int nlen , nref ;
269
270 static int nlen_old = -666 , nref_old = -666 ;
271 static double dt_old = -666.666 ;
272 static float ** ref = NULL ;
273 static double * dch = NULL ;
274
275 int ir , ii , ks,jh , j;
276 float fac , tm , val ;
277 float * fit , * tsar , *qsar ;
278
279 /*** compute how many ref functions are ordered ***/
280
281 nref = (polort+1) + nrts ;
282 for( ks=0 ; ks < nrsin ; ks++ ) nref += 2*(sinharm[ks]+1) ;
283
284 /*** do nothing if not enough data to fit ***/
285
286 nlen = nt - ignore ;
287
288 if( nlen <= nref ) return ; /* do nothing if not enough data to fit */
289
290 /** if data vectors are new length,
291 or have a new number of reference vectors,
292 or have a new time step and need sinusoids,
293 or the initialize flag is set,
294 then reinitialize reference vectors and Choleski factor **/
295
296 if( nlen != nlen_old || nref != nref_old ||
297 initialize || (dt != dt_old && nrsin > 0) ){
298
299 /* free old storage */
300
301 if( ref != NULL ){
302 for( ir=0 ; ir < nref_old ; ir++ ) if( ref[ir] != NULL ) free(ref[ir]) ;
303 free(ref) ;
304 }
305 if( dch != NULL ) free(dch) ;
306
307 /* make space for ref vectors */
308
309 ref = (float **) malloc( sizeof(float *) * nref ) ;
310 if( ref == NULL ){fprintf(stderr,"\nmalloc error in plug_lsqfit\n\a");
311 return;
312 /* EXIT(1); */
313 }
314 for( ir=0 ; ir < nref ; ir++ ){
315 ref[ir] = (float *) malloc( sizeof(float) * nlen ) ;
316 if( ref[ir] == NULL ){
317 fprintf(stderr,"\nmalloc error in plug_lsqfit\n\a");
318 for(j=0;j<=ir;j++)
319 free(ref[j]);
320 free(ref);
321 ref = NULL;
322 return;
323 /* EXIT(1); */
324 }
325 }
326 nlen_old = nlen ;
327 nref_old = nref ;
328 dt_old = dt ;
329
330 /**** fill ref vectors ****/
331
332 /* r(t) = 1 */
333
334 for( ii=0 ; ii < nlen ; ii++ ) ref[0][ii] = 1.0 ;
335
336 ir = 1 ;
337 if( polort > 0 ){
338
339 /* r(t) = t - tmid */
340
341 tm = 0.5 * (nlen-1.0) ; fac = 2.0 / nlen ;
342 for( ii=0 ; ii < nlen ; ii++ ) ref[1][ii] = (ii-tm)*fac ;
343 ir = 2 ;
344
345 /* r(t) = (t-tmid)**ir */
346
347 for( ; ir <= polort ; ir++ )
348 for( ii=0 ; ii < nlen ; ii++ )
349 ref[ir][ii] = pow( (ii-tm)*fac , (double)ir ) ;
350 }
351
352 if( dt == 0.0 ) dt = 1.0 ;
353
354 /* r(t) = sinusoids */
355
356 for( ks=0 ; ks < nrsin ; ks++ ){
357 for( jh=0 ; jh <= sinharm[ks] ; jh++ ){
358 fac = (2.0*PI) * dt * (jh+1) / sinper[ks] ;
359 for( ii=0 ; ii < nlen ; ii++ ){
360 ref[ir] [ii] = cos( fac * ii ) ;
361 ref[ir+1][ii] = sin( fac * ii ) ;
362 }
363 ir += 2 ;
364 }
365 }
366
367 /* r(t) = timeseries files */
368
369 for( ks=0 ; ks < ntsim ; ks++ ){
370 if( tsim[ks] == NULL || tsim[ks]->nx - ignore < nlen ){
371 initialize = 1 ;
372 fprintf(stderr,"Inadequate time series #%d in LSQ plugin\n\a",ks+1) ;
373 return ;
374 }
375 tsar = MRI_FLOAT_PTR(tsim[ks]) ;
376 for( jh=0 ; jh < tsim[ks]->ny ; jh++ ){
377 qsar = tsar + jh*tsim[ks]->nx ;
378 for( ii=0 ; ii < nlen ; ii++ ) ref[ir][ii] = qsar[ii+ignore] ;
379 ir++ ;
380 }
381 }
382
383 /* Cholesky-ize */
384
385 dch = startup_lsqfit( nlen , NULL , nref , ref ) ;
386 if( dch == NULL ){
387 initialize = 1 ;
388 fprintf(stderr,"Choleski error in LSQ plugin\n\a") ;
389 return ;
390 }
391
392 initialize = 0 ;
393 }
394
395 /** find least squares fit coefficients **/
396
397 fit = delayed_lsqfit( nlen , vec+ignore , nref , ref , dch ) ;
398
399 for( ii=0 ; ii < nlen ; ii++ ){
400 val = 0.0 ;
401 for( ir=0 ; ir < nref ; ir++ ) val += fit[ir] * ref[ir][ii] ;
402
403 vec[ii+ignore] = (dofit) ? val : vec[ii+ignore] - val ;
404 }
405
406 /** 22 Apr 1997: create label if desired by AFNI **/
407 /** [This is in static storage, since AFNI will copy it] **/
408
409 if( label != NULL ){ /* assemble this 1 line at a time from sbuf */
410
411 lbuf[0] = '\0' ; /* make this a 0 length string to start */
412
413 /** for each reference, make a string into sbuf **/
414
415 ir = 0 ;
416 sprintf(sbuf,"Coef of 1 = %g\n" , fit[ir++] ) ;
417 strcat(lbuf,sbuf) ;
418
419 for( ; ir <= polort ; ir++ ){
420 sprintf(sbuf,"Coef of t**%d = %g\n" , ir,fit[ir] ) ;
421 strcat(lbuf,sbuf) ;
422 }
423
424 for( ks=0 ; ks < nrsin ; ks++ ){
425 for( jh=0 ; jh <= sinharm[ks] ; jh++ ){
426 fac = sinper[ks] / (jh+1) ;
427 sprintf(sbuf,"Coef of cos(2*Pi*t/%g) = %g\n" , fac , fit[ir++] ) ;
428 strcat(lbuf,sbuf) ;
429 sprintf(sbuf,"Coef of sin(2*Pi*t/%g) = %g\n" , fac , fit[ir++] ) ;
430 strcat(lbuf,sbuf) ;
431 }
432 }
433
434 for( ks=0 ; ks < ntsim ; ks++ ){
435 for( jh=0 ; jh < tsim[ks]->ny ; jh++ ){
436 sprintf(sbuf,"Coef of %s[%d] = %g\n" , tsim[ks]->name,jh , fit[ir++] ) ;
437 strcat(lbuf,sbuf) ;
438 }
439 }
440
441 *label = lbuf ; /* send address of lbuf back in what label points to */
442 }
443
444 free(fit) ;
445 return ;
446 }
447
448 /*********************************************************************************/
449
TSGEN_init(void)450 PLUGIN_interface * TSGEN_init(void)
451 {
452 PLUGIN_interface * plint ;
453
454 /*---------------- set titles and call point ----------------*/
455
456 CHECK_IF_ALLOWED("TSGENERATE","TS Generate") ;
457
458 plint = PLUTO_new_interface( "TS Generate" ,
459 "Generate a Timeseries" ,
460 plehstring ,
461 PLUGIN_CALL_VIA_MENU , TSGEN_main ) ;
462
463 PLUTO_add_hint( plint , "Generate a 1D Timeseries" ) ;
464
465 /*----- Parameters -----*/
466
467 PLUTO_add_option( plint , "Parameters" , "Parameters" , TRUE ) ;
468 PLUTO_add_number( plint , "Delta" , 0,99999,1, 0 , TRUE ) ;
469 PLUTO_add_number( plint , "Length" , 3,9999,0,3 , TRUE ) ;
470
471 /*----- Output -----*/
472
473 PLUTO_add_option( plint , "Output" , "Output" , TRUE ) ;
474 PLUTO_add_string( plint , "Label" , 0,NULL , 19 ) ;
475
476 /*----- Polynomial -----*/
477
478 PLUTO_add_option( plint , "Polynomial" , "Polynomial" , FALSE ) ;
479 PLUTO_add_number( plint , "Order" , 2,20,0,2 , FALSE ) ;
480
481 /*----- Sinusoid -----*/
482
483 PLUTO_add_option( plint , "Sinusoid" , "Sinusoid" , FALSE ) ;
484 PLUTO_add_number( plint , "Period" , 0,99999,0,20, TRUE ) ;
485 PLUTO_add_number( plint , "Harmonics" , 1,HARM_MAX,0,1 , FALSE ) ;
486
487 return plint ;
488 }
489
TSGEN_main(PLUGIN_interface * plint)490 char * TSGEN_main( PLUGIN_interface * plint )
491 {
492 char * label , * str ;
493 int ii , jj ;
494 float * tsar ;
495 MRI_IMAGE * tsim ;
496 float delta , period=0.0 ;
497 int nx , ny=0 , npol=0 , nharm=-1 ;
498 int pp ;
499 double fac , val ;
500
501 /*--------- go to first input line ---------*/
502
503 PLUTO_next_option(plint) ;
504
505 delta = PLUTO_get_number(plint) ;
506 if( delta <= 0.0 ) return "**********************\n"
507 "Illegal value of Delta\n"
508 "**********************" ;
509
510 nx = PLUTO_get_number(plint) ;
511
512 /*----- next input line -----*/
513
514 PLUTO_next_option(plint) ;
515 label = PLUTO_get_string(plint) ;
516 if( label == NULL || strlen(label) == 0 ) return "**********************\n"
517 "Illegal value of Label\n"
518 "**********************" ;
519
520 /*----- rest of input lines -----*/
521
522 do {
523 str = PLUTO_get_optiontag(plint) ; if( str == NULL ) break ;
524
525 if( strcmp(str,"Sinusoid") == 0 ){
526
527 period = PLUTO_get_number(plint) ;
528 nharm = PLUTO_get_number(plint) - 1.0 ;
529
530 if( period <= 0.0 ) return "***********************\n"
531 "Illegal Sinusoid Period\n"
532 "***********************" ;
533
534
535 } else if( strcmp(str,"Polynomial") == 0 ){
536
537 npol = PLUTO_get_number(plint) ;
538
539 } else {
540 return "***********************\n"
541 "Illegal optiontag found\n"
542 "***********************" ;
543 }
544 } while(1) ;
545
546 /********** Make the timeseries ***********/
547
548 ny = 0 ;
549 if( npol > 0 ) ny = npol-1 ;
550 if( period > 0.0 ) ny += 2*(nharm+1) ;
551
552 if( ny < 1 ) return "***********************\n"
553 "No timeseries specified\n"
554 "***********************" ;
555
556 tsim = mri_new( nx , ny , MRI_float ) ;
557 jj = 0 ;
558
559 fac = 1.99999 / (nx-1) ;
560 for( pp=2 ; pp <= npol ; pp++,jj++ ){
561
562 tsar = MRI_FLOAT_PTR(tsim) + (jj*nx) ;
563
564 for( ii=0 ; ii < nx ; ii++ ){
565 val = fac * ii - 0.999995 ;
566 tsar[ii] = cos( pp * acos(val) ) ;
567 }
568 }
569
570 for( pp=0 ; pp <= nharm ; pp++ , jj+=2 ){
571 fac = (2.0*PI) * delta * (pp+1) / period ;
572 tsar = MRI_FLOAT_PTR(tsim) + (jj*nx) ;
573
574 for( ii=0 ; ii < nx ; ii++ ){
575 tsar[ii] = cos( fac * ii ) ;
576 tsar[ii+nx] = sin( fac * ii ) ;
577 }
578 }
579
580 PLUTO_register_timeseries( label , tsim ) ;
581 mri_free(tsim) ;
582 return NULL ;
583 }
584
585 /***************************************************************************/
586
587 #include "parser.h"
588
589 static char fredstring[] =
590 " Purpose: Control the Expr 0D Transformation function\n"
591 "\n"
592 " Variable = letter used as input to expression\n"
593 " Expression = arithmetic expression to evaluate\n"
594 ;
595
596 #define NALPHA 26
597 static char * vstring[NALPHA] = {
598 " A ", " B ", " C ", " D ", " E ",
599 " F ", " G ", " H ", " I ", " J ",
600 " K ", " L ", " M ", " N ", " O ",
601 " P ", " Q ", " R ", " S ", " T ",
602 " U ", " V ", " W ", " X ", " Y ", " Z " } ;
603
604 static int exp0d_var = 23 ;
605
606 static PARSER_code * exp0d_pc = NULL ;
607
608 static PLUGIN_interface *plint_EXP0D=NULL ;
609
EXP0D_func_init(void)610 static void EXP0D_func_init(void) /* 21 Jul 2003 */
611 {
612 PLUG_startup_plugin_CB( NULL , (XtPointer)plint_EXP0D , NULL ) ;
613 }
614
615
EXP0D_init(void)616 PLUGIN_interface * EXP0D_init(void)
617 {
618 PLUGIN_interface * plint ;
619
620 /*---------------- set titles and call point ----------------*/
621
622 plint = PLUTO_new_interface( "Expr 0D" ,
623 "Control the Expr 0D transformation" ,
624 fredstring ,
625 PLUGIN_CALL_VIA_MENU , EXP0D_main ) ;
626
627 PLUTO_add_option( plint , "Variable" , "Variable" , TRUE ) ;
628 PLUTO_add_string( plint , NULL , NALPHA,vstring , exp0d_var ) ;
629
630 PLUTO_add_option( plint , "Expression" , "Expression" , TRUE ) ;
631 PLUTO_add_string( plint , NULL , 0,NULL , 50 ) ;
632
633 PLUTO_register_0D_function( "Expr 0D" , EXP0D_worker ) ;
634
635 plint_EXP0D = plint ;
636 AFNI_register_nD_func_init( 0 , (generic_func *)EXP0D_func_init ) ; /* 21 Jul 2003 */
637
638 return plint ;
639 }
640
EXP0D_main(PLUGIN_interface * plint)641 char * EXP0D_main( PLUGIN_interface * plint )
642 {
643 char * str ;
644
645 PLUTO_next_option(plint) ;
646 str = PLUTO_get_string(plint) ;
647 exp0d_var = PLUTO_string_index( str , NALPHA , vstring ) ;
648
649 if( exp0d_pc != NULL ){ free(exp0d_pc) ; exp0d_pc = NULL ; }
650 PLUTO_next_option(plint) ;
651 str = PLUTO_get_string(plint) ;
652 exp0d_pc = PARSER_generate_code( str ) ;
653
654 if( exp0d_pc == NULL ) return "*******************************\n"
655 "Error when compiling expression\n"
656 "*******************************" ;
657
658 return NULL ;
659 }
660
661 #define VSIZE 64
662
EXP0D_worker(int num,float * vec)663 void EXP0D_worker( int num , float * vec )
664 {
665 int ii , jj , jbot,jtop ;
666
667 static int first = 1 ;
668 static double * atoz[26] ;
669 static double tvec[VSIZE] ;
670
671 if( num <= 0 || vec == NULL || exp0d_pc == NULL ) return ;
672
673 #if 0
674 fprintf(stderr,"Enter EXP0D_worker\n") ;
675 #endif
676
677 if( first ){
678 for( ii=0 ; ii < 26 ; ii++)
679 atoz[ii] = (double *) malloc(sizeof(double) * VSIZE ) ;
680 first = 0 ;
681 #if 0
682 fprintf(stderr,"Allocated atoz\n") ;
683 #endif
684 }
685
686 for( ii=0 ; ii < 26 ; ii++ )
687 for (jj=0; jj<VSIZE; jj++) atoz[ii][jj] = 0.0 ;
688
689 #if 0
690 fprintf(stderr,"Zeroed atoz\n") ;
691 #endif
692
693 for( ii=0 ; ii < num ; ii+=VSIZE ){
694 jbot = ii ;
695 jtop = MIN( ii + VSIZE , num ) ;
696
697 for( jj=jbot ; jj < jtop ; jj ++ ) atoz[exp0d_var][jj-ii] = vec[jj] ;
698
699 PARSER_evaluate_vector( exp0d_pc , atoz , jtop-jbot , tvec ) ;
700
701 for( jj=jbot ; jj < jtop ; jj ++ ) vec[jj] = tvec[jj-ii] ;
702 }
703
704 #if 0
705 fprintf(stderr,"Exit EXP0D_worker\n") ;
706 #endif
707
708 return ;
709 }
710
711 /*------------------------------------------------------------------*/
712
713 #ifdef ALLOW_LOMO
714 static int lomo_order = 3 ;
715 static int lomo_levels = 32 ;
716 int lomo_regress( int , int , int * , int * ) ;
717
LOMOR_init(void)718 PLUGIN_interface * LOMOR_init(void)
719 {
720 PLUGIN_interface * plint ;
721
722 /*---------------- set titles and call point ----------------*/
723
724 plint = PLUTO_new_interface( "Lomo Regression" ,
725 "Control the Local Monotone Regression" ,
726 NULL ,
727 PLUGIN_CALL_VIA_MENU , LOMOR_main ) ;
728
729 PLUTO_add_option( plint , "Parameters" , "Parameters" , TRUE ) ;
730 PLUTO_add_number( plint , "Order" , 3, 20,0,lomo_order , FALSE ) ;
731 PLUTO_add_number( plint , "Levels", 4,128,0,lomo_levels , FALSE ) ;
732
733 PLUTO_register_1D_function( "Lomo 1D" , LOMOR_fitter ) ;
734
735 return plint ;
736 }
737
LOMOR_main(PLUGIN_interface * plint)738 char * LOMOR_main( PLUGIN_interface * plint )
739 {
740 PLUTO_next_option(plint) ;
741 lomo_order = PLUTO_get_number(plint) ;
742 lomo_levels = PLUTO_get_number(plint) ;
743 return NULL ;
744 }
745
746 static int lnum = 0 ;
747 static int * xin = NULL ;
748 static int * yout = NULL ;
749
LOMOR_fitter(int num,double to,double dt,float * vec)750 void LOMOR_fitter( int num , double to , double dt , float * vec )
751 {
752 int ii ;
753 float top , bot , scl ;
754
755 if( num <= 3 || vec == NULL ) return ;
756
757 if( lnum < num ){
758 if( xin != NULL ){ free(xin) ; free(yout) ; }
759 xin = (int *) malloc( sizeof(int) * num ) ;
760 yout = (int *) malloc( sizeof(int) * num ) ;
761 if( xin == NULL || yout == NULL ) return ;
762 lnum = num ;
763 }
764
765 bot = top = vec[0] ;
766 for( ii=1 ; ii < num ; ii++ ){
767 if( vec[ii] < bot ) bot = vec[ii] ;
768 else if( vec[ii] > top ) top = vec[ii] ;
769 }
770
771 if( bot >= top ) return ;
772 scl = (lomo_levels - 0.01) / (top-bot) ;
773
774 fprintf(stderr,"LOMO: range %f .. %f; scl=%f\n",bot,top,scl) ;
775
776 for( ii=0 ; ii < num ; ii++ ) xin[ii] = (int)((vec[ii]-bot)*scl) ;
777
778 ii = lomo_regress( num , lomo_order , xin , yout ) ;
779 if( ii == -1 ) return ;
780
781 scl = 1.0 / scl ;
782 for( ii=0 ; ii < num ; ii++ ) vec[ii] = bot + yout[ii]*scl ;
783 return ;
784 }
785
786 /*******************************************************/
787 /* Fast Digital Locally Monotonic Regression */
788 /*******************************************************/
789 /* Copyright (c) 1995 */
790 /* University of Maryland at College Park */
791 /* All Rights Reserved */
792 /*******************************************************/
793 /* by Nicholas Sidiropoulos, Aug. 1995 */
794 /*******************************************************/
795 /* compile using -lm option */
796 /*******************************************************/
797 /* input: from stdinput.dat: standard matlab vector */
798 /* (i.e., ASCII file containing */
799 /* a long line of N vector elements */
800 /* separated by spaces) */
801 /* output: in stdoutput.dat.*: same format as input */
802 /* control: in stdcontrol_switch.dat: the */
803 /* ``effective'' M is */
804 /* the first entry here, followed by `\n` */
805 /* Rest: fixed to 1 (future option) */
806 /* So stdcontrol_switch .dat should be ASCII */
807 /* file starting with the effective M \n */
808 /* and followed by N `1''s \n, e.g., for M=10 */
809 /* 10<newline>1<newline>...1<newline> */
810 /* a total of N ones */
811 /*******************************************************/
812 /* This is just an ``Academic'' piece of software - */
813 /* it has been checked for correctness to the best */
814 /* of my ability, however, no guarantees whatsoever */
815 /* are given - all disclaimers here! */
816 /* It has been optimized for minimum development effort*/
817 /* and NOT for optimum speed and/or minimum comput. */
818 /* resources. Note, in particular, that I do not take */
819 /* advantage of trellis path merging to reduce storage */
820 /* requirements. As a result, depending on your choice */
821 /* of parameters M,N,R, the program may require */
822 /* considerable amounts of static memory. Therefore, */
823 /* if your computer lacks it you need to rlogin to */
824 /* a machine (like oxygen or apollo at SIL lab) which */
825 /* has sufficient memory */
826 /*******************************************************/
827
828 /*=====================================================*/
829 /*== Modified Feb 1997 by RWCox, to be a C function. ==*/
830 /*=====================================================*/
831
832 #include <stdio.h>
833 #include <math.h>
834 #include <string.h>
835
836 typedef struct _state {
837 int value, length, cumcost;
838 struct _state *prevstate;
839 } state;
840
841 #ifdef ABS
842 #undef ABS
843 #endif
844 #define ABS(x) (((x)<0) ? (-(x)) : (x))
845
846 #define USE_STATIC_TRELLIS
847 #ifdef USE_STATIC_TRELLIS
848 # define MMAX 35 /* (strictly) > max lomo degree */
849 # define NMAX 512 /* input data length */
850 # define RMAX 100 /* range of input: 0 -> R-1 inc. */
851 static state TRELLIS[RMAX][2*MMAX][NMAX] ;
852 # define trellis(i,j,k) TRELLIS[i][j][k]
853 #else
854 static state * TRELLIS = NULL ;
855 static int TDIM1 , TDIM2 , TDIM12 , TDIM3 ;
856 # define trellis(i,j,k) TRELLIS[i+j*TDIM1+k*TDIM12]
857 #endif
858
859 /*=====================================================*/
860 /*== N = number of input points ==*/
861 /*== m = local monotonicity order to impose ==*/
862 /*== yin = pointer to inputs (array of length N) ==*/
863 /*== x = pointer to outputs (array of length N) ==*/
864 /*== ==*/
865 /*== return value is 0 if all is OK, -1 if not ==*/
866 /*=====================================================*/
867
lomo_regress(int N,int m,int * yin,int * x)868 int lomo_regress( int N , int m , int * yin , int * x )
869 {
870 int base , R , * y ;
871 int n,v,l,pv,pl,i,j;
872 int cost, maxcost;
873 int peakval;
874 state dummy_state; /* dummy initial state */
875
876 /*== check inputs for OK-ness ==*/
877
878 if( N < 3 || m < 3 || m >= N || yin == NULL || x == NULL ) return -1 ;
879
880 /*== Compute range of input into R ==*/
881
882 y = (int *) malloc(sizeof(int)*N) ; if( y == NULL ) return -1 ;
883 base = yin[0] ;
884 for( i=1 ; i < N ; i++ ) if( yin[i] < base ) base = yin[i] ;
885 R = y[0] = yin[0] - base ;
886 for( i=1 ; i < N ; i++ ){
887 y[i] = yin[i] - base ; if( y[i] > R ) R = y[i] ;
888 }
889 R++ ;
890 if( R == 1 ){
891 for( i=0 ; i < N ; i++ ) x[i] = yin[i] ;
892 free(y) ; return 0 ;
893 }
894
895 fprintf(stderr,"LOMO: %d points; %d levels; %d order\n",N,R,m) ;
896
897 /* compute maxcost */
898
899 maxcost = 0;
900 for (n = 0; n < N; n++) maxcost += ABS(y[n]);
901
902 /* now init FLOMOR : */
903
904 dummy_state.value = -1; dummy_state.length = m ;
905 dummy_state.cumcost = 0 ; dummy_state.prevstate = NULL;
906
907 #ifndef USE_STATIC_TRELLIS
908 /*== malloc trellis space ==*/
909
910 TDIM1 = R ; TDIM2 = 2*m+2 ; TDIM12 = TDIM1*TDIM2 ; TDIM3 = N ;
911 TRELLIS = (state *) calloc( TDIM1*TDIM2*TDIM3 , sizeof(state) ) ;
912 if( TRELLIS == NULL ){ free(y) ; return -1 ; }
913 #endif
914
915 fprintf(stderr,"LOMO: init trellis\n") ;
916
917 for (n = 0; n < N; n++)
918 {
919 for (l = 1; l <= 2*m; l++)
920 {
921 for (v = 0; v < R; v++) { trellis(v,l,n).value = v;
922 trellis(v,l,n).length = l;
923 if ((n == 0) && ((l == 1) || (l == m+1)))
924 { trellis(v,l,n).cumcost = ABS((v-y[0]));
925 trellis(v,l,n).prevstate = &dummy_state;
926 }
927 else
928 { trellis(v,l,n).cumcost = maxcost;
929 trellis(v,l,n).prevstate = NULL;
930 }
931 }
932 }
933 }
934
935 /************************ main FLOMOR loop ************************************/
936 /* note that l = -1, ..., -m is mapped to l = m+1, ..., 2m resp. */
937 /******************************************************************************/
938
939 fprintf(stderr,"LOMO: compute trellis") ; fflush(stderr) ;
940
941 for (n = 1; n < N; n++)
942 {
943 /* states of the first type: (v,1) */
944
945 fprintf(stderr,".") ; fflush(stderr) ;
946
947 for (v = 0; v < R; v++)
948 {
949 for (pv = 0; pv < v; pv++)
950 {
951 for (pl = 1; pl <= m; pl++)
952 {
953 if (trellis(pv,pl,n-1).cumcost != maxcost)
954 {
955 cost = trellis(pv,pl,n-1).cumcost + ABS((v-y[n]));
956 if (cost < trellis(v,1,n).cumcost)
957 {
958 trellis(v,1,n).cumcost = cost;
959 trellis(v,1,n).prevstate = &trellis(pv,pl,n-1);
960 }
961 }
962 }
963 if (trellis(pv,2*m,n-1).cumcost != maxcost)
964 {
965 cost = trellis(pv,2*m,n-1).cumcost + ABS((v-y[n]));
966 if (cost < trellis(v,1,n).cumcost)
967 {
968 trellis(v,1,n).cumcost = cost;
969 trellis(v,1,n).prevstate = &trellis(pv,2*m,n-1);
970 }
971 }
972 }
973 }
974
975 /* states of the second type: (v,-1) */
976
977 for (v = 0; v < R; v++)
978 {
979 for (pv = R - 1; pv > v; pv--)
980 {
981 for (pl = m + 1; pl <= 2*m; pl++)
982 {
983 if (trellis(pv,pl,n-1).cumcost != maxcost)
984 {
985 cost = trellis(pv,pl,n-1).cumcost + ABS((v-y[n]));
986 if (cost < trellis(v,m+1,n).cumcost)
987 {
988 trellis(v,m+1,n).cumcost = cost;
989 trellis(v,m+1,n).prevstate = &trellis(pv,pl,n-1);
990 }
991 }
992 }
993 if (trellis(pv,m,n-1).cumcost != maxcost)
994 {
995 cost = trellis(pv,m,n-1).cumcost + ABS((v-y[n]));
996 if (cost < trellis(v,m+1,n).cumcost)
997 {
998 trellis(v,m+1,n).cumcost = cost;
999 trellis(v,m+1,n).prevstate = &trellis(pv,m,n-1);
1000 }
1001 }
1002 }
1003 }
1004
1005 /* states of the third type: (v,l), 1 < l < m */
1006
1007 for (l = 2; l < m; l++)
1008 {
1009 for (v = 0; v < R; v++)
1010 {
1011 if (trellis(v,l-1,n-1).cumcost != maxcost)
1012 {
1013 trellis(v,l,n).cumcost = trellis(v,l-1,n-1).cumcost + ABS((v-y[n]));
1014 trellis(v,l,n).prevstate = &trellis(v,l-1,n-1);
1015 }
1016 }
1017 }
1018
1019 /* states of the fourth type: (v,l), -m < l < -1 */
1020
1021 for (l = m+2; l < 2*m; l++)
1022 {
1023 for (v = 0; v < R; v++)
1024 {
1025 if (trellis(v,l-1,n-1).cumcost != maxcost)
1026 {
1027 trellis(v,l,n).cumcost = trellis(v,l-1,n-1).cumcost + ABS((v-y[n]));
1028 trellis(v,l,n).prevstate = &trellis(v,l-1,n-1);
1029 }
1030 }
1031 }
1032
1033 /* states of the fifth type: (v,l), l = m */
1034
1035 for (v = 0; v < R; v++)
1036 {
1037 if (trellis(v,m,n-1).cumcost != maxcost)
1038 {
1039 cost = trellis(v,m,n-1).cumcost + ABS((v-y[n]));
1040 if (cost < trellis(v,m,n).cumcost)
1041 {
1042 trellis(v,m,n).cumcost = cost;
1043 trellis(v,m,n).prevstate = &trellis(v,m,n-1);
1044 }
1045 }
1046
1047 if (trellis(v,m-1,n-1).cumcost != maxcost)
1048 {
1049 cost = trellis(v,m-1,n-1).cumcost + ABS((v-y[n]));
1050 if (cost < trellis(v,m,n).cumcost)
1051 {
1052 trellis(v,m,n).cumcost = cost;
1053 trellis(v,m,n).prevstate = &trellis(v,m-1,n-1);
1054 }
1055 }
1056 }
1057
1058 /* states of the sixth type: (v,l), l = -m */
1059
1060 for (v = 0; v < R; v++)
1061 {
1062 if (trellis(v,2*m,n-1).cumcost != maxcost)
1063 {
1064 cost = trellis(v,2*m,n-1).cumcost + ABS((v-y[n]));
1065 if (cost < trellis(v,2*m,n).cumcost)
1066 {
1067 trellis(v,2*m,n).cumcost = cost;
1068 trellis(v,2*m,n).prevstate = &trellis(v,2*m,n-1);
1069 }
1070 }
1071
1072 if (trellis(v,2*m-1,n-1).cumcost != maxcost)
1073 {
1074 cost = trellis(v,2*m-1,n-1).cumcost + ABS((v-y[n]));
1075 if (cost < trellis(v,2*m,n).cumcost)
1076 {
1077 trellis(v,2*m,n).cumcost = cost;
1078 trellis(v,2*m,n).prevstate = &trellis(v,2*m-1,n-1);
1079 }
1080 }
1081 }
1082
1083
1084 }
1085
1086 /************************ eof main FLOMOR loop *******************************/
1087
1088 /* now pick best path, and write it to x[n] */
1089
1090 fprintf(stderr,"\nLOMO: pick best path\n") ; fflush(stderr) ;
1091
1092 cost = maxcost;
1093 for (l = 1; l <= m; l++)
1094 {
1095 for (v = 0; v < R; v++)
1096 {
1097 if (trellis(v,l,N-1).cumcost < cost)
1098 {
1099 cost = trellis(v,l,N-1).cumcost;
1100 pl = l; pv = v;
1101 }
1102 }
1103 }
1104
1105 /* now traverse backwards */
1106
1107 if (cost < maxcost)
1108 {
1109 x[N-1] = pv;
1110 fprintf(stderr," [%d] = %d\n",N-1,pv) ;
1111
1112 for (n = N-2; n >=0; n--)
1113 {
1114 fprintf(stderr," [%d] = %d",n,trellis(pv,pl,n+1).prevstate->value) ; fflush(stderr) ;
1115 x[n] = trellis(pv,pl,n+1).prevstate->value;
1116 fprintf(stderr," pv = %d",trellis(pv,pl,n+1).prevstate->value) ; fflush(stderr) ;
1117 pv = trellis(pv,pl,n+1).prevstate->value;
1118 fprintf(stderr," pl = %d",trellis(pv,pl,n+1).prevstate->length) ; fflush(stderr) ;
1119 pl = trellis(pv,pl,n+1).prevstate->length;
1120 fprintf(stderr,"\n") ;
1121 }
1122 }
1123 else { for (n=0;n<N;n++) x[n] = 0; /* as good as any... */ }
1124
1125 /*== add base back to output ==*/
1126
1127 fprintf(stderr,"\nLOMO: add base back on\n") ;
1128
1129 for( i=0 ; i < N ; i++ ) x[n] += base ;
1130
1131 fprintf(stderr,"LOMO: exit regression\n") ;
1132
1133 free(y) ;
1134 #ifndef USE_STATIC_TRELLIS
1135 free(TRELLIS) ;
1136 #endif
1137
1138 return 0 ;
1139 }
1140 #endif /* ALLOW_LOMO */
1141