00211 {
00212    L1F_worker( nt , dt , vec , FALSE , label ) ;
00213    return ;
00214 }

00205 {
00206    L1F_worker( nt , dt , vec , TRUE , label ) ;
00207    return ;
00208 }

00130 {
00131    char * str ;
00132    int  ii ;
00133    float * tsar ;
00134 
00135    /*--------- go to first input line ---------*/
00136 
00137    PLUTO_next_option(plint) ;
00138 
00139    str    = PLUTO_get_string(plint) ;
00140    polort = PLUTO_string_index( str , NBASE , baseline_strings ) ;
00141 
00142    ignore = PLUTO_get_number(plint) ;
00143 
00144    /*------ loop over remaining options, check their tags, process them -----*/
00145 
00146    nrsin = nrts = 0 ;
00147    do {
00148       str = PLUTO_get_optiontag(plint) ; if( str == NULL ) break ;
00149 
00150       if( strcmp(str,"Sinusoid") == 0 ){
00151 
00152          sinper[nrsin]  = PLUTO_get_number(plint) ;
00153          sinharm[nrsin] = PLUTO_get_number(plint) - 1.0 ;
00154          if( sinper[nrsin] <= 0.0 )
00155             return "************************\n"
00156                    "Illegal Sinusoid Period!\n"
00157                    "************************"  ;
00158 
00159          nrsin++ ;
00160 
00161       } else if( strcmp(str,"Timeseries") == 0 ){
00162 
00163          tsim[nrts] = PLUTO_get_timeseries(plint) ;
00164 
00165          if( tsim[nrts] == NULL || tsim[nrts]->nx < 3 || tsim[nrts]->kind != MRI_float )
00166             return "*************************\n"
00167                    "Illegal Timeseries Input!\n"
00168                    "*************************"  ;
00169 
00170          tsar = MRI_FLOAT_PTR(tsim[nrts]) ;
00171          for( ii=ignore ; ii < tsim[nrts]->nx && tsar[ii] >= WAY_BIG ; ii++ ) ; /* nada */
00172          ignore = ii ;
00173          nrts++ ;
00174 
00175       } else {
00176          return "************************\n"
00177                 "Illegal optiontag found!\n"
00178                 "************************"  ;
00179       }
00180    } while(1) ;
00181 
00182    /*--- nothing left to do until data arrives ---*/
00183 
00184    initialize = 1 ;  /* force re-initialization */
00185 
00186    /*** compute how many ref functions are ordered ***/
00187 
00188    { int nref , ks ;
00189      char str[64] ;
00190 
00191      nref = (polort+1) + nrts ;
00192      for( ks=0 ; ks < nrsin ; ks++ ) nref += 2*(sinharm[ks]+1) ;
00193      sprintf(str," \nNumber of fit parameters = %d\n",nref) ;
00194      PLUTO_popup_transient( plint , str ) ;
00195    }
00196 
00197    return NULL ;
00198 }

00220 {
00221    int nlen , nref ;
00222 
00223    static int nlen_old = -666 , nref_old = -666 ;
00224    static double dt_old = -666.666 ;
00225    static float ** ref = NULL ;
00226    static float *  fit = NULL ;
00227 
00228    int ir , ii , ks,jh ;
00229    float fac , tm , val , cls ;
00230    float * tsar ;
00231 
00232    /*** compute how many ref functions are ordered ***/
00233 
00234    nref = (polort+1) + nrts ;
00235    for( ks=0 ; ks < nrsin ; ks++ ) nref += 2*(sinharm[ks]+1) ;
00236 
00237    /*** do nothing if not enough data to fit ***/
00238 
00239    nlen = nt - ignore ;
00240 
00241    if( nlen <= nref ) return ;  /* do nothing if not enough data to fit */
00242 
00243    /** if data vectors are new length,
00244        or have a new number of reference vectors,
00245        or have a new time step and need sinusoids,
00246        or the initialize flag is set,
00247        then reinitialize reference vectors and Choleski factor **/
00248 
00249    if( nlen != nlen_old || nref != nref_old ||
00250        initialize       || (dt != dt_old && nrsin > 0) ){
00251 
00252       /* free old storage */
00253 
00254       if( ref != NULL ){
00255          for( ir=0 ; ir < nref_old ; ir++ ) if( ref[ir] != NULL ) free(ref[ir]) ;
00256          free(ref) ;
00257       }
00258       if( fit != NULL ) free(fit) ;
00259 
00260       /* make space for ref vectors */
00261 
00262       ref = (float **) malloc( sizeof(float *) * nref ) ;
00263       if( ref == NULL ){fprintf(stderr,"\nmalloc error in plug_lsqfit\n\a");
00264          return;
00265          /*EXIT(1);*/
00266       }
00267       for( ir=0 ; ir < nref ; ir++ ){
00268          ref[ir] = (float *) malloc( sizeof(float) * nlen ) ;
00269          if( ref[ir] == NULL )
00270             {fprintf(stderr,"\nmalloc error in plug_lsqfit\n\a");
00271             return;
00272             /* EXIT(1); */
00273          }
00274       }
00275       nlen_old = nlen ;
00276       nref_old = nref ;
00277       dt_old   = dt ;
00278 
00279       /**** fill ref vectors ****/
00280 
00281       /* r(t) = 1 */
00282 
00283       for( ii=0 ; ii < nlen ; ii++ ) ref[0][ii] = 1.0 ;
00284 
00285       ir = 1 ;
00286       if( polort > 0 ){
00287 
00288          /* r(t) = t - tmid */
00289 
00290          tm = 0.5 * (nlen-1.0) ; fac = 2.0 / nlen ;
00291          for( ii=0 ; ii < nlen ; ii++ ) ref[1][ii] = (ii-tm)*fac ;
00292          ir = 2 ;
00293 
00294          /* r(t) = (t-tmid)**ir */
00295 
00296          for( ; ir <= polort ; ir++ )
00297             for( ii=0 ; ii < nlen ; ii++ )
00298                ref[ir][ii] = pow( (ii-tm)*fac , (double)ir ) ;
00299       }
00300 
00301       if( dt == 0.0 ) dt = 1.0 ;
00302 
00303       /* r(t) = sinusoids */
00304 
00305       for( ks=0 ; ks < nrsin ; ks++ ){
00306          for( jh=0 ; jh <= sinharm[ks] ; jh++ ){
00307             fac = (2.0*PI) * dt * (jh+1) / sinper[ks] ;
00308             for( ii=0 ; ii < nlen ; ii++ ){
00309                ref[ir]  [ii] = cos( fac * ii ) ;
00310                ref[ir+1][ii] = sin( fac * ii ) ;
00311             }
00312             ir += 2 ;
00313          }
00314       }
00315 
00316       /* r(t) = timeseries files */
00317 
00318       for( ks=0 ; ks < nrts ; ks++ ){
00319          if( tsim[ks] == NULL || tsim[ks]->nx - ignore < nlen ){
00320             initialize = 1 ;
00321             fprintf(stderr,"Inadequate time series #%d in L1F plugin\n\a",ks+1) ;
00322             return ;
00323          }
00324          tsar = MRI_FLOAT_PTR(tsim[ks]) ;
00325          for( ii=0 ; ii < nlen ; ii++ ) ref[ir][ii] = tsar[ii+ignore] ;
00326          ir++ ;
00327       }
00328 
00329       /* make space for fit vector */
00330 
00331       fit = (float *) malloc(sizeof(float)*nref) ;
00332 
00333       initialize = 0 ;
00334    }
00335 
00336    /** find L1 fit coefficients **/
00337 
00338    cls = cl1_solve( nlen , nref , vec+ignore , ref , fit,0 ) ;
00339 
00340    if( cls < 0.0 ) return ;  /* bad fit */
00341 
00342    for( ii=0 ; ii < nlen ; ii++ ){
00343       val = 0.0 ;
00344       for( ir=0 ; ir < nref ; ir++ ) val += fit[ir] * ref[ir][ii] ;
00345 
00346       vec[ii+ignore] = (dofit) ? val : vec[ii+ignore] - val ;
00347    }
00348 
00349    /** 22 Apr 1997: create label if desired by AFNI         **/
00350    /** [This is in static storage, since AFNI will copy it] **/
00351 
00352    if( label != NULL ){  /* assemble this 1 line at a time from sbuf */
00353 
00354       lbuf[0] = '\0' ;   /* make this a 0 length string to start */
00355 
00356       /** for each reference, make a string into sbuf **/
00357 
00358       ir = 0 ;
00359       sprintf(sbuf,"Coef of 1 = %g\n" , fit[ir++] ) ;
00360       strcat(lbuf,sbuf) ;
00361 
00362       for( ; ir <= polort ; ){
00363          sprintf(sbuf,"Coef of t**%d = %g\n" , ir,fit[ir++] ) ;
00364          strcat(lbuf,sbuf) ;
00365       }
00366 
00367       for( ks=0 ; ks < nrsin ; ks++ ){
00368          for( jh=0 ; jh <= sinharm[ks] ; jh++ ){
00369             fac = sinper[ks] / (jh+1) ;
00370             sprintf(sbuf,"Coef of cos(2*Pi*t/%g) = %g\n" , fac , fit[ir++] ) ;
00371             strcat(lbuf,sbuf) ;
00372             sprintf(sbuf,"Coef of sin(2*Pi*t/%g) = %g\n" , fac , fit[ir++] ) ;
00373             strcat(lbuf,sbuf) ;
00374          }
00375       }
00376 
00377       for( ks=0 ; ks < nrts ; ks++ ){
00378          sprintf(sbuf,"Coef of %s = %g\n" , tsim[ks]->name , fit[ir++] ) ;
00379          strcat(lbuf,sbuf) ;
00380       }
00381 
00382       *label = lbuf ;  /* send address of lbuf back in what label points to */
00383    }
00384 
00385    return ;
00386 }

00071 {
00072    int ii ;
00073    PLUGIN_interface * plint ;     /* will be the output of this routine */
00074 
00075    if( ncall > 0 ) return NULL ;
00076 
00077    /*---------------- set titles and call point ----------------*/
00078 
00079    plint = PLUTO_new_interface( "L1_Fit & Dtr" ,
00080                                 "Control L1_Fit and L1_Dtr Functions" ,
00081                                 helpstring ,
00082                                 PLUGIN_CALL_VIA_MENU , L1F_main ) ;
00083 
00084    global_plint = plint ;  /* make global copy */
00085 
00086    PLUTO_set_sequence( plint , "A:funcs:fitting" ) ;
00087 
00088    PLUTO_add_hint( plint , "Control L1_Fit and L1_Dtr Functions" ) ;
00089 
00090    PLUTO_set_runlabels( plint , "Set+Keep" , "Set+Close" ) ;  /* 04 Nov 2003 */
00091 
00092    /*----- Parameters -----*/
00093 
00094    PLUTO_add_option( plint , "Parameters" , "Parameters" , TRUE ) ;
00095 
00096    PLUTO_add_string( plint , "Baseline" , NBASE , baseline_strings , 1 ) ;
00097 
00098    PLUTO_add_number( plint , "Ignore" , 0,20,0,3 , FALSE ) ;
00099 
00100    /*----- Sinusoid -----*/
00101 
00102    for( ii=0 ; ii < NRMAX_SIN ; ii++ ){
00103       PLUTO_add_option( plint , "Sinusoid" , "Sinusoid" , FALSE ) ;
00104       PLUTO_add_number( plint , "Period" , 0,99999,0,20, TRUE ) ;
00105       PLUTO_add_number( plint , "Harmonics" , 1,HARM_MAX,0,1 , FALSE ) ;
00106    }
00107 
00108    /*----- Timeseries -----*/
00109 
00110    for( ii=0 ; ii < NRMAX_TS ; ii++ ){
00111       PLUTO_add_option( plint , "Timeseries" , "Timeseries" , FALSE ) ;
00112       PLUTO_add_timeseries( plint , "File" ) ;
00113    }
00114 
00115    /*--------- done with interface setup ---------*/
00116 
00117    PLUTO_register_1D_funcstr( "L1_Fit" , L1F_fitter ) ;
00118    PLUTO_register_1D_funcstr( "L1_Dtr" , L1F_detrend ) ;
00119 
00120    return plint ;
00121 }

   " Purpose: Control the 'L1_Fit' and 'L1_Dtr 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"