Changeset 786

Show
Ignore:
Timestamp:
01/16/07 10:31:41 (23 months ago)
Author:
peter
Message:

source/service.cpp

  • Bug-fix. Fix problem in e2() with negative return values for very large tau.

source/thirdparty_bevington.cpp
source/conv_eden_ioniz.cpp
source/thirdparty.cpp
source/thirdparty.h

  • Replace Bevington linfit routine with open source version.
  • Bug-fix. Replace call to y0() -> bessel_y0(), y1() -> bessel_y1() in bessel_yn().
Location:
trunk/source
Files:
5 modified

Legend:

Unmodified
Added
Removed

  • trunk/source/conv_eden_ioniz.cpp

    r785 r786  
    263263                                                if( nEval>3 && lgHistUpdate ) 
    264264                                                { 
    265                                                         double y_intercept , stdx , stdy ,reg_coef, 
    266                                                                 slope_save; 
     265                                                        double y_intercept, stdx, stdy, slope_save; 
    267266                                                        slope_save = slope_ls; 
    268267                                                        /* enough data to do least squares */ 
    269                                                         if( linfit(eden_assumed_history,  
    270                                                                 eden_error_history,  
    271                                                                 stdarray,  
    272                                                                 nEval,  
    273                                                                 /* this says same weight for all */ 
    274                                                                 0,  
    275                                                                 &y_intercept,  
    276                                                                 &stdy,  
    277                                                                 &slope_ls,  
    278                                                                 &stdx,  
    279                                                                 &reg_coef) ) 
     268                                                        if( linfit(nEval, 
     269                                                                   eden_assumed_history,  
     270                                                                   eden_error_history,  
     271                                                                   y_intercept,  
     272                                                                   stdy,  
     273                                                                   slope_ls,  
     274                                                                   stdx) ) 
    280275                                                        { 
    281276                                                                int i; 

  • trunk/source/service.cpp

    r785 r786  
    302302{ 
    303303        /* use recurrence relation */ 
    304         return exp_mt - t*ee1(t); 
     304        /* ignore exp_mt, it is *very* unreliable */ 
     305        /* guard against negative results, this can happen for very large t */ 
     306        return max(sexp(t) - t*ee1(t),0.); 
    305307} 
    306308 
    … …  
    336338        { 
    337339                /* abs. accuracy better than 2e-7 */ 
    338 /*              ans = a[0] + a[1]*x + a[2]*x*x + a[3]*powi(x,3) + a[4]*powi(x,4) +  */ 
    339 /*                a[5]*powi(x,5) - log(x); */ 
    340340                ans = ((((a[5]*x + a[4])*x + a[3])*x + a[2])*x + a[1])*x + a[0] - log(x); 
    341341        } 
    … …  
    345345        { 
    346346                /* abs. accuracy better than 2e-8 */ 
    347 /*              top = powi(x,4) + b[0]*powi(x,3) + b[1]*x*x + b[2]*x + b[3]; */ 
    348347                top = (((x + b[0])*x + b[1])*x + b[2])*x + b[3]; 
    349 /*              bot = powi(x,4) + c[0]*powi(x,3) + c[1]*x*x + c[2]*x + c[3]; */ 
    350348                bot = (((x + c[0])*x + c[1])*x + c[2])*x + c[3]; 
    351                 /*  
    352                  *>>>chng 98 dec 17 
     349                /*>>>chng 98 dec 17 
    353350                 * Jason's original implementation did not require the exp since 
    354                  * the routines that used this knew if was not present.   
    355                  * put in for safely. 
    356                   
    357                 ans = top/bot/x;*/ 
     351                 * the routines that used this knew it was not present.   
     352                 * put in for safety. */ 
    358353                ans = top/bot/x*sexp(x); 
    359354        } 

  • trunk/source/thirdparty.cpp

    r785 r786  
    88inline double p1evl(double x, const double coef[], int N); 
    99inline double chbevl(double, const double[], int); 
     10 
     11 
     12/* the routine linfit was derived from the program slopes.f */ 
     13 
     14/********************************************************************/ 
     15/*                                                                  */ 
     16/*                      PROGRAM SLOPES                              */ 
     17/*                                                                  */ 
     18/*    PROGRAM TO COMPUTE THE THEORETICAL REGRESSION SLOPES          */ 
     19/*    AND UNCERTAINTIES (VIA DELTA METHOD), AND UNCERTAINTIES       */ 
     20/*    VIA BOOTSTRAP AND BIVARIATE NORMAL SIMULATION FOR A           */ 
     21/*    (X(I),Y(I)) DATA SET WITH UNKNOWN POPULATION DISTRIBUTION.    */ 
     22/*                                                                  */ 
     23/*       WRITTEN BY ERIC FEIGELSON, PENN STATE.  JUNE 1991          */ 
     24/*                                                                  */ 
     25/*                                                                  */ 
     26/*  A full description of these methods can be found in:            */ 
     27/*    Isobe, T., Feigelson, E. D., Akritas, M. and Babu, G. J.,     */ 
     28/*       Linear Regression in Astronomy I, Astrophys. J. 364, 104   */ 
     29/*       (1990)                                                     */ 
     30/*    Babu, G. J. and Feigelson, E. D., Analytical and Monte Carlo  */ 
     31/*       Comparisons of Six Different Linear Least Squares Fits,    */ 
     32/*       Communications in Statistics, Simulation & Computation,    */ 
     33/*       21, 533 (1992)                                             */ 
     34/*    Feigelson, E. D. and Babu, G. J., Linear Regression in        */ 
     35/*       Astronomy II, Astrophys. J. 397, 55 (1992).                */ 
     36/*                                                                  */ 
     37/********************************************************************/ 
     38 
     39/* this used to be sixlin, but only the first fit has been retained */ 
     40 
     41/********************************************************************/ 
     42/************************* routine linfit ***************************/ 
     43/********************************************************************/ 
     44 
     45bool linfit( 
     46        long n, 
     47        double x[], /* x[n] */ 
     48        double y[], /* y[n] */ 
     49        double &a, 
     50        double &siga, 
     51        double &b, 
     52        double &sigb 
     53) 
     54{ 
     55 
     56/* 
     57 *                       linear regression 
     58 *     written by t. isobe, g. j. babu and e. d. feigelson 
     59 *               center for space research, m.i.t. 
     60 *                             and 
     61 *              the pennsylvania state university 
     62 * 
     63 *                   rev. 1.0,   september 1990 
     64 * 
     65 *       this subroutine provides linear regression coefficients 
     66 *    computed by six different methods described in isobe, 
     67 *    feigelson, akritas, and babu 1990, astrophysical journal 
     68 *    and babu and feigelson 1990, subm. to technometrics. 
     69 *    the methods are ols(y/x), ols(x/y), ols bisector, orthogonal, 
     70 *    reduced major axis, and mean-ols regressions. 
     71 * 
     72 *    [Modified and translated to C/C++ by Peter van Hoof, Royal 
     73 *     Observatory of Belgium; only the first method has been retained] 
     74 *      
     75 * 
     76 *    input 
     77 *         x[i] : independent variable 
     78 *         y[i] : dependent variable 
     79 *            n : number of data points 
     80 * 
     81 *    output 
     82 *         a : intercept coefficients 
     83 *         b : slope coefficients 
     84 *      siga : standard deviations of intercepts 
     85 *      sigb : standard deviations of slopes 
     86 * 
     87 *    error returns 
     88 *         returns true when division by zero will 
     89 *         occur (i.e. when sxy is zero) 
     90 */ 
     91 
     92        DEBUG_ENTRY( "linfit()" ); 
     93 
     94        /* initializations */ 
     95        a = 0.0; 
     96        siga = 0.0; 
     97        b = 0.0; 
     98        sigb = 0.0; 
     99 
     100        /* compute averages and sums */ 
     101        double s1 = 0.0; 
     102        double s2 = 0.0; 
     103        for( long i=0; i < n; i++ ) 
     104        { 
     105                s1 += x[i]; 
     106                s2 += y[i]; 
     107        } 
     108        double rn = (double)n; 
     109        double xavg = s1/rn; 
     110        double yavg = s2/rn; 
     111        double sxx = 0.0; 
     112        double sxy = 0.0; 
     113        for( long i=0; i < n; i++ ) 
     114        { 
     115                x[i] -= xavg; 
     116                y[i] -= yavg; 
     117                sxx += pow2(x[i]); 
     118                sxy += x[i]*y[i]; 
     119        } 
     120 
     121        if( sxy == 0.0 ) 
     122        { 
     123                DEBUG_EXIT( "linfit()" ); 
     124                return true; 
     125        } 
     126 
     127        /* compute the slope coefficient */ 
     128        b = sxy/sxx; 
     129 
     130        /* compute intercept coefficient */ 
     131        a = yavg - b*xavg; 
     132 
     133        /* prepare for computation of variances */ 
     134        double sum1 = 0.0; 
     135        for( long i=0; i < n; i++ ) 
     136                sum1 += pow2(x[i]*(y[i] - b*x[i])); 
     137 
     138        /* compute variance of the slope coefficient */ 
     139        sigb = sum1/pow2(sxx); 
     140 
     141        /* compute variance of the intercept coefficient */ 
     142        for( long i=0; i < n; i++ ) 
     143                siga += pow2((y[i] - b*x[i])*(1.0 - rn*xavg*x[i]/sxx)); 
     144 
     145        /* convert variances to standard deviations */ 
     146        sigb = sqrt(sigb); 
     147        siga = sqrt(siga)/rn; 
     148 
     149        /* return data arrays to their original form */ 
     150        for( long i=0; i < n; i++ ) 
     151        { 
     152                x[i] += xavg; 
     153                y[i] += yavg; 
     154        } 
     155 
     156        DEBUG_EXIT( "linfit()" ); 
     157        return false; 
     158} 
     159 
    10160 
    11161/* the routines factorial and lfactorial came originally from hydro_bauman.cpp 
    … …  
    10641214 
    10651215        /* forward recurrence on n */ 
    1066         anm2 = y0(x); 
    1067         anm1 = y1(x); 
     1216        anm2 = bessel_y0(x); 
     1217        anm1 = bessel_y1(x); 
    10681218        k = 1; 
    10691219        r = 2 * k; 

  • trunk/source/thirdparty.h

    r785 r786  
    99 * that can be represented as double is (NPRE_FACTORIAL-1)! */ 
    1010static const int NPRE_FACTORIAL = 171; 
     11 
     12bool linfit( 
     13        long n, 
     14        double x[], /* x[n] */ 
     15        double y[], /* y[n] */ 
     16        double &a, 
     17        double &siga, 
     18        double &b, 
     19        double &sigb 
     20); 
    1121 
    1222/** factorial: compute n! by lookup in table of predefined factorials */ 
    … …  
    4252double expn(int n, double x); 
    4353 
    44 double expn1(double x); 
    45 double expn2(double x); 
    46  
    4754/*============================================================================*/ 
    4855 
    … …  
    6471        double XIN,  
    6572        double *YOUT 
    66 ); 
    67  
    68 /** 
    69  linfit linear least squares 
    70  \param  x[]   
    71  \param  y[]   
    72  \param  sigmay[]   
    73  \param  npts   
    74  \param  mode   
    75  \param  *a   
    76  \param  *sigmaa   
    77  \param  *b   
    78  \param  *sigmab   
    79  \param  *r 
    80  \return 0 if ok, 1 if disaster 
    81 */  
    82 int linfit( 
    83         double x[],  
    84         double y[],  
    85         double sigmay[],  
    86         long int npts,  
    87         long int mode,  
    88         double *a,  
    89         double *sigmaa,  
    90         double *b,  
    91         double *sigmab,  
    92         double *r 
    9373); 
    9474 

  • trunk/source/thirdparty_bevington.cpp

    r785 r786  
    145145        return; 
    146146} 
    147  
    148 /* linfit do linear least squares fit to arrays, ret 0 if ok, 1 for disaster */ 
    149 int linfit(double x[],  
    150           double y[],  
    151           double sigmay[],  
    152           long int npts,  
    153           long int mode,  
    154           double *yintcp,  
    155           double *sigmaa,  
    156           double *slope,  
    157           double *sigmab,  
    158           double *r) 
    159 { 
    160         /* return error condition */ 
    161         int ierr; 
    162         /* 
    163          * subroutine linfit.f 
    164          * 
    165          * source 
    166          *   Bevington, pages 104-105. 
    167          * 
    168          * purpose 
    169          *   make a least-squares fit to a data set with a straight line 
    170          * 
    171          * usage 
    172          *   call linfit (x, y, sigmay, npts, mode, a, sigmaa, b, sigmab, r) 
    173          * 
    174          * description of parameters 
    175          *   x      - array of data points for independent variable 
    176          *   y      - array of data points for dependent variable 
    177          *   sigmay - array of standard deviations for y data points 
    178          *   npts   - number of pairs of data points 
    179          *   mode   - determines method of weighting least-squares fit 
    180          *            +1 (instrumental) weight(i) = 1./sigmay(i)**2 
    181          *             0 (no weighting) weight(i) = 1. 
    182          *            -1 (statistical)  weight(i) = 1./y(i) 
    183          *   yintcp  - y intercept of fitted straight line 
    184          *   sigmaa - standard deviation of a 
    185          *   slope  - slope of fitted straight line 
    186          *   sigmab - standard deviation of b 
    187          *   r      - linear correlation coefficient 
    188          * 
    189          * subroutines and function subprograms required 
    190          *   none 
    191          * 
    192          */ 
    193         long int i; 
    194         double c; 
    195         double delta,  
    196           sum,  
    197           sumx,  
    198           sumx2,  
    199           sumxy,  
    200           sumy,  
    201           sumy2,  
    202           varnce,  
    203           weight,  
    204           xi,  
    205           yi; 
    206  
    207         DEBUG_ENTRY( "linfit()" ); 
    208  
    209         /* accumulate weighted sums 
    210          * */ 
    211         sum = 0.; 
    212         sumx = 0.; 
    213         sumy = 0.; 
    214         sumx2 = 0.; 
    215         sumxy = 0.; 
    216         sumy2 = 0.; 
    217         for( i=0; i < npts; ++i ) 
    218         { 
    219                 xi = x[i]; 
    220                 yi = y[i]; 
    221                 weight = 1.; 
    222                 if( mode != 0 ) 
    223                 { 
    224                         if( mode > 0 ) 
    225                         { 
    226                                 weight = 1./(sigmay[i]*sigmay[i]); 
    227                         } 
    228                         else if( yi < 0 ) 
    229                         { 
    230                                 weight = 1./(-yi); 
    231                         } 
    232                         else if( yi != 0 ) 
    233                         { 
    234                                 weight = 1./yi; 
    235                         } 
    236                 } 
    237                 sum += weight; 
    238                 sumx += weight*xi; 
    239                 sumy += weight*yi; 
    240                 sumx2 += weight*xi*xi; 
    241                 sumxy += weight*xi*yi; 
    242                 sumy2 += weight*yi*yi; 
    243         } 
    244  
    245         /* calculate coefficients and standard deviations 
    246          * */ 
    247         delta = sum*sumx2 - sumx*sumx; 
    248         if( fabs( delta )>SMALLFLOAT ) 
    249         { 
    250                 /* the y intercept */ 
    251                 *yintcp = (sumx2*sumy - sumx*sumxy)/delta; 
    252                 /* the slope */ 
    253                 *slope = (sumxy*sum - sumx*sumy)/delta; 
    254                 if( mode == 0 ) 
    255                 { 
    256                         c = (double)(npts - 2); 
    257                         varnce = (sumy2 + *yintcp**yintcp*sum + *slope**slope*sumx2 - 2.*(*yintcp*sumy +  
    258                         *slope*sumxy - *yintcp**slope*sumx))/c; 
    259                 } 
    260                 else 
    261                 { 
    262                         varnce = 1.; 
    263                 } 
    264                 /* variance can be slightly negative when points very close together */ 
    265                 *sigmaa = sqrt(MAX2(0.,varnce*sumx2/delta)); 
    266                 *sigmab = sqrt(MAX2(0.,varnce*sum/delta)); 
    267                 delta = MAX2(SMALLFLOAT,delta*(sum*sumy2-sumy*sumy)); 
    268                 *r = (sum*sumxy - sumx*sumy)/sqrt(delta); 
    269                 ierr = 0; 
    270         } 
    271         else 
    272                 ierr = 1; 
    273          
    274         DEBUG_EXIT( "linfit()" ); 
    275         return ierr; 
    276 }