insane math + big numbers

Moderator
Posts: 916
Joined: 2002.10
Post: #1
ok, if you haven't heard my many rants regarding porting this fortran program before, just know that is not currently making me happy.

so here it goes, I have this function that is supposed to calculate the last optional parameters (at least that is how it is going), but I am coming up on some problems, namely regarding infinity and NaN. Perhaps yall can give me suggestions. the biggest problem being with the expm1 and the exp function


Code:
// optional parameters,
// opt = 1, do ll
// opt = 2, do dll_dmu
// opt = 4, do d2ll_dmu2
// and of course all combination of these bits
void ll_derivs(double *dna, double mu, double *seen, double *unseen,int opt,double *ll,double *dll_mu, double *d2ll_dmu2)
{
    int x;
    double max;
    double *mean=(double*)malloc(sizeof(double)*numruns);
    double *pseen=(double*)malloc(sizeof(double)*numruns);
    double *punseen=(double*)malloc(sizeof(double)*numruns);
    double *pdunseen=(double*)malloc(sizeof(double)*numruns);
    printf("mu=%lf\n",mu);
    for(x=0;x<numruns;x++)
    {
        mean[x]=mu*dna[x];
        printf("mean[%d]=%lf ",x,mean[x]);
    }
    printf("\n");
    for(x=0;x<numruns;x++)
    {
        pseen[x]=-expm1(-mean[x]);
        printf("pseen[%d]=%lf ",x,pseen[x]);
    }
    printf("\n");
    // tiny SHOULD be defined by now
    max=tiny;
    for(x=0;x<numruns;x++)
        if(pseen[x]>max);
            max=pseen[x];
    if(max==tiny);
        for(x=0;x<numruns;x++)
            pseen[x]=max;

    for(x=0;x<numruns;x++)
    {
        punseen[x]=exp(-mean[x]);
        printf("punseen[%d]=%lf ",x,punseen[x]);
    }
    printf("\n");
    for(x=0;x<numruns;x++)
        pdunseen[x]=punseen[x]/pseen[x];
    if(opt&1)
    {
        *ll=0;
        for(x=0;x<numruns;x++)
            *ll+=seen[x]*log(pseen[x])-unseen[x]*mean[x];
    }
    if(opt&2)
    {
        *dll_mu=0;
        for(x=0;x<numruns;x++)
            *dll_mu+=expge[x]*(seen[x]*pdunseen[x]-unseen[x]);
    }
    if(opt&4)
    {
        *d2ll_dmu2=0;
        for(x=0;x<numruns;x++)
            *d2ll_dmu2+=seen[x]*pow(dna[x],2)*pdunseen[x]/pseen[x];
    
    }
}
Quote this message in a reply
Nibbie
Posts: 2
Joined: 2006.10
Post: #2
skyhawk Wrote:ok, if you haven't heard my many rants regarding porting this fortran program before, just know that is not currently making me happy.

Can you just use f2c instead? Fink can probably install it for you, or look at the bottom of this page:

http://homepage.mac.com/persquare/octave10_0_4.html

--phillip
Quote this message in a reply
Moderator
Posts: 916
Joined: 2002.10
Post: #3
I doubt it would work with fortran90
Quote this message in a reply
Post Reply