00602 {
00603    PLUG_startup_plugin_CB( NULL , (XtPointer)plint_EXP0D , NULL ) ;
00604 }

00608 {
00609    PLUGIN_interface * plint ;
00610 
00611    /*---------------- set titles and call point ----------------*/
00612 
00613    plint = PLUTO_new_interface( "Expr 0D" ,
00614                                 "Control the Expr 0D transformation" ,
00615                                 fredstring ,
00616                                 PLUGIN_CALL_VIA_MENU , EXP0D_main ) ;
00617 
00618    PLUTO_add_option( plint , "Variable" , "Variable" , TRUE ) ;
00619    PLUTO_add_string( plint , NULL , NALPHA,vstring , exp0d_var ) ;
00620 
00621    PLUTO_add_option( plint , "Expression" , "Expression" , TRUE ) ;
00622    PLUTO_add_string( plint , NULL , 0,NULL , 50 ) ;
00623 
00624    PLUTO_register_0D_function( "Expr 0D" , EXP0D_worker ) ;
00625 
00626    plint_EXP0D = plint ;
00627    AFNI_register_nD_func_init( 0 , (generic_func *)EXP0D_func_init ) ;  /* 21 Jul 2003 */
00628 
00629    return plint ;
00630 }

00633 {
00634    char * str ;
00635 
00636    PLUTO_next_option(plint) ;
00637    str       = PLUTO_get_string(plint) ;
00638    exp0d_var = PLUTO_string_index( str , NALPHA , vstring ) ;
00639 
00640    if( exp0d_pc != NULL ){ free(exp0d_pc) ; exp0d_pc = NULL ; }
00641    PLUTO_next_option(plint) ;
00642    str       = PLUTO_get_string(plint) ;
00643    exp0d_pc  = PARSER_generate_code( str ) ;
00644 
00645    if( exp0d_pc == NULL ) return "*******************************\n"
00646                                  "Error when compiling expression\n"
00647                                  "*******************************"  ;
00648 
00649    return NULL ;
00650 }

00655 {
00656    int ii , jj , jbot,jtop ;
00657 
00658    static int first = 1 ;
00659    static double * atoz[26] ;
00660    static double   tvec[VSIZE] ;
00661 
00662    if( num <= 0 || vec == NULL || exp0d_pc == NULL ) return ;
00663 
00664 #if 0
00665 fprintf(stderr,"Enter EXP0D_worker\n") ;
00666 #endif
00667 
00668    if( first ){
00669       for( ii=0 ; ii < 26 ; ii++)
00670         atoz[ii] = (double *) malloc(sizeof(double) * VSIZE ) ;
00671       first = 0 ;
00672 #if 0
00673 fprintf(stderr,"Allocated atoz\n") ;
00674 #endif
00675    }
00676 
00677    for( ii=0 ; ii < 26 ; ii++ )
00678       for (jj=0; jj<VSIZE; jj++) atoz[ii][jj] = 0.0 ;
00679 
00680 #if 0
00681 fprintf(stderr,"Zeroed atoz\n") ;
00682 #endif
00683 
00684    for( ii=0 ; ii < num ; ii+=VSIZE ){
00685       jbot = ii ;
00686       jtop = MIN( ii + VSIZE , num ) ;
00687 
00688       for( jj=jbot ; jj < jtop ; jj ++ ) atoz[exp0d_var][jj-ii] = vec[jj] ;
00689 
00690       PARSER_evaluate_vector( exp0d_pc , atoz , jtop-jbot , tvec ) ;
00691 
00692       for( jj=jbot ; jj < jtop ; jj ++ ) vec[jj] = tvec[jj-ii] ;
00693    }
00694 
00695 #if 0
00696 fprintf(stderr,"Exit EXP0D_worker\n") ;
00697 #endif
00698 
00699    return ;
00700 }

00256 {
00257    LSQ_worker( nt , dt , vec , FALSE , label ) ;
00258    return ;
00259 }

00250 {
00251    LSQ_worker( nt , dt , vec , TRUE , label ) ;
00252    return ;
00253 }

00175 {
00176    char * str ;
00177    int  ii ;
00178    float * tsar ;
00179 
00180    /*--------- go to first input line ---------*/
00181 
00182    PLUTO_next_option(plint) ;
00183 
00184    str    = PLUTO_get_string(plint) ;
00185    polort = PLUTO_string_index( str , NBASE , baseline_strings ) ;
00186 
00187    ignore = PLUTO_get_number(plint) ;
00188 
00189    /*------ loop over remaining options, check their tags, process them -----*/
00190 
00191    nrsin = nrts = 0 ;
00192    do {
00193       str = PLUTO_get_optiontag(plint) ; if( str == NULL ) break ;
00194 
00195       if( strcmp(str,"Sinusoid") == 0 ){
00196 
00197          sinper[nrsin]  = PLUTO_get_number(plint) ;
00198          sinharm[nrsin] = PLUTO_get_number(plint) - 1.0 ;
00199          if( sinper[nrsin] <= 0.0 )
00200             return "************************\n"
00201                    "Illegal Sinusoid Period!\n"
00202                    "************************"  ;
00203 
00204          nrsin++ ;
00205 
00206       } else if( strcmp(str,"Timeseries") == 0 ){
00207 
00208          tsim[nrts] = PLUTO_get_timeseries(plint) ;
00209 
00210          if( tsim[nrts] == NULL || tsim[nrts]->nx < 3 || tsim[nrts]->kind != MRI_float )
00211             return "*************************\n"
00212                    "Illegal Timeseries Input!\n"
00213                    "*************************"  ;
00214 
00215          tsar = MRI_FLOAT_PTR(tsim[nrts]) ;
00216          for( ii=ignore ; ii < tsim[nrts]->nx && tsar[ii] >= WAY_BIG ; ii++ ) ; /* nada */
00217          ignore = ii ;
00218          nrts++ ;
00219 
00220       } else {
00221          return "************************\n"
00222                 "Illegal optiontag found!\n"
00223                 "************************"  ;
00224       }
00225    } while(1) ;
00226 
00227    /*--- nothing left to do until data arrives ---*/
00228 
00229    initialize = 1 ;  /* force re-initialization */
00230 
00231    /*** compute how many ref functions are ordered ***/
00232 
00233    { int nref , ks ;
00234      char str[64] ;
00235 
00236      nref = (polort+1) + nrts ;
00237      for( ks=0 ; ks < nrsin ; ks++ ) nref += 2*(sinharm[ks]+1) ;
00238      sprintf(str," \nNumber of fit parameters = %d\n",nref) ;
00239      PLUTO_popup_transient( plint , str ) ;
00240    }
00241 
00242    return NULL ;
00243 }

00265 {
00266    int nlen , nref ;
00267 
00268    static int nlen_old = -666 , nref_old = -666 ;
00269    static double dt_old = -666.666 ;
00270    static float ** ref = NULL ;
00271    static double * dch = NULL ;
00272 
00273    int ir , ii , ks,jh , j;
00274    float fac , tm , val ;
00275    float * fit , * tsar ;
00276 
00277    /*** compute how many ref functions are ordered ***/
00278 
00279     nref = (polort+1) + nrts ;
00280     for( ks=0 ; ks < nrsin ; ks++ ) nref += 2*(sinharm[ks]+1) ;
00281 
00282    /*** do nothing if not enough data to fit ***/
00283 
00284    nlen = nt - ignore ;
00285 
00286    if( nlen <= nref ) return ;  /* do nothing if not enough data to fit */
00287 
00288    /** if data vectors are new length,
00289        or have a new number of reference vectors,
00290        or have a new time step and need sinusoids,
00291        or the initialize flag is set,
00292        then reinitialize reference vectors and Choleski factor **/
00293 
00294    if( nlen != nlen_old || nref != nref_old ||
00295        initialize       || (dt != dt_old && nrsin > 0) ){
00296 
00297       /* free old storage */
00298 
00299       if( ref != NULL ){
00300          for( ir=0 ; ir < nref_old ; ir++ ) if( ref[ir] != NULL ) free(ref[ir]) ;
00301          free(ref) ;
00302       }
00303       if( dch != NULL ) free(dch) ;
00304 
00305       /* make space for ref vectors */
00306 
00307       ref = (float **) malloc( sizeof(float *) * nref ) ;
00308       if( ref == NULL ){fprintf(stderr,"\nmalloc error in plug_lsqfit\n\a");
00309          return;
00310          /* EXIT(1); */
00311       }
00312       for( ir=0 ; ir < nref ; ir++ ){
00313          ref[ir] = (float *) malloc( sizeof(float) * nlen ) ;
00314          if( ref[ir] == NULL ){
00315             fprintf(stderr,"\nmalloc error in plug_lsqfit\n\a");
00316             for(j=0;j<=ir;j++)
00317               free(ref[j]);
00318             free(ref);
00319             ref = NULL;
00320             return;
00321             /* EXIT(1); */
00322          }
00323       }
00324       nlen_old = nlen ;
00325       nref_old = nref ;
00326       dt_old   = dt ;
00327 
00328       /**** fill ref vectors ****/
00329 
00330       /* r(t) = 1 */
00331 
00332       for( ii=0 ; ii < nlen ; ii++ ) ref[0][ii] = 1.0 ;
00333 
00334       ir = 1 ;
00335       if( polort > 0 ){
00336 
00337          /* r(t) = t - tmid */
00338 
00339          tm = 0.5 * (nlen-1.0) ; fac = 2.0 / nlen ;
00340          for( ii=0 ; ii < nlen ; ii++ ) ref[1][ii] = (ii-tm)*fac ;
00341          ir = 2 ;
00342 
00343          /* r(t) = (t-tmid)**ir */
00344 
00345          for( ; ir <= polort ; ir++ )
00346             for( ii=0 ; ii < nlen ; ii++ )
00347                ref[ir][ii] = pow( (ii-tm)*fac , (double)ir ) ;
00348       }
00349 
00350       if( dt == 0.0 ) dt = 1.0 ;
00351 
00352       /* r(t) = sinusoids */
00353 
00354       for( ks=0 ; ks < nrsin ; ks++ ){
00355          for( jh=0 ; jh <= sinharm[ks] ; jh++ ){
00356             fac = (2.0*PI) * dt * (jh+1) / sinper[ks] ;
00357             for( ii=0 ; ii < nlen ; ii++ ){
00358                ref[ir]  [ii] = cos( fac * ii ) ;
00359                ref[ir+1][ii] = sin( fac * ii ) ;
00360             }
00361             ir += 2 ;
00362          }
00363       }
00364 
00365       /* r(t) = timeseries files */
00366 
00367       for( ks=0 ; ks < nrts ; ks++ ){
00368          if( tsim[ks] == NULL || tsim[ks]->nx - ignore < nlen ){
00369             initialize = 1 ;
00370             fprintf(stderr,"Inadequate time series #%d in LSQ plugin\n\a",ks+1) ;
00371             return ;
00372          }
00373          tsar = MRI_FLOAT_PTR(tsim[ks]) ;
00374          for( ii=0 ; ii < nlen ; ii++ ) ref[ir][ii] = tsar[ii+ignore] ;
00375          ir++ ;
00376       }
00377 
00378       /* Cholesky-ize */
00379 
00380       dch = startup_lsqfit( nlen , NULL , nref , ref ) ;
00381       if( dch == NULL ){
00382          initialize = 1 ;
00383          fprintf(stderr,"Choleski error in LSQ plugin\n\a") ;
00384          return ;
00385       }
00386 
00387       initialize = 0 ;
00388    }
00389 
00390    /** find least squares fit coefficients **/
00391 
00392    fit = delayed_lsqfit( nlen , vec+ignore , nref , ref , dch ) ;
00393 
00394    for( ii=0 ; ii < nlen ; ii++ ){
00395       val = 0.0 ;
00396       for( ir=0 ; ir < nref ; ir++ ) val += fit[ir] * ref[ir][ii] ;
00397 
00398       vec[ii+ignore] = (dofit) ? val : vec[ii+ignore] - val ;
00399    }
00400 
00401    /** 22 Apr 1997: create label if desired by AFNI         **/
00402    /** [This is in static storage, since AFNI will copy it] **/
00403 
00404    if( label != NULL ){  /* assemble this 1 line at a time from sbuf */
00405 
00406       lbuf[0] = '\0' ;   /* make this a 0 length string to start */
00407 
00408       /** for each reference, make a string into sbuf **/
00409 
00410       ir = 0 ;
00411       sprintf(sbuf,"Coef of 1 = %g\n" , fit[ir++] ) ;
00412       strcat(lbuf,sbuf) ;
00413 
00414       for( ; ir <= polort ; ){
00415          sprintf(sbuf,"Coef of t**%d = %g\n" , ir,fit[ir++] ) ;
00416          strcat(lbuf,sbuf) ;
00417       }
00418 
00419       for( ks=0 ; ks < nrsin ; ks++ ){
00420          for( jh=0 ; jh <= sinharm[ks] ; jh++ ){
00421             fac = sinper[ks] / (jh+1) ;
00422             sprintf(sbuf,"Coef of cos(2*Pi*t/%g) = %g\n" , fac , fit[ir++] ) ;
00423             strcat(lbuf,sbuf) ;
00424             sprintf(sbuf,"Coef of sin(2*Pi*t/%g) = %g\n" , fac , fit[ir++] ) ;
00425             strcat(lbuf,sbuf) ;
00426          }
00427       }
00428 
00429       for( ks=0 ; ks < nrts ; ks++ ){
00430          sprintf(sbuf,"Coef of %s = %g\n" , tsim[ks]->name , fit[ir++] ) ;
00431          strcat(lbuf,sbuf) ;
00432       }
00433 
00434       *label = lbuf ;  /* send address of lbuf back in what label points to */
00435    }
00436 
00437    free(fit) ;
00438    return ;
00439 }

00106 {
00107    int ii ;
00108    PLUGIN_interface * plint ;     /* will be the output of this routine */
00109 
00110    if( ncall > 3 ) return NULL ;  /* generate interfaces for ncall 0-3 */
00111 
00112    if( ncall == 1 ) return TSGEN_init() ;  /* interface # 1 */
00113    if( ncall == 2 ) return EXP0D_init() ;  /* interface # 2 */
00114 #ifdef ALLOW_LOMO
00115    if( ncall == 3 ) return LOMOR_init() ;  /* interface # 3 */
00116 #else
00117    if( ncall == 3 ) return NULL ;
00118 #endif
00119 
00120    /***** otherwise, do interface # 0 *****/
00121 
00122    /*---------------- set titles and call point ----------------*/
00123 
00124    plint = PLUTO_new_interface( "LSqFit & Dtr" ,
00125                                 "Control LSqFit and LSqDtr Functions" ,
00126                                 helpstring ,
00127                                 PLUGIN_CALL_VIA_MENU , LSQ_main ) ;
00128 
00129    global_plint = plint ;  /* make global copy */
00130 
00131    PLUTO_set_sequence( plint , "A:funcs:fitting" ) ;
00132 
00133    PLUTO_add_hint( plint , "Control LSqFit and LSqDtr Functions" ) ;
00134 
00135    PLUTO_set_runlabels( plint , "Set+Keep" , "Set+Close" ) ;  /* 04 Nov 2003 */
00136 
00137    /*----- Parameters -----*/
00138 
00139    PLUTO_add_option( plint , "Parameters" , "Parameters" , TRUE ) ;
00140 
00141    PLUTO_add_string( plint , "Baseline" , NBASE , baseline_strings , 1 ) ;
00142 
00143    PLUTO_add_number( plint , "Ignore" , 0,20,0,3 , FALSE ) ;
00144 
00145    /*----- Sinusoid -----*/
00146 
00147    for( ii=0 ; ii < NRMAX_SIN ; ii++ ){
00148       PLUTO_add_option( plint , "Sinusoid" , "Sinusoid" , FALSE ) ;
00149       PLUTO_add_number( plint , "Period" , 0,99999,0,20, TRUE ) ;
00150       PLUTO_add_number( plint , "Harmonics" , 1,HARM_MAX,0,1 , FALSE ) ;
00151    }
00152 
00153    /*----- Timeseries -----*/
00154 
00155    for( ii=0 ; ii < NRMAX_TS ; ii++ ){
00156       PLUTO_add_option( plint , "Timeseries" , "Timeseries" , FALSE ) ;
00157       PLUTO_add_timeseries( plint , "File" ) ;
00158    }
00159 
00160    /*--------- done with interface setup ---------*/
00161 
00162    PLUTO_register_1D_funcstr( "LSqFit" , LSQ_fitter ) ;
00163    PLUTO_register_1D_funcstr( "LSqDtr" , LSQ_detrend ) ;
00164 
00165    return plint ;
00166 }

00444 {
00445    PLUGIN_interface * plint ;
00446 
00447    /*---------------- set titles and call point ----------------*/
00448 
00449    plint = PLUTO_new_interface( "TS Generate" ,
00450                                 "Generate a Timeseries" ,
00451                                 plehstring ,
00452                                 PLUGIN_CALL_VIA_MENU , TSGEN_main ) ;
00453 
00454    PLUTO_add_hint( plint , "Generate a 1D Timeseries" ) ;
00455 
00456    /*----- Parameters -----*/
00457 
00458    PLUTO_add_option( plint , "Parameters" , "Parameters" , TRUE ) ;
00459    PLUTO_add_number( plint , "Delta"  , 0,99999,1, 0 , TRUE ) ;
00460    PLUTO_add_number( plint , "Length" , 3,9999,0,3   , TRUE ) ;
00461 
00462    /*----- Output -----*/
00463 
00464    PLUTO_add_option( plint , "Output" , "Output" , TRUE ) ;
00465    PLUTO_add_string( plint , "Label" , 0,NULL , 19 ) ;
00466 
00467    /*----- Polynomial -----*/
00468 
00469    PLUTO_add_option( plint , "Polynomial" , "Polynomial" , FALSE ) ;
00470    PLUTO_add_number( plint , "Order" , 2,20,0,2 , FALSE ) ;
00471 
00472    /*----- Sinusoid -----*/
00473 
00474    PLUTO_add_option( plint , "Sinusoid" , "Sinusoid" , FALSE ) ;
00475    PLUTO_add_number( plint , "Period" , 0,99999,0,20, TRUE ) ;
00476    PLUTO_add_number( plint , "Harmonics" , 1,HARM_MAX,0,1 , FALSE ) ;
00477 
00478    return plint ;
00479 }

00482 {
00483    char * label , * str ;
00484    int  ii , jj ;
00485    float * tsar ;
00486    MRI_IMAGE * tsim ;
00487    float delta , period=0.0 ;
00488    int   nx , ny=0 , npol=0 , nharm=-1 ;
00489    int pp ;
00490    double fac , val ;
00491 
00492    /*--------- go to first input line ---------*/
00493 
00494    PLUTO_next_option(plint) ;
00495 
00496    delta = PLUTO_get_number(plint) ;
00497    if( delta <= 0.0 ) return "**********************\n"
00498                              "Illegal value of Delta\n"
00499                              "**********************"  ;
00500 
00501    nx = PLUTO_get_number(plint) ;
00502 
00503    /*----- next input line -----*/
00504 
00505    PLUTO_next_option(plint) ;
00506    label = PLUTO_get_string(plint) ;
00507    if( label == NULL || strlen(label) == 0 ) return "**********************\n"
00508                                                     "Illegal value of Label\n"
00509                                                     "**********************"  ;
00510 
00511    /*----- rest of input lines -----*/
00512 
00513    do {
00514       str = PLUTO_get_optiontag(plint) ; if( str == NULL ) break ;
00515 
00516       if( strcmp(str,"Sinusoid") == 0 ){
00517 
00518          period = PLUTO_get_number(plint) ;
00519          nharm  = PLUTO_get_number(plint) - 1.0 ;
00520 
00521          if( period <= 0.0 ) return "***********************\n"
00522                                     "Illegal Sinusoid Period\n"
00523                                     "***********************"  ;
00524 
00525 
00526       } else if( strcmp(str,"Polynomial") == 0 ){
00527 
00528          npol = PLUTO_get_number(plint) ;
00529 
00530       } else {
00531          return "***********************\n"
00532                 "Illegal optiontag found\n"
00533                 "***********************"  ;
00534       }
00535    } while(1) ;
00536 
00537    /********** Make the timeseries ***********/
00538 
00539    ny = 0 ;
00540    if( npol > 0 )     ny  = npol-1 ;
00541    if( period > 0.0 ) ny += 2*(nharm+1) ;
00542 
00543    if( ny < 1 ) return "***********************\n"
00544                        "No timeseries specified\n"
00545                        "***********************"  ;
00546 
00547    tsim = mri_new( nx , ny , MRI_float ) ;
00548    jj   = 0 ;
00549 
00550    fac  = 1.99999 / (nx-1) ;
00551    for( pp=2 ; pp <= npol ; pp++,jj++ ){
00552 
00553       tsar = MRI_FLOAT_PTR(tsim) + (jj*nx) ;
00554 
00555       for( ii=0 ; ii < nx ; ii++ ){
00556          val = fac * ii - 0.999995 ;
00557          tsar[ii] = cos( pp * acos(val) ) ;
00558       }
00559    }
00560 
00561    for( pp=0 ; pp <= nharm ; pp++ , jj+=2 ){
00562       fac  = (2.0*PI) * delta * (pp+1) / period ;
00563       tsar = MRI_FLOAT_PTR(tsim) + (jj*nx) ;
00564 
00565       for( ii=0 ; ii < nx ; ii++ ){
00566          tsar[ii]    = cos( fac * ii ) ;
00567          tsar[ii+nx] = sin( fac * ii ) ;
00568       }
00569    }
00570 
00571    PLUTO_register_timeseries( label , tsim ) ;
00572    mri_free(tsim) ;
00573    return NULL ;
00574 }

  " Purpose: Control the Expr 0D Transformation function\n"
  "\n"
  " Variable   = letter used as input to expression\n"
  " Expression = arithmetic expression to evaluate\n"

   " Purpose: Control the 'LSqFit' and 'LSqDtr 1D functions.\n"
   "\n"
   " Parameters:  Baseline = 'Constant' or 'Linear'\n"
   "                           Is the baseline 'a' or 'a+b*t'?\n"
   "              Ignore   = Number of points to ignore at\n"
   "                           start of each timeseries.\n"
   " \n"
   " Sinusoids:   Period    = Fundamental period to use.\n"
   "              Harmonics = Number of overtones to use.\n"
   " \n"
   " Timeseries:  File      = Input timeseries file to use.\n"

   " Purpose: Generate a timeseries and store in AFNI list\n"
   "\n"
   " Parameters:  Delta  = time step between points\n"
   "              Length = number of points\n"
   "\n"
   " Output:      Label  = String to label timeseries by\n"
   "                        in AFNI choosers\n"
   "\n"
   " Polynomial:  Order  = Maximum power to include\n"
   "                       (Chebyshev polynomials are used)\n"
   "\n"
   " Sinusoid:    Period    = Fundamental period to use.\n"
   "              Harmonics = Number of overtones to use.\n"

 {
  " A ",  " B ",  " C ",  " D ",  " E ",
  " F ",  " G ",  " H ",  " I ",  " J ",
  " K ",  " L ",  " M ",  " N ",  " O ",
  " P ",  " Q ",  " R ",  " S ",  " T ",
  " U ",  " V ",  " W ",  " X ",  " Y ",  " Z " }