--- imach/src/imach.c 2006/03/14 17:16:22 1.117 +++ imach/src/imach.c 2006/03/22 17:13:53 1.124 @@ -1,6 +1,52 @@ -/* $Id: imach.c,v 1.117 2006/03/14 17:16:22 brouard Exp $ +/* $Id: imach.c,v 1.124 2006/03/22 17:13:53 lievre Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.124 2006/03/22 17:13:53 lievre + Parameters are printed with %lf instead of %f (more numbers after the comma). + The log-likelihood is printed in the log file + + Revision 1.123 2006/03/20 10:52:43 brouard + * imach.c (Module): changed, corresponds to .htm file + name. <head> headers where missing. + + * imach.c (Module): Weights can have a decimal point as for + English (a comma might work with a correct LC_NUMERIC environment, + otherwise the weight is truncated). + Modification of warning when the covariates values are not 0 or + 1. + Version 0.98g + + Revision 1.122 2006/03/20 09:45:41 brouard + (Module): Weights can have a decimal point as for + English (a comma might work with a correct LC_NUMERIC environment, + otherwise the weight is truncated). + Modification of warning when the covariates values are not 0 or + 1. + Version 0.98g + + Revision 1.121 2006/03/16 17:45:01 lievre + * imach.c (Module): Comments concerning covariates added + + * imach.c (Module): refinements in the computation of lli if + status=-2 in order to have more reliable computation if stepm is + not 1 month. Version 0.98f + + Revision 1.120 2006/03/16 15:10:38 lievre + (Module): refinements in the computation of lli if + status=-2 in order to have more reliable computation if stepm is + not 1 month. Version 0.98f + + Revision 1.119 2006/03/15 17:42:26 brouard + (Module): Bug if status = -2, the loglikelihood was + computed as likelihood omitting the logarithm. Version O.98e + + Revision 1.118 2006/03/14 18:20:07 brouard + (Module): varevsij Comments added explaining the second + table of variances if popbased=1 . + (Module): Covariances of eij, ekl added, graphs fixed, new html link. + (Module): Function pstamp added + (Module): Version 0.98d + Revision 1.117 2006/03/14 17:16:22 brouard (Module): varevsij Comments added explaining the second table of variances if popbased=1 . @@ -315,11 +361,11 @@ extern int errno; #define ODIRSEPARATOR '/' #endif -/* $Id: imach.c,v 1.117 2006/03/14 17:16:22 brouard Exp $ */ +/* $Id: imach.c,v 1.124 2006/03/22 17:13:53 lievre Exp $ */ /* $State: Exp $ */ -char version[]="Imach version 0.98d, March 2006, INED-EUROREVES-Institut de longevite "; -char fullversion[]="$Revision: 1.117 $ $Date: 2006/03/14 17:16:22 $"; +char version[]="Imach version 0.98g, March 2006, INED-EUROREVES-Institut de longevite "; +char fullversion[]="$Revision: 1.124 $ $Date: 2006/03/22 17:13:53 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -393,6 +439,7 @@ char strcurr[80], strfor[80]; char *endptr; long lval; +double dval; #define NR_END 1 #define FREE_ARG char* @@ -958,9 +1005,8 @@ void powell(double p[], double **xi, int last_time=curr_time; (void) gettimeofday(&curr_time,&tzp); printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec);fflush(stdout); - /* fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec); - fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tv_sec-start_time.tv_sec); - */ + fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec); fflush(ficlog); +/* fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tv_sec-start_time.tv_sec); */ for (i=1;i<=n;i++) { printf(" %d %.12f",i, p[i]); fprintf(ficlog," %d %.12lf",i, p[i]); @@ -1375,23 +1421,23 @@ double func( double *x) } else if (s2==-2) { for (j=1,survp=0. ; j<=nlstate; j++) - survp += out[s1][j]; - lli= survp; - } - - else if (s2==-4) { - for (j=3,survp=0. ; j<=nlstate; j++) - survp += out[s1][j]; - lli= survp; + survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; + /*survp += out[s1][j]; */ + lli= log(survp); } - else if (s2==-5) { - for (j=1,survp=0. ; j<=2; j++) - survp += out[s1][j]; - lli= survp; - } - + else if (s2==-4) { + for (j=3,survp=0. ; j<=nlstate; j++) + survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; + lli= log(survp); + } + else if (s2==-5) { + for (j=1,survp=0. ; j<=2; j++) + survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; + lli= log(survp); + } + else{ lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ /* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */ @@ -1586,7 +1632,11 @@ double funcone( double *x) */ if( s2 > nlstate && (mle <5) ){ /* Jackson */ lli=log(out[s1][s2] - savm[s1][s2]); - } else if (mle==1){ + } else if (s2==-2) { + for (j=1,survp=0. ; j<=nlstate; j++) + survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; + lli= log(survp); + }else if (mle==1){ lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ } else if(mle==2){ lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ @@ -1602,8 +1652,8 @@ double funcone( double *x) ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; /* printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */ if(globpr){ - fprintf(ficresilk,"%9d %6d %1d %1d %1d %1d %3d %10.6f %6.4f\ - %10.6f %10.6f %10.6f ", \ + fprintf(ficresilk,"%9d %6d %2d %2d %1d %1d %3d %11.6f %8.4f\ + %11.6f %11.6f %11.6f ", \ num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i], 2*weight[i]*lli,out[s1][s2],savm[s1][s2]); for(k=1,llt=0.,l=0.; k<=nlstate; k++){ @@ -2047,6 +2097,7 @@ void freqsummary(char fileres[], int ia for(i=iagemin; i <= iagemax+3; i++){ if(i==iagemax+3){ fprintf(ficlog,"Total"); + fprintf(fichtm,"<br>Total<br>"); }else{ if(first==1){ first=0; @@ -2468,23 +2519,21 @@ void evsij(char fileres[], double ***eij for (age=bage; age<=fage; age ++){ /* If stepm=6 months */ /* Computed by stepm unit matrices, product of hstepm matrices, stored in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */ - + hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij); - + hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ - - /* Computing Variances of health expectancies */ - /* Gradient is computed with plus gp and minus gm. Code is duplicated in order to - decrease memory allocation */ - printf("%d|",(int)age);fflush(stdout); - fprintf(ficlog,"%d|",(int)age);fflush(ficlog); + + printf("%d|",(int)age);fflush(stdout); + fprintf(ficlog,"%d|",(int)age);fflush(ficlog); + /* Computing expectancies */ for(i=1; i<=nlstate;i++) for(j=1; j<=nlstate;j++) for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){ eij[i][j][(int)age] += (p3mat[i][j][h]+p3mat[i][j][h+1])/2.0*hf; -/* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/ + /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/ } @@ -2498,12 +2547,12 @@ void evsij(char fileres[], double ***eij fprintf(ficreseij,"%9.4f", eip ); } fprintf(ficreseij,"\n"); - + } free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); printf("\n"); fprintf(ficlog,"\n"); - + } void cvevsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,double delti[],double **matcov,char strstart[] ) @@ -3435,7 +3484,7 @@ void printinghtml(char fileres[], char t - Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n", subdirf2(fileres,"pl"),subdirf2(fileres,"pl")); fprintf(fichtm,"\ - - (a) Life expectancies by health status at initial age, (b) health expectancies by health status at initial age: ei., eij (estepm=%2d months): \ + - (a) Life expectancies by health status at initial age, (b) health expectancies by health status at initial age: ei., eij . If one or more covariate are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \ <a href=\"%s\">%s</a> <br>\n</li>", estepm,subdirf2(fileres,"e"),subdirf2(fileres,"e")); @@ -4404,7 +4453,7 @@ int main(int argc, char *argv[]) while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0'); printf("val= |%s| pathr=%s\n",val,pathr); strcpy (pathtot, val); - if(pathr[0] == '\0') break; /* Un peu sale */ + if(pathr[0] == '\0') break; /* Dirty */ } } else{ @@ -4423,7 +4472,9 @@ int main(int argc, char *argv[]) /* Split argv[1]=pathtot, parameter file name to get path, optionfile, extension and name */ split(pathtot,path,optionfile,optionfilext,optionfilefiname); printf("\npathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname); - chdir(path); + chdir(path); /* Can be a relative path */ + if(getcwd(pathcd,MAXLINE) > 0) /* So pathcd is the full path */ + printf("Current directory %s!\n",pathcd); strcpy(command,"mkdir "); strcat(command,optionfilefiname); if((outcmd=system(command)) != 0){ @@ -4532,6 +4583,7 @@ int main(int argc, char *argv[]) free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); fclose (ficparo); fclose (ficlog); + goto end; exit(0); } else if(mle==-3) { @@ -4562,7 +4614,9 @@ int main(int argc, char *argv[]) j++; fscanf(ficpar,"%1d%1d",&i1,&j1); if ((i1 != i) && (j1 != j)){ - printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1); + printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n \ +It might be a problem of design; if ncovcol and the model are correct\n \ +run imach with mle=-1 to get a correct template of the parameter file.\n",numlinepar, i,j, i1, j1); exit(1); } fprintf(ficparo,"%1d%1d",i1,j1); @@ -4689,8 +4743,8 @@ int main(int argc, char *argv[]) /*-------- data file ----------*/ if((fic=fopen(datafile,"r"))==NULL) { - printf("Problem with datafile: %s\n", datafile);goto end; - fprintf(ficlog,"Problem with datafile: %s\n", datafile);goto end; + printf("Problem while opening datafile: %s\n", datafile);goto end; + fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);goto end; } n= lastobs; @@ -4786,12 +4840,12 @@ int main(int argc, char *argv[]) cutv(stra, strb,line,' '); errno=0; - lval=strtol(strb,&endptr,10); + dval=strtod(strb,&endptr); if( strb[0]=='\0' || (*endptr != '\0')){ - printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a weight. Exiting.\n",lval, i,line,linei); + printf("Error reading data around '%f' at line number %ld, \"%s\" for individual %d\nShould be a weight. Exiting.\n",dval, i,line,linei); exit(1); } - weight[i]=(double)(lval); + weight[i]=dval; strcpy(line,stra); for (j=ncovcol;j>=1;j--){ @@ -4803,7 +4857,15 @@ int main(int argc, char *argv[]) exit(1); } if(lval <-1 || lval >1){ - printf("Error reading data around '%d' at line number %ld %s for individual %d, '%s'\nShould be a value of the %d covar (meaning 0 for the reference or 1. IMaCh does not build design variables, do it your self). Exiting.\n",lval,linei, i,line,j); + printf("Error reading data around '%d' at line number %ld for individual %d, '%s'\n \ + Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \ + for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \ + For example, for multinomial values like 1, 2 and 3,\n \ + build V1=0 V2=0 for the reference value (1),\n \ + V1=1 V2=0 for (2) \n \ + and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \ + output of IMaCh is often meaningless.\n \ + Exiting.\n",lval,linei, i,line,j); exit(1); } covar[j][i]=(double)(lval); @@ -5099,13 +5161,13 @@ int main(int argc, char *argv[]) printf("Problem with %s \n",optionfilehtmcov), exit(0); } else{ - fprintf(fichtmcov,"<body>\n<title>IMaCh Cov %s\n %s
%s \ + fprintf(fichtmcov,"\nIMaCh Cov %s\n %s
%s \
\n\ Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
\n",\ - fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); + optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); } - fprintf(fichtm,"\nIMaCh %s\n %s
%s \ + fprintf(fichtm,"\nIMaCh %s\n %s
%s \
\n\ Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
\n\ \n\ @@ -5116,7 +5178,7 @@ Title=%s
Datafile=%s Firstpass=%d La - Log file of the run: %s
\n\ - Gnuplot file name: %s
\n\ - Date and time at start: %s\n",\ - fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model,\ + optionfilehtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model,\ optionfilefiname,optionfilext,optionfilefiname,optionfilext,\ fileres,fileres,\ filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart); @@ -5256,7 +5318,7 @@ Interval (in months) between two waves: printf("%d %.0lf %lf %.0lf %.0lf %.0lf %lf\n",k,lsurv[k],p[1]*exp(p[2]*(k-agegomp)),(p[1]*exp(p[2]*(k-agegomp)))*lsurv[k],lpop[k],tpop[k],tpop[k]/lsurv[k]); - replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */ + replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */ printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \ @@ -5299,9 +5361,9 @@ Interval (in months) between two waves: fprintf(ficlog,"%d%d ",i,k); fprintf(ficres,"%1d%1d ",i,k); for(j=1; j <=ncovmodel; j++){ - printf("%f ",p[jk]); - fprintf(ficlog,"%f ",p[jk]); - fprintf(ficres,"%f ",p[jk]); + printf("%lf ",p[jk]); + fprintf(ficlog,"%lf ",p[jk]); + fprintf(ficres,"%lf ",p[jk]); jk++; } printf("\n"); @@ -5493,7 +5555,7 @@ Interval (in months) between two waves: /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint);*/ /*,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ - replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */ + replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */ printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\ @@ -5748,7 +5810,7 @@ Interval (in months) between two waves: } pstamp(ficrest); - fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state ( e.. (std) "); + fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# Age ( e.. (std) "); for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); fprintf(ficrest,"\n"); @@ -5881,14 +5943,20 @@ Interval (in months) between two waves: fprintf(ficlog,"Total time was %d Sec.\n", end_time.tv_sec -start_time.tv_sec); /* printf("Total time was %d uSec.\n", total_usecs);*/ /* if(fileappend(fichtm,optionfilehtm)){ */ - fprintf(fichtm,"
Local time at start %s
Local time at end %s
",strstart, strtend); + fprintf(fichtm,"
Local time at start %s
Local time at end %s
\n",strstart, strtend); fclose(fichtm); + fprintf(fichtmcov,"
Local time at start %s
Local time at end %s
\n",strstart, strtend); fclose(fichtmcov); fclose(ficgp); fclose(ficlog); /*------ End -----------*/ - chdir(path); + + printf("Before Current directory %s!\n",pathcd); + if(chdir(pathcd) != 0) + printf("Can't move to directory %s!\n",path); + if(getcwd(pathcd,MAXLINE) > 0) + printf("Current directory %s!\n",pathcd); /*strcat(plotcmd,CHARSEPARATOR);*/ sprintf(plotcmd,"gnuplot"); #ifndef UNIX