--- imach/src/imach.c 2017/07/17 08:53:49 1.277 +++ imach/src/imach.c 2018/02/21 07:58:13 1.280 @@ -1,6 +1,17 @@ -/* $Id: imach.c,v 1.277 2017/07/17 08:53:49 brouard Exp $ +/* $Id: imach.c,v 1.280 2018/02/21 07:58:13 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.280 2018/02/21 07:58:13 brouard + Summary: 0.99r15 + + New Makefile with recent VirtualBox 5.26. Bug in sqrt negatve in imach.c + + Revision 1.279 2017/07/20 13:35:01 brouard + Summary: temporary working + + Revision 1.278 2017/07/19 14:09:02 brouard + Summary: Bug for mobil_average=0 and prevforecast fixed(?) + Revision 1.277 2017/07/17 08:53:49 brouard Summary: BOM files can be read now @@ -1034,12 +1045,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.277 2017/07/17 08:53:49 brouard Exp $ */ +/* $Id: imach.c,v 1.280 2018/02/21 07:58:13 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; char copyright[]="February 2016,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018"; -char fullversion[]="$Revision: 1.277 $ $Date: 2017/07/17 08:53:49 $"; +char fullversion[]="$Revision: 1.280 $ $Date: 2018/02/21 07:58:13 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -2513,15 +2524,18 @@ void powell(double p[], double **xi, int double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int *ncvyear, int ij, int nres) { - /* Computes the prevalence limit in each live state at age x and for covariate combination ij - (and selected quantitative values in nres) - 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 */ + /**< Computes the prevalence limit in each live state at age x and for covariate combination ij + * (and selected quantitative values in nres) + * 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. + * Inputs are the parameter, age, a tolerance for the prevalence limit ftolpl. + * Output is prlim. + * Initial matrix pimij + */ /* {0.85204250825084937, 0.13044499163996345, 0.017512500109187184, */ /* 0.090851990222114765, 0.88271245433047185, 0.026435555447413338, */ /* 0, 0 , 1} */ @@ -5769,10 +5783,11 @@ void concatwav(int wav[], int **dh, int /************ Variance ******************/ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[], int nres) { - /* Variance of health expectancies */ - /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/ - /* double **newm;*/ - /* int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav)*/ + /** Variance of health expectancies + * double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl); + * double **newm; + * int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav) + */ /* int movingaverage(); */ double **dnewm,**doldm; @@ -5780,11 +5795,11 @@ void concatwav(int wav[], int **dh, int int i, j, nhstepm, hstepm, h, nstepm ; int k; double *xp; - double **gp, **gm; /* for var eij */ - double ***gradg, ***trgradg; /*for var eij */ - double **gradgp, **trgradgp; /* for var p point j */ - double *gpp, *gmp; /* for var p point j */ - double **varppt; /* for var p point j nlstate to nlstate+ndeath */ + double **gp, **gm; /**< for var eij */ + double ***gradg, ***trgradg; /**< for var eij */ + double **gradgp, **trgradgp; /**< for var p point j */ + double *gpp, *gmp; /**< for var p point j */ + double **varppt; /**< for var p point j nlstate to nlstate+ndeath */ double ***p3mat; double age,agelim, hf; /* double ***mobaverage; */ @@ -5845,7 +5860,7 @@ void concatwav(int wav[], int **dh, int /* fprintf(fichtm, "#Local time at start: %s", strstart);*/ fprintf(fichtm,"\n
  • Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)

  • \n"); fprintf(fichtm,"\n
    %s
    \n",digitp); - /* } */ + varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); pstamp(ficresvij); fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are "); @@ -5900,9 +5915,12 @@ void concatwav(int wav[], int **dh, int for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/ xp[i] = x[i] + (i==theta ?delti[theta]:0); } - + /**< Computes the prevalence limit with parameter theta shifted of delta up to ftolpl precision and + * returns into prlim . + */ prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij, nres); - + + /* If popbased = 1 we use crossection prevalences. Previous step is useless but prlim is created */ if (popbased==1) { if(mobilav ==0){ for(i=1; i<=nlstate;i++) @@ -5912,23 +5930,28 @@ void concatwav(int wav[], int **dh, int prlim[i][i]=mobaverage[(int)age][i][ij]; } } - - hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres); /* Returns p3mat[i][j][h] for h=1 to nhstepm */ + /**< Computes the shifted transition matrix \f$ {}{h}_p^{ij}_x\f$ at horizon h. + */ + hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres); /* Returns p3mat[i][j][h] for h=0 to nhstepm */ + /**< And for each alive state j, sums over i \f$ w^i_x {}{h}_p^{ij}_x\f$, which are the probability + * at horizon h in state j including mortality. + */ for(j=1; j<= nlstate; j++){ for(h=0; h<=nhstepm; h++){ for(i=1, gp[h][j]=0.;i<=nlstate;i++) gp[h][j] += prlim[i][i]*p3mat[i][j][h]; } } - /* Next for computing probability of death (h=1 means + /* Next for computing shifted+ probability of death (h=1 means computed over hstepm matrices product = hstepm*stepm months) - as a weighted average of prlim. + as a weighted average of prlim(i) * p(i,j) p.3=w1*p13 + w2*p23 . */ for(j=nlstate+1;j<=nlstate+ndeath;j++){ for(i=1,gpp[j]=0.; i<= nlstate; i++) gpp[j] += prlim[i][i]*p3mat[i][j][1]; - } - /* end probability of death */ + } + + /* Again with minus shift */ for(i=1; i<=npar; i++) /* Computes gradient x - delta */ xp[i] = x[i] - (i==theta ?delti[theta]:0); @@ -5961,19 +5984,23 @@ void concatwav(int wav[], int **dh, int for(i=1,gmp[j]=0.; i<= nlstate; i++) gmp[j] += prlim[i][i]*p3mat[i][j][1]; } - /* end probability of death */ - + /* end shifting computations */ + + /**< Computing gradient matrix at horizon h + */ for(j=1; j<= nlstate; j++) /* vareij */ for(h=0; h<=nhstepm; h++){ gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta]; } - - for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu */ + /**< Gradient of overall mortality p.3 (or p.j) + */ + for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu mortality from j */ gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta]; } } /* End theta */ - + + /* We got the gradient matrix for each theta and state j */ trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); /* veij */ for(h=0; h<=nhstepm; h++) /* veij */ @@ -5984,13 +6011,19 @@ void concatwav(int wav[], int **dh, int for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */ for(theta=1; theta <=npar; theta++) trgradgp[j][theta]=gradgp[theta][j]; - + /**< as well as its transposed matrix + */ hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ for(i=1;i<=nlstate;i++) for(j=1;j<=nlstate;j++) vareij[i][j][(int)age] =0.; - + + /* Computing trgradg by matcov by gradg at age and summing over h + * and k (nhstepm) formula 15 of article + * Lievre-Brouard-Heathcote + */ + for(h=0;h<=nhstepm;h++){ for(k=0;k<=nhstepm;k++){ matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov); @@ -6001,7 +6034,11 @@ void concatwav(int wav[], int **dh, int } } - /* pptj */ + /* pptj is p.3 or p.j = trgradgp by cov by gradgp, variance of + * p.j overall mortality formula 49 but computed directly because + * we compute the grad (wix pijx) instead of grad (pijx),even if + * wix is independent of theta. + */ matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov); matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp); for(j=nlstate+1;j<=nlstate+ndeath;j++) @@ -6610,7 +6647,12 @@ To be simple, these graphs help to under } /* Eigen vectors */ - v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12)); + if(1+(v1-lc1)*(v1-lc1)/cv12/cv12 <1.e-5){ + printf(" Error sqrt of a negative number: %lf\n",1+(v1-lc1)*(v1-lc1)/cv12/cv12); + fprintf(ficlog," Error sqrt of a negative number: %lf\n",1+(v1-lc1)*(v1-lc1)/cv12/cv12); + v11=(1./sqrt(fabs(1+(v1-lc1)*(v1-lc1)/cv12/cv12))); + }else + v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12)); /*v21=sqrt(1.-v11*v11); *//* error */ v21=(lc1-v1)/cv12*v11; v12=-v21; @@ -6641,8 +6683,8 @@ To be simple, these graphs help to under fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not", \ - mu1,std,v11,sqrt(lc1),v12,sqrt(fabs(lc2)), \ - mu2,std,v21,sqrt(lc1),v22,sqrt(fabs(lc2))); /* For gnuplot only */ + mu1,std,v11,sqrt(fabs(lc1)),v12,sqrt(fabs(lc2)), \ + mu2,std,v21,sqrt(fabs(lc1)),v22,sqrt(fabs(lc2))); /* For gnuplot only */ }else{ first=0; fprintf(fichtmcov," %d (%.3f),",(int) age, c12); @@ -7945,8 +7987,8 @@ set ter svg size 640, 480\nunset log y\n double *agemingoodr, *agemaxgoodr; - /* modcovmax=2*cptcoveff;/\* Max number of modalities. We suppose */ - /* a covariate has 2 modalities, should be equal to ncovcombmax *\/ */ + /* modcovmax=2*cptcoveff; Max number of modalities. We suppose */ + /* a covariate has 2 modalities, should be equal to ncovcombmax */ sumnewp = vector(1,ncovcombmax); sumnewm = vector(1,ncovcombmax); @@ -8270,11 +8312,11 @@ set ter svg size 640, 480\nunset log y\n for(j=1; j<=nlstate+ndeath;j++) { ppij=0.; for(i=1; i<=nlstate;i++) { - /* if (mobilav>=1) */ - ppij=ppij+p3mat[i][j][h]*prev[(int)agec][i][k]; - /* else { */ /* even if mobilav==-1 we use mobaverage */ - /* ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */ - /* } */ + if (mobilav>=1) + ppij=ppij+p3mat[i][j][h]*prev[(int)agec][i][k]; + else { /* even if mobilav==-1 we use mobaverage, probs may not sums to 1 */ + ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; + } fprintf(ficresf," %.3f", p3mat[i][j][h]); } /* end i */ fprintf(ficresf," %.3f", ppij); @@ -10842,7 +10884,20 @@ int main(int argc, char *argv[]) fflush(ficlog); goto end; } + /*-------- Rewriting parameter file ----------*/ + strcpy(rfileres,"r"); /* "Rparameterfile */ + strcat(rfileres,optionfilefiname); /* Parameter file first name */ + strcat(rfileres,"."); /* */ + strcat(rfileres,optionfilext); /* Other files have txt extension */ + if((ficres =fopen(rfileres,"w"))==NULL) { + printf("Problem writing new parameter file: %s\n", rfileres);goto end; + fprintf(ficlog,"Problem writing new parameter file: %s\n", rfileres);goto end; + fflush(ficlog); + goto end; + } + fprintf(ficres,"#IMaCh %s\n",version); + /* Reads comments: lines beginning with '#' */ numlinepar=0; /* Is it a BOM UTF-8 Windows file? */ @@ -10874,6 +10929,7 @@ int main(int argc, char *argv[]) numlinepar++; fputs(line,stdout); fputs(line,ficparo); + fputs(line,ficres); fputs(line,ficlog); continue; }else @@ -10894,6 +10950,7 @@ int main(int argc, char *argv[]) numlinepar++; fputs(line,stdout); fputs(line,ficparo); + fputs(line,ficres); fputs(line,ficlog); continue; }else @@ -10916,20 +10973,16 @@ int main(int argc, char *argv[]) numlinepar++; fputs(line,stdout); fputs(line,ficparo); + fputs(line,ficres); fputs(line,ficlog); continue; }else break; } if((num_filled=sscanf(line,"model=1+age%[^.\n]", model)) !=EOF){ - if (num_filled == 0){ - printf("ERROR %d: Model should be at minimum 'model=1+age.' WITHOUT space:'%s'\n",num_filled, line); - fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age.' WITHOUT space:'%s'\n",num_filled, line); - model[0]='\0'; - goto end; - } else if (num_filled != 1){ - printf("ERROR %d: Model should be at minimum 'model=1+age.' %s\n",num_filled, line); - fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age.' %s\n",num_filled, line); + if (num_filled != 1){ + printf("ERROR %d: Model should be at minimum 'model=1+age' %s\n",num_filled, line); + fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age' %s\n",num_filled, line); model[0]='\0'; goto end; } @@ -10951,11 +11004,11 @@ int main(int argc, char *argv[]) fflush(ficlog); /* if(model[0]=='#'|| model[0]== '\0'){ */ if(model[0]=='#'){ - printf("Error in 'model' line: model should start with 'model=1+age+' and end with '.' \n \ - 'model=1+age+.' or 'model=1+age+V1.' or 'model=1+age+age*age+V1+V1*age.' or \n \ - 'model=1+age+V1+V2.' or 'model=1+age+V1+V2+V1*V2.' etc. \n"); \ + printf("Error in 'model' line: model should start with 'model=1+age+' and end without space \n \ + 'model=1+age+' or 'model=1+age+V1.' or 'model=1+age+age*age+V1+V1*age' or \n \ + 'model=1+age+V1+V2' or 'model=1+age+V1+V2+V1*V2' etc. \n"); \ if(mle != -1){ - printf("Fix the model line and run imach with mle=-1 to get a correct template of the parameter file.\n"); + printf("Fix the model line and run imach with mle=-1 to get a correct template of the parameter vectors and subdiagonal covariance matrix.\n"); exit(1); } } @@ -11177,16 +11230,6 @@ Please run with mle=-1 to get a correct fflush(ficlog); - /*-------- Rewriting parameter file ----------*/ - strcpy(rfileres,"r"); /* "Rparameterfile */ - strcat(rfileres,optionfilefiname); /* Parameter file first name*/ - strcat(rfileres,"."); /* */ - strcat(rfileres,optionfilext); /* Other files have txt extension */ - if((ficres =fopen(rfileres,"w"))==NULL) { - printf("Problem writing new parameter file: %s\n", rfileres);goto end; - fprintf(ficlog,"Problem writing new parameter file: %s\n", rfileres);goto end; - } - fprintf(ficres,"#%s\n",version); } /* End of mle != -3 */ /* Main data