--- imach/src/imach.c 2015/07/16 16:49:02 1.192 +++ imach/src/imach.c 2015/08/18 23:17:52 1.196 @@ -1,6 +1,23 @@ -/* $Id: imach.c,v 1.192 2015/07/16 16:49:02 brouard Exp $ +/* $Id: imach.c,v 1.196 2015/08/18 23:17:52 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.196 2015/08/18 23:17:52 brouard + Summary: 0.98q5 + + Revision 1.195 2015/08/18 16:28:39 brouard + Summary: Adding a hack for testing purpose + + After reading the title, ftol and model lines, if the comment line has + a q, starting with #q, the answer at the end of the run is quit. It + permits to run test files in batch with ctest. The former workaround was + $ echo q | imach foo.imach + + Revision 1.194 2015/08/18 13:32:00 brouard + Summary: Adding error when the covariance matrix doesn't contain the exact number of lines required by the model line. + + Revision 1.193 2015/08/04 07:17:42 brouard + Summary: 0.98q4 + Revision 1.192 2015/07/16 16:49:02 brouard Summary: Fixing some outputs @@ -682,6 +699,7 @@ typedef struct { #define YEARM 12. /**< Number of months per year */ #define AGESUP 130 #define AGEBASE 40 +#define AGEOVERFLOW 1.e20 #define AGEGOMP 10 /**< Minimal age for Gompertz adjustment */ #ifdef _WIN32 #define DIRSEPARATOR '\\' @@ -693,11 +711,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.192 2015/07/16 16:49:02 brouard Exp $ */ +/* $Id: imach.c,v 1.196 2015/08/18 23:17:52 brouard Exp $ */ /* $State: Exp $ */ - -char version[]="Imach version 0.98q3, July 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015"; -char fullversion[]="$Revision: 1.192 $ $Date: 2015/07/16 16:49:02 $"; +#include "version.h" +char version[]=__IMACH_VERSION__; +char copyright[]="August 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015"; +char fullversion[]="$Revision: 1.196 $ $Date: 2015/08/18 23:17:52 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -2891,7 +2910,7 @@ void lubksb(double **a, int n, int *indx void pstamp(FILE *fichier) { - fprintf(fichier,"# %s.%s\n#%s\n#%s\n# %s", optionfilefiname,optionfilext,version,fullversion,strstart); + fprintf(fichier,"# %s.%s\n#IMaCh version %s, %s\n#%s\n# %s", optionfilefiname,optionfilext,version,copyright, fullversion, strstart); } /************ Frequencies ********************/ @@ -4536,12 +4555,12 @@ fprintf(fichtm," \n"); - fprintf(fichtm,"\ \n

  • Result files (second order: variances)

    \n\ - - Parameter file with estimated parameters and covariance matrix: %s
    \n", rfileres,rfileres); + - Parameter file with estimated parameters and covariance matrix: %s
    \ + - 95%% confidence intervals and T statistics are in the log file.
    \n", rfileres,rfileres); - fprintf(fichtm," - Variance of one-step probabilities: %s
    \n", + fprintf(fichtm," - Standard deviation of one-step probabilities: %s
    \n", subdirf2(fileres,"prob"),subdirf2(fileres,"prob")); fprintf(fichtm,"\ - Variance-covariance of one-step probabilities: %s
    \n", @@ -5433,8 +5452,8 @@ int readdata(char datafile[], int firsto if((fic=fopen(datafile,"r"))==NULL) { - printf("Problem while opening datafile: %s\n", datafile);return 1; - fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);return 1; + printf("Problem while opening datafile: %s\n", datafile);fflush(stdout); + fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);fflush(ficlog);return 1; } i=1; @@ -6269,7 +6288,8 @@ int main(int argc, char *argv[]) /* FILE *ficgp;*/ /*Gnuplot File */ struct stat info; double agedeb=0.; - double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20; + + double ageminpar=AGEOVERFLOW,agemin=AGEOVERFLOW, agemaxpar=-AGEOVERFLOW, agemax=-AGEOVERFLOW; double fret; double dum=0.; /* Dummy variable */ @@ -6281,9 +6301,11 @@ int main(int argc, char *argv[]) char pathr[MAXLINE], pathimach[MAXLINE]; char *tok, *val; /* pathtot */ int firstobs=1, lastobs=10; - int c, h , cpt; + int c, h , cpt, c2; int jl=0; int i1, j1, jk, stepsize=0; + int count=0; + int *tab; int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */ int mobilav=0,popforecast=0; @@ -6355,7 +6377,7 @@ int main(int argc, char *argv[]) getcwd(pathcd, size); #endif syscompilerinfo(0); - printf("\n%s\n%s",version,fullversion); + printf("\nIMaCh version %s, %s\n%s",version, copyright, fullversion); if(argc <=1){ printf("\nEnter the parameter file name: "); fgets(pathr,FILENAMELENGTH,stdin); @@ -6474,7 +6496,7 @@ int main(int argc, char *argv[]) ungetc(c,ficpar); fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); - numlinepar++; + numlinepar=numlinepar+3; /* In general */ printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); if(model[strlen(model)-1]=='.') /* Suppressing leading dot in the model */ model[strlen(model)-1]='\0'; @@ -6495,6 +6517,10 @@ int main(int argc, char *argv[]) ungetc(c,ficpar); fgets(line, MAXLINE, ficpar); numlinepar++; + if(line[1]=='q'){ /* This #q will quit imach (the answer is q) */ + z[0]=line[1]; + } + /* printf("****line [1] = %c \n",line[1]); */ fputs(line, stdout); //puts(line); fputs(line,ficparo); @@ -6561,7 +6587,7 @@ int main(int argc, char *argv[]) if(jj==i) continue; j++; fscanf(ficpar,"%1d%1d",&i1,&j1); - if ((i1 != i) && (j1 != j)){ + if ((i1 != i) || (j1 != jj)){ 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); @@ -6569,8 +6595,8 @@ run imach with mle=-1 to get a correct t } fprintf(ficparo,"%1d%1d",i1,j1); if(mle==1) - printf("%1d%1d",i,j); - fprintf(ficlog,"%1d%1d",i,j); + printf("%1d%1d",i,jj); + fprintf(ficlog,"%1d%1d",i,jj); for(k=1; k<=ncovmodel;k++){ fscanf(ficpar," %lf",¶m[i][j][k]); if(mle==1){ @@ -6651,12 +6677,22 @@ run imach with mle=-1 to get a correct t for(i=1; i <=npar; i++) for(j=1; j <=npar; j++) matcov[i][j]=0.; + /* Scans npar lines */ for(i=1; i <=npar; i++){ - fscanf(ficpar,"%s",str); + count=fscanf(ficpar,"%1d%1d%1d",&i1,&j1,&jk); + if(count != 3){ + printf("Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\ +This is probably because your covariance matrix doesn't \n contain exactly %d lines corresponding to your model line '1+age+%s'.\n\ +Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model); + fprintf(ficlog,"Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\ +This is probably because your covariance matrix doesn't \n contain exactly %d lines corresponding to your model line '1+age+%s'.\n\ +Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model); + exit(1); + }else if(mle==1) - printf("%s",str); - fprintf(ficlog,"%s",str); - fprintf(ficparo,"%s",str); + printf("%1d%1d%1d",i1,j1,jk); + fprintf(ficlog,"%1d%1d%1d",i1,j1,jk); + fprintf(ficparo,"%1d%1d%1d",i1,j1,jk); for(j=1; j <=i; j++){ fscanf(ficpar," %le",&matcov[i][j]); if(mle==1){ @@ -6672,6 +6708,7 @@ run imach with mle=-1 to get a correct t fprintf(ficlog,"\n"); fprintf(ficparo,"\n"); } + /* End of read covariance matrix npar lines */ for(i=1; i <=npar; i++) for(j=i+1;j<=npar;j++) matcov[i][j]=matcov[j][i]; @@ -7111,9 +7148,10 @@ Interval (in months) between two waves: } printf("iter=%d MLE=%f Eq=%lf*exp(%lf*(age-%d))\n",iter,-gompertz(p),p[1],p[2],agegomp); - for (i=1;i<=NDIM;i++) + for (i=1;i<=NDIM;i++) { printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); - + fprintf(ficlog,"%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); + } lsurv=vector(1,AGESUP); lpop=vector(1,AGESUP); tpop=vector(1,AGESUP); @@ -7145,8 +7183,15 @@ Interval (in months) between two waves: replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */ - printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); - + if(ageminpar == AGEOVERFLOW ||agemaxpar == AGEOVERFLOW){ + printf("Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ +This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ +Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); + fprintf(ficlog,"Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ +This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ +Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); + }else + printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \ stepm, weightopt,\ model,imx,p,matcov,agemortsup); @@ -7210,6 +7255,24 @@ Interval (in months) between two waves: ftolhess=ftol; /* Usually correct */ hesscov(matcov, p, npar, delti, ftolhess, func); } + printf("Parameters and 95%% confidence intervals\n"); + fprintf(ficlog, "Parameters, T and confidence intervals\n"); + for(i=1,jk=1; i <=nlstate; i++){ + for(k=1; k <=(nlstate+ndeath); k++){ + if (k != i) { + printf("%d%d ",i,k); + fprintf(ficlog,"%d%d ",i,k); + for(j=1; j <=ncovmodel; j++){ + printf("%12.7f T=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-2*sqrt(matcov[jk][jk]),p[jk]+2*sqrt(matcov[jk][jk])); + fprintf(ficlog,"%12.7f T=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-2*sqrt(matcov[jk][jk]),p[jk]+2*sqrt(matcov[jk][jk])); + jk++; + } + printf("\n"); + fprintf(ficlog,"\n"); + } + } + } + fprintf(ficres,"# Scales (for hessian or gradient estimation)\n"); printf("# Scales (for hessian or gradient estimation)\n"); fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n"); @@ -7366,6 +7429,7 @@ Interval (in months) between two waves: dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/365.; fscanf(ficpar,"pop_based=%d\n",&popbased); + fprintf(ficlog,"pop_based=%d\n",popbased); fprintf(ficparo,"pop_based=%d\n",popbased); fprintf(ficres,"pop_based=%d\n",popbased); @@ -7390,7 +7454,15 @@ Interval (in months) between two waves: /* ,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); */ replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */ - printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); + if(ageminpar == AGEOVERFLOW ||agemaxpar == -AGEOVERFLOW){ + printf("Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ +This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ +Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); + fprintf(ficlog,"Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ +This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ +Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); + }else + printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\ model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\ @@ -7793,7 +7865,7 @@ Interval (in months) between two waves: } end: while (z[0] != 'q') { - printf("\nType q for exiting: "); + printf("\nType q for exiting: "); fflush(stdout); scanf("%s",z); } }