--- imach/src/imach.c 2015/10/01 16:20:26 1.204 +++ imach/src/imach.c 2015/10/23 15:50:53 1.205 @@ -1,6 +1,9 @@ -/* $Id: imach.c,v 1.204 2015/10/01 16:20:26 brouard Exp $ +/* $Id: imach.c,v 1.205 2015/10/23 15:50:53 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.205 2015/10/23 15:50:53 brouard + Summary: 0.98r3 some clarification for graphs on likelihood contributions + Revision 1.204 2015/10/01 16:20:26 brouard Summary: Some new graphs of contribution to likelihood @@ -748,12 +751,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.204 2015/10/01 16:20:26 brouard Exp $ */ +/* $Id: imach.c,v 1.205 2015/10/23 15:50:53 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; char copyright[]="October 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.204 $ $Date: 2015/10/01 16:20:26 $"; +char fullversion[]="$Revision: 1.205 $ $Date: 2015/10/23 15:50:53 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -943,7 +946,7 @@ static int split( char *path, char *dirc } /* got dirc from getcwd*/ printf(" DIRC = %s \n",dirc); - } else { /* strip direcotry from path */ + } else { /* strip directory from path */ ss++; /* after this, the filename */ l2 = strlen( ss ); /* length of filename */ if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH ); @@ -2603,9 +2606,9 @@ 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,"%9ld %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f\ + fprintf(ficresilk,"%9ld %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f %8.3f\ %11.6f %11.6f %11.6f ", \ - num[i], agexact, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i], + num[i], agexact, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw, 2*weight[i]*lli,out[s1][s2],savm[s1][s2]); for(k=1,llt=0.,l=0.; k<=nlstate; k++){ llt +=ll[k]*gipmx/gsw; @@ -2643,8 +2646,8 @@ void likelione(FILE *ficres,double p[], printf("Problem with resultfile: %s\n", fileresilk); fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk); } - fprintf(ficresilk, "#individual(line's_record) s1 s2 wave# effective_wave# number_of_matrices_product pij weight -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n"); - fprintf(ficresilk, "#num_i age i s1 s2 mi mw dh likeli weight 2wlli out sav "); + fprintf(ficresilk, "#individual(line's_record) count age s1 s2 wave# effective_wave# number_of_matrices_product pij weight weight/gpw -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n"); + fprintf(ficresilk, "#num_i age i s1 s2 mi mw dh likeli weight %weight 2wlli out sav "); /* 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; k<=nlstate; k++) fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k); @@ -2654,19 +2657,23 @@ void likelione(FILE *ficres,double p[], *fretone=(*funcone)(p); if(*globpri !=0){ fclose(ficresilk); - fprintf(fichtm,"\n
File of contributions to the likelihood computed with initial parameters and mle >= 1. You should at least run with mle >= 1 and starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: %s
\n",subdirf(fileresilk),subdirf(fileresilk)); - fprintf(fichtm,"
- The function drawn is -2Log(L) in Log scale: by state of origin %s-ori.png
\ + if (mle ==0) + fprintf(fichtm,"\n
File of contributions to the likelihood computed with initial parameters and mle = %d.",mle); + else if(mle >=1) + fprintf(fichtm,"\n
File of contributions to the likelihood computed with optimized parameters mle = %d.",mle); + fprintf(fichtm," You should at least run with mle >= 1 to get starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: %s
\n",subdirf(fileresilk),subdirf(fileresilk)); + + fprintf(fichtm,"
- The function drawn is -2Log(L) in Log scale: by state of origin %s-ori.png
\ ",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_")); - fprintf(fichtm,"
- and by state of destination %s-dest.png
\ + fprintf(fichtm,"
- and by state of destination %s-dest.png
\ ",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_")); - fflush(fichtm); - - for (k=1; k<= nlstate ; k++) { - fprintf(fichtm,"
- Probability p%dj by origin %d and destination j %s-p%dj.png
\ + fflush(fichtm); + + for (k=1; k<= nlstate ; k++) { + fprintf(fichtm,"
- Probability p%dj by origin %d and destination j %s-p%dj.png
\ ",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k); - - } - } + } + } return; } @@ -4275,7 +4282,7 @@ void cvevsij(double ***eij, double x[], /************ Variance of prevlim ******************/ void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyear, int ij, char strstart[]) { - /* Variance of prevalence limit */ + /* Variance of prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/ /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ double **dnewm,**doldm; @@ -4334,8 +4341,13 @@ void cvevsij(double ***eij, double x[], for(i=1;i<=nlstate;i++) varpl[i][(int)age] =0.; + if((int)age==67 ||(int)age== 66 ){ + matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov); + matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg); + }else{ matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov); matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg); + } for(i=1;i<=nlstate;i++) varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */ @@ -4766,7 +4778,7 @@ divided by h: hPij/h : ", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1); } for(cpt=1; cpt<=nlstate;cpt++) { - fprintf(fichtm,"\n
- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): %s%d%d.svg
\ + fprintf(fichtm,"\n
- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): %s_%d%d.svg
\ ",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1); } /* } /\* end i1 *\/ */ @@ -4836,15 +4848,15 @@ See page 'Matrix of variance-covariance } for(cpt=1; cpt<=nlstate;cpt++) { fprintf(fichtm,"
- Observed (cross-sectional) and period (incidence based) \ -prevalence (with 95%% confidence interval) in state (%d): %s%d_%d.svg
\ -",cpt,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1); +prevalence (with 95%% confidence interval) in state (%d): %s_%d-%d.svg
\ +",cpt,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1); } fprintf(fichtm,"\n
- Total life expectancy by age and \ health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \ true period expectancies (those weighted with period prevalences are also\ drawn in addition to the population based expectancies computed using\ - observed and cahotic prevalences: %s_%d.svg
\ -",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1); + observed and cahotic prevalences: %s_%d.svg
\ +",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1); /* } /\* end i1 *\/ */ }/* End k1 */ fprintf(fichtm,""); @@ -4873,21 +4885,21 @@ void printinggnuplot(char fileresu[], ch fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n"); fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Likelihood (-2Log(L))\";"); /* fprintf(ficgp,"\nset ter svg size 640, 480"); */ /* Too big for svg */ - fprintf(ficgp,"\nset ter png size 640, 480"); + fprintf(ficgp,"\nset ter pngcairo size 640, 480"); /* nice for mle=4 plot by number of matrix products. replot "rrtest1/toto.txt" u 2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with point lc 1 */ /* replot exp(p1+p2*x)/(1+exp(p1+p2*x)+exp(p3+p4*x)+exp(p5+p6*x)) t "p12(x)" */ /* fprintf(ficgp,"\nset out \"%s.svg\";",subdirf2(optionfilefiname,"ILK_")); */ fprintf(ficgp,"\nset out \"%s-dest.png\";",subdirf2(optionfilefiname,"ILK_")); - fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$11):5 t \"All sample, transitions colored by destination\" with dots lc variable; set out;\n",subdirf(fileresilk)); + fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$12):5 t \"All sample, transitions colored by destination\" with dots lc variable; set out;\n",subdirf(fileresilk)); fprintf(ficgp,"\nset out \"%s-ori.png\";",subdirf2(optionfilefiname,"ILK_")); - fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$11):4 t \"All sample, transitions colored by origin\" with dots lc variable; set out;\n\n",subdirf(fileresilk)); + fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$12):4 t \"All sample, transitions colored by origin\" with dots lc variable; set out;\n\n",subdirf(fileresilk)); for (i=1; i<= nlstate ; i ++) { fprintf(ficgp,"\nset out \"%s-p%dj.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i); - fprintf(ficgp,"unset log;\n plot \"%s\"",subdirf(fileresilk)); - fprintf(ficgp," u 2:($4 == %d && $5==%d ? $9 : 1/0):5 t \"p%d%d\" with points lc variable \\\n",i,1,i,1); + fprintf(ficgp,"unset log;\n# plot weighted, mean weight should have point size of 0.5\n plot \"%s\"",subdirf(fileresilk)); + fprintf(ficgp," u 2:($4 == %d && $5==%d ? $9 : 1/0):($11/4.):5 t \"p%d%d\" with points pointtype 7 ps variable lc variable \\\n",i,1,i,1); for (j=2; j<= nlstate+ndeath ; j ++) { - fprintf(ficgp,", \"\" u 2:($4 == %d && $5==%d ? $9 : 1/0):5 t \"p%d%d\" with points lc variable ",i,j,i,j); + fprintf(ficgp,",\\\n \"\" u 2:($4 == %d && $5==%d ? $9 : 1/0):($11/4.):5 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j); } fprintf(ficgp,";\nset out; unset ylabel;\n"); } @@ -6742,14 +6754,23 @@ int main(int argc, char *argv[]) printf("\nIMaCh version %s, %s\n%s",version, copyright, fullversion); if(argc <=1){ printf("\nEnter the parameter file name: "); - fgets(pathr,FILENAMELENGTH,stdin); + if(!fgets(pathr,FILENAMELENGTH,stdin)){ + printf("ERROR Empty parameter file name\n"); + goto end; + } i=strlen(pathr); if(pathr[i-1]=='\n') pathr[i-1]='\0'; i=strlen(pathr); - if(pathr[i-1]==' ') /* This may happen when dragging on oS/X! */ + if(i >= 1 && pathr[i-1]==' ') {/* This may happen when dragging on oS/X! */ pathr[i-1]='\0'; - for (tok = pathr; tok != NULL; ){ + } + i=strlen(pathr); + if( i==0 ){ + printf("ERROR Empty parameter file name\n"); + goto end; + } + for (tok = pathr; tok != NULL; ){ printf("Pathr |%s|\n",pathr); while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0'); printf("val= |%s| pathr=%s\n",val,pathr); @@ -7645,24 +7666,30 @@ Please run with mle=-1 to get a correct free_matrix(ximort,1,NDIM,1,NDIM); #endif } /* Endof if mle==-3 mortality only */ - /* Standard maximisation */ - else{ /* For mle !=- 3 */ - globpr=0;/* debug */ - /* Computes likelihood for initial parameters */ + /* Standard */ + else{ /* For mle !=- 3, could be 0 or 1 or 4 etc. */ + globpr=0;/* Computes sum of likelihood for globpr=1 and funcone */ + /* Computes likelihood for initial parameters, uses funcone to compute gpimx and gsw */ likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); for (k=1; k<=npar;k++) printf(" %d %8.5f",k,p[k]); printf("\n"); - globpr=1; /* again, to print the contributions */ + if(mle>=1){ /* Could be 1 or 2, Real Maximization */ + /* mlikeli uses func not funcone */ + mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func); + } + if(mle==0) {/* No optimization, will print the likelihoods for the datafile */ + globpr=0;/* Computes sum of likelihood for globpr=1 and funcone */ + /* Computes likelihood for initial parameters, uses funcone to compute gpimx and gsw */ + likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ + } + globpr=1; /* again, to print the individual contributions using computed gpimx and gsw */ likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); for (k=1; k<=npar;k++) printf(" %d %8.5f",k,p[k]); printf("\n"); - if(mle>=1){ /* Could be 1 or 2, Real Maximisation */ - mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func); - } /*--------- results files --------------*/ fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model); @@ -8115,6 +8142,7 @@ Please run with mle=-1 to get a correct } fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav); + /* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */ /* printf(" age %4.0f ",age); */ for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){ for(i=1, epj[j]=0.;i <=nlstate;i++) {