skyhawk
2004.08.16, 02:53 PM
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
// 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];
}
}
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
// 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];
}
}