00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #define  NEED_PARSER_INTERNALS
00008 #include "parser.h"
00009 #include "Amalloc.h"
00010 
00011 /***** C routines to interface to the f2c generated parser code *****/
00012 
00013 static int printout = 0 ;
00014 
00015 void PARSER_set_printout( int p ){ printout = p ; }
00016 
00017 /*------------------------------------------------------------------
00018    Input  = expression string
00019    Output = structure containing information about how to
00020             evaluate the expression; should be free-d when unneeded;
00021             if NULL is returned, an error occurred.
00022 --------------------------------------------------------------------*/
00023 
00024 PARSER_code * PARSER_generate_code( char * expression )
00025 {
00026    logical pr ;
00027    integer num_code ;
00028    int nexp ;
00029    PARSER_code * pc ;
00030    char *exp,cc ; int ii,jj ;  /* 22 Jul 2003 */
00031 
00032    if( expression == NULL ) return NULL ;
00033    nexp = strlen( expression ) ;
00034    if( nexp == 0 ) return NULL ;
00035 
00036    /* 22 Jul 2003: copy into local string, tossing bad stuff */
00037 
00038    exp = AFMALL(char, nexp+4) ;
00039    for( ii=jj=0 ; ii < nexp ; ii++ ){
00040      cc = expression[ii] ;
00041      if( !isspace(cc) && !iscntrl(cc) ) exp[jj++] = cc ;
00042    }
00043    exp[jj] = '\0' ;
00044    nexp = strlen(exp) ; if( nexp == 0 ) return NULL ;
00045 
00046    pc = (PARSER_code *) malloc( sizeof(PARSER_code) ) ;
00047 
00048    pr = (printout) ? TRUE_ : FALSE_ ;
00049 
00050    parser_( exp, &pr, &num_code, pc->c_code, (ftnlen) nexp, (ftnlen) 8 ) ;
00051 
00052    free(exp) ;  /* 22 Jul 2003 */
00053 
00054    if( num_code <= 0 ){ free(pc) ; return NULL ; }
00055 
00056    pc->num_code = (int) num_code ;
00057    return pc ;
00058 }
00059 
00060 /*---------------------------------------------------------------
00061    pc   = code to evaluate expression
00062    atoz = double [26] containing values for variables A,B,...,Z
00063 -----------------------------------------------------------------*/
00064 
00065 double PARSER_evaluate_one( PARSER_code * pc , double atoz[] )
00066 {
00067    integer num_code ;
00068    double  value ;
00069 
00070    if( pc == NULL || pc->num_code <= 0 ) return 0.0 ;
00071 
00072    num_code = (integer) pc->num_code ;
00073 
00074    value = (double) pareval_( &num_code, pc->c_code,
00075                               (doublereal *) atoz , (ftnlen) 8 ) ;
00076    return value ;
00077 }
00078 
00079 /*----------------------------------------------------------------------
00080    Return 1 if the given code uses the symbol given by the
00081    first character of sym, otherwise return 0 - 15 Sep 1999 - RWCox.
00082 ------------------------------------------------------------------------*/
00083 
00084 extern integer hassym_(char *sym, integer *num_code__, char *c_code__, ftnlen
00085         sym_len, ftnlen c_code_len) ;
00086 
00087 int PARSER_has_symbol( char * sym , PARSER_code * pc )
00088 {
00089    int hh ;
00090    char sss[8] ;
00091    integer num_code ;
00092 
00093    if( !isalpha(sym[0]) ) return 0 ;          /* not alphabetic */
00094 
00095    sss[0] = toupper(sym[0]) ; sss[1] = '\0' ; /* uppercase it */
00096 
00097    num_code = (integer) pc->num_code ;
00098 
00099    hh = (int) hassym_( sss , &num_code , pc->c_code ,
00100                        (ftnlen) 8 , (ftnlen) 8       ) ;
00101 
00102    return hh ;
00103 }
00104 
00105 void PARSER_mark_symbols( PARSER_code * pc , int * sl )
00106 {
00107    int ii ;
00108    static char abet[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZ" ;
00109 
00110    if( pc == NULL || sl == NULL ) return ;
00111 
00112    for( ii=0 ; ii < 26 ; ii++ )
00113       sl[ii] = PARSER_has_symbol( abet+ii , pc ) ;
00114 
00115    return ;
00116 }
00117 
00118 /*----------------------------------------------------------------------
00119    pc    = code to evaluate expression
00120    atoz  = double * [26] containing values for variables A..Z
00121    nv    = length of vectors
00122    vout  = double [nv]; will get the output of the expression
00123              evaluated using the values in atoz
00124 
00125    The only reason for calling this routine instead of
00126    PARSER_evaluate_one nv different times is efficiency.
00127 ------------------------------------------------------------------------*/
00128 
00129 void PARSER_evaluate_vector( PARSER_code * pc , double* atoz[] ,
00130                              int nv , double vout[] )
00131 {
00132    integer num_code , nvar , ivar , lvec , ldvec ;
00133 
00134    if( pc == NULL || pc->num_code <= 0 ) return ;
00135 
00136    num_code = (integer) pc->num_code ;
00137    lvec     = (integer) nv ;
00138 
00139    parevec_( &num_code , pc->c_code ,
00140              (doublereal *) atoz[0]  , (doublereal *) atoz[1] ,
00141              (doublereal *) atoz[2]  , (doublereal *) atoz[3] ,
00142              (doublereal *) atoz[4]  , (doublereal *) atoz[5] ,
00143              (doublereal *) atoz[6]  , (doublereal *) atoz[7] ,
00144              (doublereal *) atoz[8]  , (doublereal *) atoz[9] ,
00145              (doublereal *) atoz[10] , (doublereal *) atoz[11] ,
00146              (doublereal *) atoz[12] , (doublereal *) atoz[13] ,
00147              (doublereal *) atoz[14] , (doublereal *) atoz[15] ,
00148              (doublereal *) atoz[16] , (doublereal *) atoz[17] ,
00149              (doublereal *) atoz[18] , (doublereal *) atoz[19] ,
00150              (doublereal *) atoz[20] , (doublereal *) atoz[21] ,
00151              (doublereal *) atoz[22] , (doublereal *) atoz[23] ,
00152              (doublereal *) atoz[24] , (doublereal *) atoz[25] ,
00153          &lvec , (doublereal *) vout , (ftnlen) 8 ) ;
00154 
00155    return ;
00156 }
00157 
00158 /*----------------------------------------------------------------------
00159    Evaluate an expression at a set of evenly spaced points.
00160      expr = expression - first symbol found will be the variable
00161      nt   = number of points
00162      tz   = value of first point (the variable)
00163      dt   = spacing between points
00164      vec  = pointer to pre-allocated output location of length nt
00165    Return value is 1 for good results, 0 for errors.
00166 
00167    17 Nov 1999 - RW Cox - adapted from 1deval.c [hence the name]
00168 ------------------------------------------------------------------------*/
00169 
00170 int PARSER_1deval( char * expr, int nt, float tz, float dt, float * vec )
00171 {
00172    PARSER_code * pcode = NULL ;
00173    char sym[4] ;
00174    double atoz[26] ;
00175    int ii , kvar ;
00176 
00177    if( expr == NULL || nt <= 0 || vec == NULL ) return 0 ;  /* bad */
00178 
00179    pcode = PARSER_generate_code( expr ) ;        /* compile */
00180    if( pcode == NULL ) return 0 ;                /* bad news */
00181 
00182    kvar = -1 ;                                   /* find symbol */
00183    for( ii=0 ; ii < 26 ; ii++ ){
00184       sym[0] = 'A' + ii ; sym[1] = '\0' ;
00185       if( PARSER_has_symbol(sym,pcode) ){ kvar = ii ; break ; }
00186    }
00187 
00188    for( ii=0 ; ii < 26 ; ii++ ) atoz[ii] = 0.0 ; /* initialize */
00189 
00190    if( kvar >= 0 ){                              /* the normal case */
00191       for( ii=0 ; ii < nt ; ii++ ){
00192          atoz[kvar] = tz + ii*dt ;
00193          vec[ii]    = PARSER_evaluate_one( pcode , atoz ) ;
00194       }
00195    } else {                                      /* no variable found! */
00196       vec[0] = PARSER_evaluate_one( pcode , atoz ) ;
00197       for( ii=1 ; ii < nt ; ii++ ) vec[ii] = vec[0] ;
00198    }
00199 
00200    free(pcode) ; return 1 ;
00201 }
00202 
00203 /*------------------------------------------------------------------------*/
00204 /*! Sort of like strtod(), but with arithmetic -- 03 Sep 2002 - RWCox.
00205 --------------------------------------------------------------------------*/
00206 
00207 double PARSER_strtod( char *expr )
00208 {
00209    PARSER_code * pcode = NULL ;
00210    double atoz[26] , val ;
00211    int ii ;
00212 
00213    if( expr == NULL ) return 0 ;                 /* bad */
00214 
00215    pcode = PARSER_generate_code( expr ) ;        /* compile */
00216    if( pcode == NULL ) return 0 ;                /* bad news */
00217 
00218    for( ii=0 ; ii < 26 ; ii++ ) atoz[ii] = 0.0 ; /* initialize */
00219 
00220    val = PARSER_evaluate_one( pcode , atoz ) ;
00221 
00222    free(pcode) ; return val ;
00223 }
00224 
00225 /********************************************************************/
00226 /*** use the C math library to provide Bessel and error functions ***/
00227 
00228 doublereal dbesj0_( doublereal * x )
00229 { return (doublereal) j0( (double) *x ) ; }
00230 
00231 doublereal dbesj1_( doublereal * x )
00232 { return (doublereal) j1( (double) *x ) ; }
00233 
00234 doublereal dbesy0_( doublereal * x )
00235 { return (doublereal) (*x>0) ? y0( (double) *x ) : 0.0 ; }
00236 
00237 doublereal dbesy1_( doublereal * x )
00238 { return (doublereal) (*x>0) ? y1( (double) *x ) : 0.0 ; }
00239 
00240 doublereal derf_ ( doublereal * x )
00241 { return (doublereal) erf( (double) *x ) ; }
00242 
00243 doublereal derfc_( doublereal * x )
00244 { return (doublereal) erfc( (double) *x ) ; }
00245 
00246 #include <time.h>
00247 
00248 doublereal unif_( doublereal * x )  /* 04 Feb 2000 */
00249 {
00250    static first=1 ;
00251    doublereal val ;
00252    if( first ){ srand48((long)time(NULL)); first=0; }
00253    val = (doublereal) drand48() ;
00254    return val ;
00255 }
00256 
00257 /********************************************************************/
00258 /*** Legendre Polynomials (0 <= mm <= 20 ; -1 <= xx < 1)          ***/
00259 
00260 doublereal legendre_( doublereal *mm , doublereal *xx )
00261 {
00262    double x = *xx ;
00263    int    m = (int)(*mm) ;
00264 
00265    if( m < 0 ) return 1.0 ;    /* bad input */
00266 
00267    switch( m ){
00268     case 0: return 1.0 ;
00269     case 1: return x ;
00270     case 2: return (3.0*x*x-1.0)/2.0 ;
00271     case 3: return (5.0*x*x-3.0)*x/2.0 ;
00272     case 4: return ((35.0*x*x-30.0)*x*x+3.0)/8.0 ;
00273     case 5: return ((63.0*x*x-70.0)*x*x+15.0)*x/8.0 ;
00274     case 6: return (((231.0*x*x-315.0)*x*x+105.0)*x*x-5.0)/16.0 ;
00275     case 7: return (((429.0*x*x-693.0)*x*x+315.0)*x*x-35.0)*x/16.0 ;
00276     case 8: return ((((6435.0*x*x-12012.0)*x*x+6930.0)*x*x-1260.0)*x*x+35.0)/128.0;
00277 
00278            /** 07 Feb 2005: this part generated by Maple, then hand massaged **/
00279 
00280     case 9:
00281       return (0.24609375e1 + (-0.3609375e2 + (0.140765625e3 + (-0.20109375e3
00282               + 0.949609375e2 * x * x) * x * x) * x * x) * x * x) * x;
00283 
00284     case 10:
00285       return -0.24609375e0 + (0.1353515625e2 + (-0.1173046875e3 +
00286               (0.3519140625e3 + (-0.42732421875e3 + 0.18042578125e3 * x * x)
00287              * x * x) * x * x) * x * x) * x * x;
00288 
00289     case 11:
00290       return (-0.270703125e1 + (0.5865234375e2 + (-0.3519140625e3 +
00291              (0.8546484375e3 + (-0.90212890625e3 + 0.34444921875e3 * x * x)
00292              * x * x) * x * x) * x * x) * x * x) * x;
00293 
00294     case 12:
00295       return 0.2255859375e0 + (-0.17595703125e2 + (0.2199462890625e3 +
00296              (-0.99708984375e3 + (0.20297900390625e4 + (-0.1894470703125e4
00297              + 0.6601943359375e3 * x * x) * x * x) * x * x) * x * x) * x * x)
00298              * x * x;
00299 
00300     case 13:
00301       return (0.29326171875e1 + (-0.87978515625e2 + (0.7478173828125e3 +
00302              (-0.270638671875e4 + (0.47361767578125e4 + (-0.3961166015625e4
00303              + 0.12696044921875e4 * x * x) * x * x) * x * x) * x * x) * x * x)
00304             * x * x) * x;
00305 
00306     case 14:
00307       return -0.20947265625e0 + (0.2199462890625e2 + (-0.37390869140625e3 +
00308              (0.236808837890625e4 + (-0.710426513671875e4 +
00309              (0.1089320654296875e5 + (-0.825242919921875e4 +
00310             0.244852294921875e4 * x * x) * x * x) * x * x) * x * x) * x * x)
00311            * x * x) * x * x;
00312 
00313     case 15:
00314       return (-0.314208984375e1 + (0.12463623046875e3 + (-0.142085302734375e4
00315             + (0.710426513671875e4 + (-0.1815534423828125e5 +
00316               (0.2475728759765625e5 + (-0.1713966064453125e5 +
00317                0.473381103515625e4 * x * x) * x * x) * x * x) * x * x)
00318              * x * x) * x * x) * x * x) * x;
00319 
00320     case 16:
00321       return 0.196380615234375e0 + (-0.26707763671875e2 + (0.5920220947265625e3
00322             + (-0.4972985595703125e4 + (0.2042476226806641e5 +
00323               (-0.4538836059570312e5 + (0.5570389709472656e5 +
00324                (-0.3550358276367188e5 + 0.9171758880615234e4 * x * x) * x * x)
00325             * x * x) * x * x) * x * x) * x * x) * x * x) * x * x;
00326 
00327     case 17:
00328       return (0.3338470458984375e1 + (-0.169149169921875e3 +
00329              (0.2486492797851562e4 + (-0.1633980981445312e5 +
00330              (0.5673545074462891e5 + (-0.1114077941894531e6 +
00331              (0.1242625396728516e6 + (-0.7337407104492188e5 +
00332               0.1780400253295898e5 * x * x) * x * x) * x * x) * x * x)
00333            * x * x) * x * x) * x * x) * x * x) * x;
00334 
00335     case 18:
00336       return -0.1854705810546875e0 + (0.3171546936035156e2 +
00337             (-0.8880331420898438e3 + (0.9531555725097656e4 +
00338             (-0.5106190567016602e5 + (0.153185717010498e6 +
00339             (-0.2692355026245117e6 + (0.275152766418457e6 +
00340             (-0.1513340215301514e6 + 0.3461889381408691e5 * x * x) * x * x)
00341            * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x * x;
00342 
00343     case 19:
00344       return (-0.3523941040039062e1 + (0.2220082855224609e3 +
00345              (-0.4084952453613281e4 + (0.3404127044677734e5 +
00346              (-0.153185717010498e6 + (0.4038532539367676e6 +
00347              (-0.6420231216430664e6 + (0.6053360861206055e6 +
00348              (-0.3115700443267822e6 + 0.6741574058532715e5 * x * x) * x * x)
00349           * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x;
00350 
00351     case 20:
00352       return 0.1761970520019531e0 + (-0.3700138092041016e2 +
00353             (0.127654764175415e4 + (-0.1702063522338867e5 +
00354             (0.1148892877578735e6 + (-0.4442385793304443e6 +
00355             (0.1043287572669983e7 + (-0.1513340215301514e7 +
00356             (0.1324172688388824e7 + (-0.6404495355606079e6 +
00357              0.1314606941413879e6 * x * x) * x * x) * x * x) * x * x) * x * x)
00358             * x * x) * x * x) * x * x) * x * x) * x * x;
00359    }
00360 
00361         /** if here, m > 20 ==> use recurrence relation **/
00362 
00363    { double k , pk, pkm1, pkm2 ;
00364 
00365      k = 19.0; pkm2 = legendre_( &k , xx ) ;
00366      k = 20.0; pkm1 = legendre_( &k , xx ) ;
00367 
00368      while( k < m ){
00369        k += 1.0 ;
00370        pk = ((2.0*k-1.0)*x*pkm1 - (k-1.0)*pkm2)/k ;
00371        pkm2 = pkm1 ; pkm1 = pk ;
00372      }
00373      return pk ;
00374    }
00375 }
00376 
00377 /**** statistic conversion routines ****/
00378 
00379 /*--- macros to create functions ---*/
00380 
00381 #undef FUNC3
00382 #undef FUNC2
00383 #undef FUNC1
00384 #undef FUNC0
00385 
00386 #define FUNC4(fname,fcode,sfunc)                                      \
00387    doublereal fname( doublereal * x, doublereal * a, doublereal * b,  \
00388                                                      doublereal * c ) \
00389    {  float aux[3] , xx , val ;                                       \
00390       xx     = (float) (*x) ;                                         \
00391       aux[0] = (float) (*a) ;                                         \
00392       aux[1] = (float) (*b) ;                                         \
00393       aux[2] = (float) (*c) ;                                         \
00394       val = sfunc( xx , fcode , aux ) ;                               \
00395       return (doublereal) val ;                                       \
00396    }
00397 
00398 #define FUNC3(fname,fcode,sfunc)                                      \
00399    doublereal fname( doublereal * x, doublereal * a, doublereal * b ) \
00400    {  float aux[2] , xx , val ;                                       \
00401       xx     = (float) (*x) ;                                         \
00402       aux[0] = (float) (*a) ;                                         \
00403       aux[1] = (float) (*b) ;                                         \
00404       val = sfunc( xx , fcode , aux ) ;                               \
00405       return (doublereal) val ;                                       \
00406    }
00407 
00408 #define FUNC2(fname,fcode,sfunc)                       \
00409    doublereal fname( doublereal * x , doublereal * a ) \
00410    {  float aux[1] , xx , val ;                        \
00411       xx     = (float) (*x) ;                          \
00412       aux[0] = (float) (*a) ;                          \
00413       val = sfunc( xx , fcode , aux ) ;                \
00414       return (doublereal) val ;                        \
00415    }
00416 
00417 #define FUNC1(fname,fcode,sfunc)         \
00418    doublereal fname( doublereal * x )    \
00419    {  float  xx , val ;                  \
00420       xx     = (float) (*x) ;            \
00421       val = sfunc( xx , fcode , NULL ) ; \
00422       return (doublereal) val ;          \
00423    }
00424 
00425 #define FUNC_COR_TYPE   2
00426 #define FUNC_TT_TYPE    3
00427 #define FUNC_FT_TYPE    4
00428 #define FUNC_ZT_TYPE    5
00429 #define FUNC_CT_TYPE    6
00430 #define FUNC_BT_TYPE    7
00431 #define FUNC_BN_TYPE    8
00432 #define FUNC_GT_TYPE    9
00433 #define FUNC_PT_TYPE   10
00434 
00435 extern float THD_stat_to_pval(float,int,float *) ;
00436 extern float THD_pval_to_stat(float,int,float *) ;
00437 extern float THD_stat_to_zscore(float,int,float *) ;
00438 
00439 /*--- now create the functions ---*/
00440 
00441 FUNC4(ficotp_,FUNC_COR_TYPE,THD_stat_to_pval)
00442 FUNC4(ficopt_,FUNC_COR_TYPE,THD_pval_to_stat)
00443 FUNC4(ficotz_,FUNC_COR_TYPE,THD_stat_to_zscore)
00444 
00445 FUNC2(fitttp_,FUNC_TT_TYPE,THD_stat_to_pval)
00446 FUNC2(fittpt_,FUNC_TT_TYPE,THD_pval_to_stat)
00447 FUNC2(fitttz_,FUNC_TT_TYPE,THD_stat_to_zscore)
00448 
00449 FUNC3(fifttp_,FUNC_FT_TYPE,THD_stat_to_pval)
00450 FUNC3(fiftpt_,FUNC_FT_TYPE,THD_pval_to_stat)
00451 FUNC3(fifttz_,FUNC_FT_TYPE,THD_stat_to_zscore)
00452 
00453 FUNC1(fizttp_,FUNC_ZT_TYPE,THD_stat_to_pval)
00454 FUNC1(fiztpt_,FUNC_ZT_TYPE,THD_pval_to_stat)
00455 FUNC1(fizttz_,FUNC_ZT_TYPE,THD_stat_to_zscore)
00456 
00457 FUNC2(ficttp_,FUNC_CT_TYPE,THD_stat_to_pval)
00458 FUNC2(fictpt_,FUNC_CT_TYPE,THD_pval_to_stat)
00459 FUNC2(ficttz_,FUNC_CT_TYPE,THD_stat_to_zscore)
00460 
00461 FUNC3(fibttp_,FUNC_BT_TYPE,THD_stat_to_pval)
00462 FUNC3(fibtpt_,FUNC_BT_TYPE,THD_pval_to_stat)
00463 FUNC3(fibttz_,FUNC_BT_TYPE,THD_stat_to_zscore)
00464 
00465 FUNC3(fibntp_,FUNC_BN_TYPE,THD_stat_to_pval)
00466 FUNC3(fibnpt_,FUNC_BN_TYPE,THD_pval_to_stat)
00467 FUNC3(fibntz_,FUNC_BN_TYPE,THD_stat_to_zscore)
00468 
00469 FUNC3(figttp_,FUNC_GT_TYPE,THD_stat_to_pval)
00470 FUNC3(figtpt_,FUNC_GT_TYPE,THD_pval_to_stat)
00471 FUNC3(figttz_,FUNC_GT_TYPE,THD_stat_to_zscore)
00472 
00473 FUNC2(fipttp_,FUNC_PT_TYPE,THD_stat_to_pval)
00474 FUNC2(fiptpt_,FUNC_PT_TYPE,THD_pval_to_stat)
00475 FUNC2(fipttz_,FUNC_PT_TYPE,THD_stat_to_zscore)