--- imach/src/imach.c 2015/09/30 17:45:14 1.203 +++ imach/src/imach.c 2015/10/27 17:36:57 1.207 @@ -1,6 +1,18 @@ -/* $Id: imach.c,v 1.203 2015/09/30 17:45:14 brouard Exp $ +/* $Id: imach.c,v 1.207 2015/10/27 17:36:57 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.207 2015/10/27 17:36:57 brouard + *** empty log message *** + + Revision 1.206 2015/10/24 07:14:11 brouard + *** empty log message *** + + 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 + Revision 1.203 2015/09/30 17:45:14 brouard Summary: looking at better estimation of the hessian @@ -745,12 +757,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.203 2015/09/30 17:45:14 brouard Exp $ */ +/* $Id: imach.c,v 1.207 2015/10/27 17:36:57 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; -char copyright[]="September 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.203 $ $Date: 2015/09/30 17:45:14 $"; +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.207 $ $Date: 2015/10/27 17:36:57 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -940,7 +952,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 ); @@ -1954,7 +1966,24 @@ double **prevalim(double **prlim, int nl { /* Computes the prevalence limit in each live state at age x by left multiplying the unit matrix by transitions matrix until convergence is reached with precision ftolpl */ - + /* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1 = Wx-n Px-n ... Px-2 Px-1 I */ + /* Wx is row vector: population in state 1, population in state 2, population dead */ + /* or prevalence in state 1, prevalence in state 2, 0 */ + /* newm is the matrix after multiplications, its rows are identical at a factor */ + /* Initial matrix pimij */ + /* {0.85204250825084937, 0.13044499163996345, 0.017512500109187184, */ + /* 0.090851990222114765, 0.88271245433047185, 0.026435555447413338, */ + /* 0, 0 , 1} */ + /* + * and after some iteration: */ + /* {0.45504275246439968, 0.42731458730878791, 0.11764266022681241, */ + /* 0.45201005341706885, 0.42865420071559901, 0.11933574586733192, */ + /* 0, 0 , 1} */ + /* And prevalence by suppressing the deaths are close to identical rows in prlim: */ + /* {0.51571254859325999, 0.4842874514067399, */ + /* 0.51326036147820708, 0.48673963852179264} */ + /* If we start from prlim again, prlim tends to a constant matrix */ + int i, ii,j,k; double min, max, maxmin, maxmax,sumnew=0.; /* double **matprod2(); */ /* test */ @@ -2010,13 +2039,13 @@ double **prevalim(double **prlim, int nl prlim[i][j]= newm[i][j]/(1-sumnew); max=FMAX(max,prlim[i][j]); min=FMIN(min,prlim[i][j]); - /* printf(" age= %d prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d max=%f min=%f\n", (int)age, i, j, i, j, prlim[i][j],(int)agefin, max, min); */ + printf(" age= %d prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d max=%f min=%f\n", (int)age, i, j, i, j, prlim[i][j],(int)agefin, max, min); } maxmin=(max-min)/(max+min)*2; maxmax=FMAX(maxmax,maxmin); } /* j loop */ *ncvyear= (int)age- (int)agefin; - /* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */ + printf("maxmax=%lf maxmin=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); if(maxmax < ftolpl){ /* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */ return prlim; @@ -2600,9 +2629,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; @@ -2640,8 +2669,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); @@ -2651,11 +2680,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 first 3 individuals are drawn with lines. The function drawn is -2Log(L) in log scale: %s.png
\ -",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_")); + 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
\ +",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
\ +",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k); + } + } return; } @@ -4264,7 +4305,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; @@ -4323,8 +4364,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 */ @@ -4755,7 +4801,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 *\/ */ @@ -4825,15 +4871,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,""); @@ -4862,14 +4908,27 @@ 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"); -/* good for mle=4 plot by number of matrix products. + 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.png\";",subdirf2(optionfilefiname,"ILK_")); - fprintf(ficgp,"\nplot \"%s\" u 2:(-$11):3 t \"All sample, all transitions\" with dots lc variable",subdirf(fileresilk)); - fprintf(ficgp,"\nreplot \"%s\" u 2:($3 <= 3 ? -$11 : 1/0):3 t \"First 3 individuals\" with line lc variable", subdirf(fileresilk)); + fprintf(ficgp,"\nset out \"%s-dest.png\";",subdirf2(optionfilefiname,"ILK_")); + 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:(-$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 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,",\\\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"); + } + /* unset log; plot "rrtest1_sorted_4/ILK_rrtest1_sorted_4.txt" u 2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with points lc variable */ + /* fprintf(ficgp,"\nset log y;plot \"%s\" u 2:(-$11):3 t \"All sample, all transitions\" with dots lc variable",subdirf(fileresilk)); */ + /* fprintf(ficgp,"\nreplot \"%s\" u 2:($3 <= 3 ? -$11 : 1/0):3 t \"First 3 individuals\" with line lc variable", subdirf(fileresilk)); */ fprintf(ficgp,"\nset out;unset log\n"); /* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */ @@ -6718,14 +6777,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); @@ -7319,7 +7387,7 @@ Please run with mle=-1 to get a correct printf("Problem with file %s",optionfilegnuplot); } else{ - fprintf(ficgp,"\n# %s\n", version); + fprintf(ficgp,"\n# IMaCh-%s\n", version); fprintf(ficgp,"# %s\n", optionfilegnuplot); //fprintf(ficgp,"set missing 'NaNq'\n"); fprintf(ficgp,"set datafile missing 'NaNq'\n"); @@ -7346,13 +7414,15 @@ Please run with mle=-1 to get a correct else{ fprintf(fichtmcov,"\nIMaCh Cov %s\n %s
%s \
\n\ -Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
\n",\ +Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n",\ optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); } - fprintf(fichtm,"\nIMaCh %s\n %s
%s \ + fprintf(fichtm,"\n\n\nIMaCh %s\n IMaCh for Interpolated Markov Chain
\nSponsored by Copyright (C) 2002-2015 INED-EUROREVES-Institut de longévité-Japan Society for the Promotion of Sciences 日本学術振興会 (Grant-in-Aid for Scientific Research 25293121) - Intel Software 2015
\ +
\n\ +IMaCh-%s
%s \
\n\ -Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s
\n\ +Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s
\n\ \n\
\