--- imach/src/imach.c 2017/05/24 05:45:29 1.270 +++ imach/src/imach.c 2019/05/09 14:17:20 1.292 @@ -1,6 +1,76 @@ -/* $Id: imach.c,v 1.270 2017/05/24 05:45:29 brouard Exp $ +/* $Id: imach.c,v 1.292 2019/05/09 14:17:20 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.292 2019/05/09 14:17:20 brouard + Summary: Some updates + + Revision 1.291 2019/05/09 13:44:18 brouard + Summary: Before ncovmax + + Revision 1.290 2019/05/09 13:39:37 brouard + Summary: 0.99r18 unlimited number of individuals + + The number n which was limited to 20,000 cases is now unlimited, from firstobs to lastobs. If the number is too for the virtual memory, probably an error will occur. + + Revision 1.289 2018/12/13 09:16:26 brouard + Summary: Bug for young ages (<-30) will be in r17 + + Revision 1.288 2018/05/02 20:58:27 brouard + Summary: Some bugs fixed + + Revision 1.287 2018/05/01 17:57:25 brouard + Summary: Bug fixed by providing frequencies only for non missing covariates + + Revision 1.286 2018/04/27 14:27:04 brouard + Summary: some minor bugs + + Revision 1.285 2018/04/21 21:02:16 brouard + Summary: Some bugs fixed, valgrind tested + + Revision 1.284 2018/04/20 05:22:13 brouard + Summary: Computing mean and stdeviation of fixed quantitative variables + + Revision 1.283 2018/04/19 14:49:16 brouard + Summary: Some minor bugs fixed + + Revision 1.282 2018/02/27 22:50:02 brouard + *** empty log message *** + + Revision 1.281 2018/02/27 19:25:23 brouard + Summary: Adding second argument for quitting + + 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 + + Revision 1.276 2017/06/30 15:48:31 brouard + Summary: Graphs improvements + + Revision 1.275 2017/06/30 13:39:33 brouard + Summary: Saito's color + + Revision 1.274 2017/06/29 09:47:08 brouard + Summary: Version 0.99r14 + + Revision 1.273 2017/06/27 11:06:02 brouard + Summary: More documentation on projections + + Revision 1.272 2017/06/27 10:22:40 brouard + Summary: Color of backprojection changed from 6 to 5(yellow) + + Revision 1.271 2017/06/27 10:17:50 brouard + Summary: Some bug with rint + Revision 1.270 2017/05/24 05:45:29 brouard *** empty log message *** @@ -990,14 +1060,15 @@ typedef struct { #define NINTERVMAX 8 #define NLSTATEMAX 8 /**< Maximum number of live states (for func) */ #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */ -#define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */ +#define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */ #define codtabm(h,k) (1 & (h-1) >> (k-1))+1 /*#define decodtabm(h,k,cptcoveff)= (h <= (1<> (k-1)) & 1) +1 : -1)*/ #define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1 -#define MAXN 20000 +/*#define MAXN 20000 */ /* Should by replaced by nobs, real number of observations and unlimited */ #define YEARM 12. /**< Number of months per year */ /* #define AGESUP 130 */ -#define AGESUP 150 +/* #define AGESUP 150 */ +#define AGESUP 200 #define AGEINF 0 #define AGEMARGE 25 /* Marge for agemin and agemax for(iage=agemin-AGEMARGE; iage <= agemax+3+AGEMARGE; iage++) */ #define AGEBASE 40 @@ -1013,12 +1084,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.270 2017/05/24 05:45:29 brouard Exp $ */ +/* $Id: imach.c,v 1.292 2019/05/09 14:17:20 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.270 $ $Date: 2017/05/24 05:45:29 $"; +char copyright[]="April 2018,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.292 $ $Date: 2019/05/09 14:17:20 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -1041,6 +1112,7 @@ int nqfveff=0; /**< nqfveff Number of Qu int ntveff=0; /**< ntveff number of effective time varying variables */ int nqtveff=0; /**< ntqveff number of effective time varying quantitative variables */ int cptcov=0; /* Working variable */ +int nobs=10; /* Number of observations in the data lastobs-firstobs */ int ncovcombmax=NCOVMAX; /* Maximum calculated number of covariate combination = pow(2, cptcoveff) */ int npar=NPARMAX; int nlstate=2; /* Number of live states */ @@ -2492,15 +2564,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} */ @@ -2521,6 +2596,7 @@ void powell(double p[], double **xi, int double **newm; double agefin, delaymax=200. ; /* 100 Max number of years to converge */ int ncvloop=0; + int first=0; min=vector(1,nlstate); max=vector(1,nlstate); @@ -2617,10 +2693,14 @@ void powell(double p[], double **xi, int free_vector(meandiff,1,nlstate); return prlim; } - } /* age loop */ + } /* agefin loop */ /* After some age loop it doesn't converge */ - printf("Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\ -Earliest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); + if(!first){ + first=1; + printf("Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.d years and %d loops. Try to lower 'ftolpl'. Youngest age to start was %d=(%d-%d). Others in log file only...\n", (int)age, maxmax, ftolpl, *ncvyear, ncvloop, (int)(agefin+stepm/YEARM), (int)(age-stepm/YEARM), (int)delaymax); + } + fprintf(ficlog, "Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.d years and %d loops. Try to lower 'ftolpl'. Youngest age to start was %d=(%d-%d).\n", (int)age, maxmax, ftolpl, *ncvyear, ncvloop, (int)(agefin+stepm/YEARM), (int)(age-stepm/YEARM), (int)delaymax); + /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */ free_vector(min,1,nlstate); free_vector(max,1,nlstate); @@ -2686,7 +2766,8 @@ Earliest age to start was %d-%d=%d, ncvl /* Even if hstepm = 1, at least one multiplication by the unit matrix */ /* Start at agefin= age, computes the matrix of passage and loops decreasing agefin until convergence is reached */ /* for(agefin=age+stepm/YEARM; agefin<=age+delaymax; agefin=agefin+stepm/YEARM){ /\* A changer en age *\/ */ - for(agefin=age; agefin ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. Others in log file only...\n\ Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear); @@ -2914,7 +2995,7 @@ double **pmij(double **ps, double *cov, double **out, **pmij(); double sumnew=0.; double agefin; - double k3=0.; /* constant of the w_x diagonal matrixe (in order for B to sum to 1 even for death state) */ + double k3=0.; /* constant of the w_x diagonal matrix (in order for B to sum to 1 even for death state) */ double **dnewm, **dsavm, **doldm; double **bbmij; @@ -2933,7 +3014,7 @@ double **pmij(double **ps, double *cov, /* outputs pmmij which is a stochastic matrix in row */ /* Diag(w_x) */ - /* Problem with prevacurrent which can be zero */ + /* Rescaling the cross-sectional prevalence: Problem with prevacurrent which can be zero */ sumnew=0.; /*for (ii=1;ii<=nlstate+ndeath;ii++){*/ for (ii=1;ii<=nlstate;ii++){ /* Only on live states */ @@ -2960,10 +3041,10 @@ double **pmij(double **ps, double *cov, } /* End doldm, At the end doldm is diag[(w_i)] */ - /* left Product of this diag matrix by pmmij=Px (dnewm=dsavm*doldm) */ - bbmij=matprod2(dnewm, doldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, pmmij); /* Bug Valgrind */ + /* Left product of this diag matrix by pmmij=Px (dnewm=dsavm*doldm): diag[(w_i)*Px */ + bbmij=matprod2(dnewm, doldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, pmmij); /* was a Bug Valgrind */ - /* Diag(Sum_i w^i_x p^ij_x */ + /* Diag(Sum_i w^i_x p^ij_x, should be the prevalence at age x+stepm */ /* w1 p11 + w2 p21 only on live states N1./N..*N11/N1. + N2./N..*N21/N2.=(N11+N21)/N..=N.1/N.. */ for (j=1;j<=nlstate+ndeath;j++){ sumnew=0.; @@ -2981,7 +3062,7 @@ double **pmij(double **ps, double *cov, } /*End ii */ } /* End j, At the end dsavm is diag[1/(w_1p1i+w_2 p2i)] for ALL states even if the sum is only for live states */ - ps=matprod2(ps, dnewm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dsavm); /* Bug Valgrind */ + ps=matprod2(ps, dnewm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dsavm); /* was a Bug Valgrind */ /* ps is now diag[w_i] * Px * diag [1/(w_1p1i+w_2 p2i)] */ /* end bmij */ return ps; /*pointer is unchanged */ @@ -3243,7 +3324,7 @@ double ***hbxij(double ***po, int nhstep newm=savm; /* Covariates have to be included here again */ cov[1]=1.; - agexact=age-((h-1)*hstepm + (d))*stepm/YEARM; /* age just before transition, d or d-1? */ + agexact=age-( (h-1)*hstepm + (d) )*stepm/YEARM; /* age just before transition, d or d-1? */ /* agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /\* age just before transition *\/ */ cov[2]=agexact; if(nagesqr==1) @@ -3813,7 +3894,7 @@ return -l; /*************** function likelione ***********/ -void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpri, long *ipmx, double *sw, double *fretone, double (*funcone)(double [])) +void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpri, long *ipmx, double *sw, double *fretone, double (*func)(double [])) { /* This routine should help understanding what is done with the selection of individuals/waves and @@ -3837,7 +3918,7 @@ void likelione(FILE *ficres,double p[], fprintf(ficresilk," -2*gipw/gsw*weight*ll(total)\n"); } - *fretone=(*funcone)(p); + *fretone=(*func)(p); if(*globpri !=0){ fclose(ficresilk); if (mle ==0) @@ -3845,7 +3926,7 @@ void likelione(FILE *ficres,double p[], 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,"\n
Equation of the model: model=1+age+%s
\n",model); for (k=1; k<= nlstate ; k++) { fprintf(fichtm,"
- Probability p%dj by origin %d and destination j. Dot's sizes are related to corresponding weight: %s-p%dj.png
\ @@ -4327,7 +4408,7 @@ void freqsummary(char fileres[], double double ***freq; /* Frequencies */ double *x, *y, a=0.,b=0.,r=1., sa=0., sb=0.; /* for regression, y=b+m*x and r is the correlation coefficient */ int no=0, linreg(int ifi, int ila, int *no, const double x[], const double y[], double* a, double* b, double* r, double* sa, double * sb); - double *meanq; + double *meanq, *stdq, *idq; double **meanqt; double *pp, **prop, *posprop, *pospropt; double pos=0., posproptt=0., pospropta=0., k2, dateintsum=0,k2cpt=0; @@ -4340,6 +4421,8 @@ void freqsummary(char fileres[], double pospropt=vector(1,nlstate); /* Counting the number of transition starting from a live state */ /* prop=matrix(1,nlstate,iagemin,iagemax+3); */ meanq=vector(1,nqfveff); /* Number of Quantitative Fixed Variables Effective */ + stdq=vector(1,nqfveff); /* Number of Quantitative Fixed Variables Effective */ + idq=vector(1,nqfveff); /* Number of Quantitative Fixed Variables Effective */ meanqt=matrix(1,lastpass,1,nqtveff); strcpy(fileresp,"P_"); strcat(fileresp,fileresu); @@ -4448,13 +4531,16 @@ Title=%s
Datafile=%s Firstpass=%d La posprop[i]=0; pospropt[i]=0; } - /* for (z1=1; z1<= nqfveff; z1++) { */ - /* meanq[z1]+=0.; */ + for (z1=1; z1<= nqfveff; z1++) { /* zeroing for each combination j1 as well as for the total */ + idq[z1]=0.; + meanq[z1]=0.; + stdq[z1]=0.; + } + /* for (z1=1; z1<= nqtveff; z1++) { */ /* for(m=1;m<=lastpass;m++){ */ - /* meanqt[m][z1]=0.; */ - /* } */ - /* } */ - + /* meanqt[m][z1]=0.; */ + /* } */ + /* } */ /* dateintsum=0; */ /* k2cpt=0; */ @@ -4464,9 +4550,6 @@ Title=%s
Datafile=%s Firstpass=%d La if(j !=0){ if(anyvaryingduminmodel==0){ /* If All fixed covariates */ if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ - /* for (z1=1; z1<= nqfveff; z1++) { */ - /* meanq[z1]+=coqvar[Tvar[z1]][iind]; /\* Computes mean of quantitative with selected filter *\/ */ - /* } */ for (z1=1; z1<=cptcoveff; z1++) { /* loops on covariates in the model */ /* if(Tvaraff[z1] ==-20){ */ /* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */ @@ -4487,7 +4570,7 @@ Title=%s
Datafile=%s Firstpass=%d La }/* end j==0 */ if (bool==1){ /* We selected an individual iind satisfying combination j1 (V4=1 V3=0) or all fixed covariates */ /* for(m=firstpass; m<=lastpass; m++){ */ - for(mi=1; miDatafile=%s Firstpass=%d La }/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop */ } /* end j==0 */ /* bool =0 we keep that guy which corresponds to the combination of dummy values */ - if(bool==1){ + if(bool==1){ /*Selected */ /* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind] and mw[mi+1][iind]. dh depends on stepm. */ agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/ @@ -4525,6 +4608,11 @@ Title=%s
Datafile=%s Firstpass=%d La if(s[m][iind]==-1) printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.)); freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */ + for (z1=1; z1<= nqfveff; z1++) { /* Quantitative variables, calculating mean */ + idq[z1]=idq[z1]+weight[iind]; + meanq[z1]+=covar[ncovcol+z1][iind]*weight[iind]; /* Computes mean of quantitative with selected filter */ + stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][iind]*weight[iind]*weight[iind]; /* *weight[iind];*/ /* Computes mean of quantitative with selected filter */ + } /* if((int)agev[m][iind] == 55) */ /* printf("j=%d, j1=%d Age %d, iind=%d, num=%09ld m=%d\n",j,j1,(int)agev[m][iind],iind, num[iind],m); */ /* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */ @@ -4540,6 +4628,11 @@ Title=%s
Datafile=%s Firstpass=%d La bool=1; }/* end bool 2 */ } /* end m */ + /* for (z1=1; z1<= nqfveff; z1++) { /\* Quantitative variables, calculating mean *\/ */ + /* idq[z1]=idq[z1]+weight[iind]; */ + /* meanq[z1]+=covar[ncovcol+z1][iind]*weight[iind]; /\* Computes mean of quantitative with selected filter *\/ */ + /* stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][iind]*weight[iind]*weight[iind]; /\* *weight[iind];*\/ /\* Computes mean of quantitative with selected filter *\/ */ + /* } */ } /* end bool */ } /* end iind = 1 to imx */ /* prop[s][age] is feeded for any initial and valid live state as well as @@ -4577,6 +4670,27 @@ Title=%s
Datafile=%s Firstpass=%d La fprintf(ficresphtmfr, "**********\n"); fprintf(ficlog, "**********\n"); } + /* + Printing means of quantitative variables if any + */ + for (z1=1; z1<= nqfveff; z1++) { + fprintf(ficlog,"Mean of fixed quantitative variable V%d on %.0f individuals sum=%f", ncovcol+z1, idq[z1], meanq[z1]); + fprintf(ficlog,", mean=%.3g\n",meanq[z1]/idq[z1]); + if(weightopt==1){ + printf(" Weighted mean and standard deviation of"); + fprintf(ficlog," Weighted mean and standard deviation of"); + fprintf(ficresphtmfr," Weighted mean and standard deviation of"); + } + printf(" fixed quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1])); + fprintf(ficlog," fixed quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1])); + fprintf(ficresphtmfr," fixed quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)

\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1])); + } + /* for (z1=1; z1<= nqtveff; z1++) { */ + /* for(m=1;m<=lastpass;m++){ */ + /* fprintf(ficresphtmfr,"V quantitative id %d, pass id=%d, mean=%f

\n", z1, m, meanqt[m][z1]); */ + /* } */ + /* } */ + fprintf(ficresphtm,""); if((cptcoveff==0 && nj==1)|| nj==2 ) /* no covariate and first pass */ fprintf(ficresp, " Age"); @@ -4811,7 +4925,7 @@ Title=%s
Datafile=%s Firstpass=%d La fprintf(ficlog,"\n"); } } - } + } /* end of state i */ printf("#Freqsummary\n"); fprintf(ficlog,"\n"); for(s1=-1; s1 <=nlstate+ndeath; s1++){ @@ -4857,7 +4971,9 @@ Title=%s
Datafile=%s Firstpass=%d La fclose(ficresp); fclose(ficresphtm); fclose(ficresphtmfr); + free_vector(idq,1,nqfveff); free_vector(meanq,1,nqfveff); + free_vector(stdq,1,nqfveff); free_matrix(meanqt,1,lastpass,1,nqtveff); free_vector(x, iagemin-AGEMARGE, iagemax+4+AGEMARGE); free_vector(y, iagemin-AGEMARGE, iagemax+4+AGEMARGE); @@ -4964,7 +5080,7 @@ void prevalence(double ***probs, double /*j=cptcoveff;*/ if (cptcovn<1) {j=1;ncodemax[1]=1;} - first=1; + first=0; for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){ /* For each combination of covariate */ for (i=1; i<=nlstate; i++) for(iage=iagemin-AGEMARGE; iage <= iagemax+4+AGEMARGE; iage++) @@ -5022,12 +5138,11 @@ void prevalence(double ***probs, double if(posprop>=1.e-5){ probs[i][jk][j1]= prop[jk][i]/posprop; } else{ - if(first==1){ - first=0; + if(!first){ + first=1; printf("Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,jk, j1,probs[i][jk][j1]); - fprintf(ficlog,"Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,jk, j1,probs[i][jk][j1]); }else{ - fprintf(ficlog,"Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,jk, j1,probs[i][jk][j1]); + fprintf(ficlog,"Warning Observed prevalence doesn't sum to 1 for state %d: probs[%d][%d][%d]=%lf because of lack of cases.\n",jk,i,jk, j1,probs[i][jk][j1]); } } } @@ -5273,9 +5388,11 @@ void concatwav(int wav[], int **dh, int /* *cptcov=0; */ for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */ + for (k=1; k <= maxncov; k++) + for(j=1; j<=2; j++) + nbcode[k][j]=0; /* Valgrind */ /* Loop on covariates without age and products and no quantitative variable */ - /* for (j=1; j<=(cptcovs); j++) { /\* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only *\/ */ for (k=1; k<=cptcovt; k++) { /* From model V1 + V2*age + V3 + V3*V4 keeps V1 + V3 = 2 only */ for (j=-1; (j < maxncov); j++) Ndum[j]=0; if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ @@ -5293,7 +5410,11 @@ void concatwav(int wav[], int **dh, int modmaxcovj=ij; else if (ij < modmincovj) modmincovj=ij; - if ((ij < -1) && (ij > NCOVMAX)){ + if (ij <0 || ij >1 ){ + printf("Information, IMaCh doesn't treat covariate with missing values (-1), individual %d will be skipped.\n",i); + fprintf(ficlog,"Information, currently IMaCh doesn't treat covariate with missing values (-1), individual %d will be skipped.\n",i); + } + if ((ij < -1) || (ij > NCOVMAX)){ printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX ); exit(1); }else @@ -5339,12 +5460,18 @@ void concatwav(int wav[], int **dh, int /* nbcode[Tvar[j]][3]=2; */ /* To be continued (not working yet). */ ij=0; /* ij is similar to i but can jump over null modalities */ - for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/ + + /* for (i=modmincovj; i<=modmaxcovj; i++) { */ /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/ + /* Skipping the case of missing values by reducing nbcode to 0 and 1 and not -1, 0, 1 */ + /* model=V1+V2+V3, if V2=-1, 0 or 1, then nbcode[2][1]=0 and nbcode[2][2]=1 instead of + * nbcode[2][1]=-1, nbcode[2][2]=0 and nbcode[2][3]=1 */ + /*, could be restored in the future */ + for (i=0; i<=1; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/ if (Ndum[i] == 0) { /* If nobody responded to this modality k */ break; } ij++; - nbcode[Tvar[k]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality. nbcode[1][1]=0 nbcode[1][2]=1*/ + nbcode[Tvar[k]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality. nbcode[1][1]=0 nbcode[1][2]=1 . Could be -1*/ cptcode = ij; /* New max modality for covar j */ } /* end of loop on modality i=-1 to 1 or more */ break; @@ -5360,21 +5487,7 @@ void concatwav(int wav[], int **dh, int break; } /* end switch */ } /* end dummy test */ - - /* for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 *\//\* Could be 1 to 4 *\//\* cptcode=modmaxcovj *\/ */ - /* /\*recode from 0 *\/ */ - /* k is a modality. If we have model=V1+V1*sex */ - /* then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */ - /* But if some modality were not used, it is recoded from 0 to a newer modmaxcovj=cptcode *\/ */ - /* } */ - /* /\* cptcode = ij; *\/ /\* New max modality for covar j *\/ */ - /* if (ij > ncodemax[j]) { */ - /* printf( " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */ - /* fprintf(ficlog, " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */ - /* break; */ - /* } */ - /* } /\* end of loop on modality k *\/ */ - } /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/ + } /* end of loop on model-covariate k. nbcode[Tvark][1]=-1, nbcode[Tvark][1]=0 and nbcode[Tvark][2]=1 sets the value of covariate k*/ for (k=-1; k< maxncov; k++) Ndum[k]=0; /* Look at fixed dummy (single or product) covariates to check empty modalities */ @@ -5748,22 +5861,24 @@ 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; double **dnewmp,**doldmp; int i, j, nhstepm, hstepm, h, nstepm ; + int first=0; 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; */ @@ -5824,7 +5939,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 "); @@ -5879,9 +5994,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++) @@ -5891,27 +6009,32 @@ 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); - + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij, nres); if (popbased==1) { @@ -5940,19 +6063,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 */ @@ -5963,13 +6090,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); @@ -5980,7 +6113,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++) @@ -6083,7 +6220,7 @@ void concatwav(int wav[], int **dh, int int theta; pstamp(ficresvpl); - fprintf(ficresvpl,"# Standard deviation of period (stable) prevalences \n"); + fprintf(ficresvpl,"# Standard deviation of period (forward stable) prevalences \n"); fprintf(ficresvpl,"# Age "); if(nresult >=1) fprintf(ficresvpl," Result# "); @@ -6112,20 +6249,20 @@ void concatwav(int wav[], int **dh, int for(i=1; i<=npar; i++){ /* Computes gradient */ xp[i] = x[i] + (i==theta ?delti[theta]:0); } - if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); - else - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); + /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) */ + /* prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); */ + /* else */ + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); for(i=1;i<=nlstate;i++){ gp[i] = prlim[i][i]; mgp[theta][i] = prlim[i][i]; } for(i=1; i<=npar; i++) /* Computes gradient */ xp[i] = x[i] - (i==theta ?delti[theta]:0); - if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); - else - prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); + /* if((int)age==79 ||(int)age== 80 ||(int)age== 81 ) */ + /* prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); */ + /* else */ + prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij,nres); for(i=1;i<=nlstate;i++){ gm[i] = prlim[i][i]; mgm[theta][i] = prlim[i][i]; @@ -6174,8 +6311,11 @@ void concatwav(int wav[], int **dh, int fprintf(ficresvpl,"%.0f ",age ); if(nresult >=1) fprintf(ficresvpl,"%d ",nres ); - for(i=1; i<=nlstate;i++) + for(i=1; i<=nlstate;i++){ fprintf(ficresvpl," %.5f (%.5f)",prlim[i][i],sqrt(varpl[i][(int)age])); + /* for(j=1;j<=nlstate;j++) */ + /* fprintf(ficresvpl," %d %.5f ",j,prlim[j][i]); */ + } fprintf(ficresvpl,"\n"); free_vector(gp,1,nlstate); free_vector(gm,1,nlstate); @@ -6392,7 +6532,7 @@ void varprob(char optionfilefiname[], do fprintf(fichtm,"\n

  • Computing and drawing one step probabilities with their confidence intervals

    • \n"); fprintf(fichtm,"\n"); - fprintf(fichtm,"\n

  • Matrix of variance-covariance of one-step probabilities (drawings)

    this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back. %s
    • \n",optionfilehtmcov,optionfilehtmcov); + fprintf(fichtm,"\n

  • Matrix of variance-covariance of one-step probabilities (drawings)

    this page is important in order to visualize confidence intervals and especially correlation between disability and recovery, or more generally, way in and way back. File %s
    • \n",optionfilehtmcov,optionfilehtmcov); fprintf(fichtmcov,"Current page is file %s
    \n\n

    Matrix of variance-covariance of pairs of step probabilities

    \n",optionfilehtmcov, optionfilehtmcov); fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (pij, pkl) are estimated \ and drawn. It helps understanding how is the covariance between two incidences.\ @@ -6589,7 +6729,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; @@ -6620,8 +6765,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); @@ -6658,8 +6803,8 @@ void printinghtml(char fileresu[], char int lastpass, int stepm, int weightopt, char model[],\ int imx,int jmin, int jmax, double jmeanint,char rfileres[],\ int popforecast, int mobilav, int prevfcast, int mobilavproj, int backcast, int estepm , \ - double jprev1, double mprev1,double anprev1, double dateprev1, \ - double jprev2, double mprev2,double anprev2, double dateprev2){ + double jprev1, double mprev1,double anprev1, double dateprev1, double dateproj1, double dateback1, \ + double jprev2, double mprev2,double anprev2, double dateprev2, double dateproj2, double dateback2){ int jj1, k1, i1, cpt, k4, nres; fprintf(fichtm,"
    • Result files (first order: no variance)\n \ @@ -6680,10 +6825,10 @@ void printinghtml(char fileresu[], char - Estimated back transition probabilities over %d (stepm) months: %s
      \n ", stepm,subdirf2(fileresu,"PIJB_"),subdirf2(fileresu,"PIJB_")); fprintf(fichtm,"\ - - Period (stable) prevalence in each health state: %s
      \n", + - Period (forward) prevalence in each health state: %s
      \n", subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_")); fprintf(fichtm,"\ - - Period (stable) back prevalence in each health state: %s
      \n", + - Backward prevalence in each health state: %s
      \n", subdirf2(fileresu,"PLB_"),subdirf2(fileresu,"PLB_")); fprintf(fichtm,"\ - (a) Life expectancies by health status at initial age, ei. (b) health expectancies by health status at initial age, eij . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \ @@ -6789,39 +6934,42 @@ divided by h: hPij ",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres); /* Survival functions (period) in state j */ for(cpt=1; cpt<=nlstate;cpt++){ - fprintf(fichtm,"
      \n- Survival functions in state %d. Or probability to survive in state %d being in state (1 to %d) at different ages. %s_%d-%d-%d.svg
      \ + fprintf(fichtm,"
      \n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages. %s_%d-%d-%d.svg
      \ ", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres); } /* State specific survival functions (period) */ for(cpt=1; cpt<=nlstate;cpt++){ - fprintf(fichtm,"
      \n- Survival functions from state %d in each live state and total.\ - Or probability to survive in various states (1 to %d) being in state %d at different ages. \ - %s_%d%d-%d.svg
      ", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres); + fprintf(fichtm,"
      \n- Survival functions in state %d and in any other live state (total).\ + And probability to be observed in various states (up to %d) being in state %d at different ages. \ + %s_%d-%d-%d.svg
      ", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres); } - /* Period (stable) prevalence in each health state */ + /* Period (forward stable) prevalence in each health state */ for(cpt=1; cpt<=nlstate;cpt++){ fprintf(fichtm,"
      \n- Convergence to period (stable) prevalence in state %d. Or probability for a person being in state (1 to %d) at different ages, to be in state %d some years after. %s_%d-%d-%d.svg
      \ ", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres); } if(backcast==1){ - /* Period (stable) back prevalence in each health state */ + /* Backward prevalence in each health state */ for(cpt=1; cpt<=nlstate;cpt++){ fprintf(fichtm,"
      \n- Convergence to mixed (stable) back prevalence in state %d. Or probability for a person to be in state %d at a younger age, knowing that she/he was in state (1 to %d) at different older ages. %s_%d-%d-%d.svg
      \ ", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres); } } if(prevfcast==1){ - /* Projection of prevalence up to period (stable) prevalence in each health state */ + /* Projection of prevalence up to period (forward stable) prevalence in each health state */ for(cpt=1; cpt<=nlstate;cpt++){ - fprintf(fichtm,"
      \n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d) up to period (stable) prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). %s_%d-%d-%d.svg
      \ -", dateprev1, dateprev2, mobilavproj, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres); + fprintf(fichtm,"
      \n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), from year %.1f up to year %.1f tending to period (stable) forward prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). %s_%d-%d-%d.svg
      \ +", dateprev1, dateprev2, mobilavproj, dateproj1, dateproj2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres); } } if(backcast==1){ /* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */ for(cpt=1; cpt<=nlstate;cpt++){ - fprintf(fichtm,"
      \n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d) up to stable (mixed) back prevalence in state %d. Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) with weights corresponding to observed prevalence at different ages. %s_%d-%d-%d.svg
      \ -", dateprev1, dateprev2, mobilavproj, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres); + fprintf(fichtm,"
      \n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), \ + from year %.1f up to year %.1f (probably close to stable [mixed] back prevalence in state %d (randomness in cross-sectional prevalence is not taken into \ + account but can visually be appreciated). Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) \ +with weights corresponding to observed prevalence at different ages. %s_%d-%d-%d.svg
      \ + ", dateprev1, dateprev2, mobilavproj, dateback1, dateback2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres); } } @@ -6863,13 +7011,13 @@ See page 'Matrix of variance-covariance %s
      \n
      • ", estepm,subdirf2(fileresu,"STDE_"),subdirf2(fileresu,"STDE_")); fprintf(fichtm,"\ - - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), eij are weighted by the period prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): %s
      \n", + - Variances and covariances of health expectancies by age. Status (i) based health expectancies (in state j), eij are weighted by the forward (period) prevalences in each state i (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): %s
      \n", estepm, subdirf2(fileresu,"V_"),subdirf2(fileresu,"V_")); fprintf(fichtm,"\ - Total life expectancy and total health expectancies to be spent in each health state e.j with their standard errors (if popbased=1, an additional computation is done using the cross-sectional prevalences, i.e population based) (estepm=%d months): %s
      \n", estepm, subdirf2(fileresu,"T_"),subdirf2(fileresu,"T_")); fprintf(fichtm,"\ - - Standard deviation of period (stable) prevalences: %s
      \n",\ + - Standard deviation of forward (period) prevalences: %s
      \n",\ subdirf2(fileresu,"VPL_"),subdirf2(fileresu,"VPL_")); /* if(popforecast==1) fprintf(fichtm,"\n */ @@ -6951,6 +7099,20 @@ void printinggnuplot(char fileresu[], ch /*#endif */ m=pow(2,cptcoveff); + /* diagram of the model */ + fprintf(ficgp,"\n#Diagram of the model \n"); + fprintf(ficgp,"\ndelta=0.03;delta2=0.07;unset arrow;\n"); + fprintf(ficgp,"yoff=(%d > 2? 0:1);\n",nlstate); + fprintf(ficgp,"\n#Peripheral arrows\nset for [i=1:%d] for [j=1:%d] arrow i*10+j from cos(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d))-(i!=j?(i-j)/abs(i-j)*delta:0), yoff +sin(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) rto -0.95*(cos(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d))+(i!=j?(i-j)/abs(i-j)*delta:0) - cos(pi*((1-(%d/2)*2./%d)/2+(j-1)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta2:0)), -0.95*(sin(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) - sin(pi*((1-(%d/2)*2./%d)/2+(j-1)*2./%d))+( i!=j?(i-j)/abs(i-j)*delta2:0)) ls (i < j? 1:2)\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate); + + fprintf(ficgp,"\n#Centripete arrows (turning in other direction (1-i) instead of (i-1)) \nset for [i=1:%d] arrow (%d+1)*10+i from cos(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d))-(i!=j?(i-j)/abs(i-j)*delta:0), yoff +sin(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) rto -0.80*(cos(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d))+(i!=j?(i-j)/abs(i-j)*delta:0) ), -0.80*(sin(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) + yoff ) ls 4\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate); + fprintf(ficgp,"\n#show arrow\nunset label\n"); + fprintf(ficgp,"\n#States labels, starting from 2 (2-i) instead of (1-i), was (i-1)\nset for [i=1:%d] label i sprintf(\"State %%d\",i) center at cos(pi*((1-(%d/2)*2./%d)/2+(2-i)*2./%d)), yoff+sin(pi*((1-(%d/2)*2./%d)/2+(2-i)*2./%d)) font \"helvetica, 16\" tc rgbcolor \"blue\"\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate); + fprintf(ficgp,"\nset label %d+1 sprintf(\"State %%d\",%d+1) center at 0.,0. font \"helvetica, 16\" tc rgbcolor \"red\"\n",nlstate,nlstate); + fprintf(ficgp,"\n#show label\nunset border;unset xtics; unset ytics;\n"); + fprintf(ficgp,"\n\nset ter svg size 640, 480;set out \"%s_.svg\" \n",subdirf2(optionfilefiname,"D_")); + fprintf(ficgp,"unset log y; plot [-1.2:1.2][yoff-1.2:1.2] 1/0 not; set out;reset;\n"); + /* Contribution to likelihood */ /* Plot the probability implied in the likelihood */ fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n"); @@ -6991,7 +7153,7 @@ void printinggnuplot(char fileresu[], ch continue; /* We are interested in selected combination by the resultline */ /* printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); */ - fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); + fprintf(ficgp,"\n# 1st: Forward (stable period) prevalence with CI: 'VPL_' files and live state =%d ", cpt); strcpy(gplotlabel,"("); for (k=1; k<=cptcoveff; k++){ /* For each covariate k get corresponding value lv for combination k1 */ lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */ @@ -7020,7 +7182,8 @@ void printinggnuplot(char fileresu[], ch fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1,nres); fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres); - fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel); + /* fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel); */ + fprintf(ficgp,"set title \"Alive state %d %s\" font \"Helvetica,12\"\n",cpt,gplotlabel); fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres); /* fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); */ /* k1-1 error should be nres-1*/ @@ -7028,7 +7191,7 @@ void printinggnuplot(char fileresu[], ch if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); else fprintf(ficgp," %%*lf (%%*lf)"); } - fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres); + fprintf(ficgp,"\" t\"Forward prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres); for (i=1; i<= nlstate ; i ++) { if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); else fprintf(ficgp," %%*lf (%%*lf)"); @@ -7042,7 +7205,7 @@ void printinggnuplot(char fileresu[], ch fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" u 1:((",subdirf2(fileresu,"P_")); if(cptcoveff ==0){ - fprintf(ficgp,"$%d)) t 'Observed prevalence in state %d' with line lt 3", 2+(cpt-1), cpt ); + fprintf(ficgp,"$%d)) t 'Observed prevalence in state %d' with line lt 3", 2+3*(cpt-1), cpt ); }else{ kl=0; for (k=1; k<=cptcoveff; k++){ /* For each combination of covariate */ @@ -7100,20 +7263,21 @@ void printinggnuplot(char fileresu[], ch if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); else fprintf(ficgp," %%*lf (%%*lf)"); } - fprintf(ficgp,"\" t\"Backward (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres); + fprintf(ficgp,"\" t\"Backward (stable) prevalence\" w l lt 6 dt 3,\"%s\" every :::%d::%d u 1:($2==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres); for (i=1; i<= nlstate ; i ++) { if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); else fprintf(ficgp," %%*lf (%%*lf)"); } - fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres); + fprintf(ficgp,"\" t\"95%% CI\" w l lt 4,\"%s\" every :::%d::%d u 1:($2==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres); for (i=1; i<= nlstate ; i ++) { if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); else fprintf(ficgp," %%*lf (%%*lf)"); } - fprintf(ficgp,"\" t\"\" w l lt 1"); + fprintf(ficgp,"\" t\"\" w l lt 4"); } /* end if backprojcast */ } /* end if backcast */ - fprintf(ficgp,"\nset out ;unset label;\n"); + /* fprintf(ficgp,"\nset out ;unset label;\n"); */ + fprintf(ficgp,"\nset out ;unset title;\n"); } /* nres */ } /* k1 */ } /* cpt */ @@ -7356,7 +7520,7 @@ set ter svg size 640, 480\nunset log y\n continue; for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state of arrival */ strcpy(gplotlabel,"("); - fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt); + fprintf(ficgp,"\n#\n#\n#CV preval stable (forward): 'pij' files, covariatecombination#=%d state=%d",k1, cpt); for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ @@ -7400,14 +7564,14 @@ set ter svg size 640, 480\nunset log y\n /* 7eme */ if(backcast == 1){ - /* CV back preval stable (period) for each covariate */ + /* CV backward prevalence for each covariate */ for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */ for(nres=1; nres <= nresult; nres++){ /* For each resultline */ if(m != 1 && TKresult[nres]!= k1) continue; for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life origin state */ strcpy(gplotlabel,"("); - fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pijb' files, covariatecombination#=%d state=%d",k1, cpt); + fprintf(ficgp,"\n#\n#\n#CV Backward stable prevalence: 'pijb' files, covariatecombination#=%d state=%d",k1, cpt); for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ @@ -7455,7 +7619,7 @@ set ter svg size 640, 480\nunset log y\n /* 8eme */ if(prevfcast==1){ - /* Projection from cross-sectional to stable (period) for each covariate */ + /* Projection from cross-sectional to forward stable (period) prevalence for each covariate */ for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */ for(nres=1; nres <= nresult; nres++){ /* For each resultline */ @@ -7463,7 +7627,7 @@ set ter svg size 640, 480\nunset log y\n continue; for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ strcpy(gplotlabel,"("); - fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt); + fprintf(ficgp,"\n#\n#\n#Projection of prevalence to forward stable prevalence (period): 'PROJ_' files, covariatecombination#=%d state=%d",k1, cpt); for (k=1; k<=cptcoveff; k++){ /* For each correspondig covariate value */ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ @@ -7727,7 +7891,7 @@ set ter svg size 640, 480\nunset log y\n continue; fprintf(ficgp,"\n\n# Combination of dummy k1=%d which is ",k1); strcpy(gplotlabel,"("); - sprintf(gplotlabel+strlen(gplotlabel)," Dummy combination %d ",k1); + /*sprintf(gplotlabel+strlen(gplotlabel)," Dummy combination %d ",k1);*/ for (k=1; k<=cptcoveff; k++){ /* For each correspondig covariate value */ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ @@ -7744,7 +7908,9 @@ set ter svg size 640, 480\nunset log y\n strcpy(gplotlabel+strlen(gplotlabel),")"); fprintf(ficgp,"\n#\n"); fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),k1,ng,nres); - fprintf(ficgp,"\nset label \"%s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",gplotlabel); + fprintf(ficgp,"\nset key outside "); + /* fprintf(ficgp,"\nset label \"%s\" at graph 1.2,0.5 center rotate font \"Helvetica,12\"\n",gplotlabel); */ + fprintf(ficgp,"\nset title \"%s\" font \"Helvetica,12\"\n",gplotlabel); fprintf(ficgp,"\nset ter svg size 640, 480 "); if (ng==1){ fprintf(ficgp,"\nset ylabel \"Value of the logit of the model\"\n"); /* exp(a12+b12*x) could be nice */ @@ -7864,12 +8030,12 @@ set ter svg size 640, 480\nunset log y\n } fprintf(ficgp,")"); if(ng ==2) - fprintf(ficgp," t \"p%d%d\" ", k2,k); + fprintf(ficgp," w l lw 2 lt (%d*%d+%d)%%%d+1 dt %d t \"p%d%d\" ", nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k); else /* ng= 3 */ - fprintf(ficgp," t \"i%d%d\" ", k2,k); + fprintf(ficgp," w l lw 2 lt (%d*%d+%d)%%%d+1 dt %d t \"i%d%d\" ", nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k); }else{ /* end ng <> 1 */ if( k !=k2) /* logit p11 is hard to draw */ - fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k); + fprintf(ficgp," w l lw 2 lt (%d*%d+%d)%%%d+1 dt %d t \"logit(p%d%d)\" ", nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k); } if ((k+k2)!= (nlstate*2+ndeath) && ng != 1) fprintf(ficgp,","); @@ -7878,7 +8044,8 @@ set ter svg size 640, 480\nunset log y\n i=i+ncovmodel; } /* end k */ } /* end k2 */ - fprintf(ficgp,"\n set out; unset label;\n"); + /* fprintf(ficgp,"\n set out; unset label;set key default;\n"); */ + fprintf(ficgp,"\n set out; unset title;set key default;\n"); } /* end k1 */ } /* end ng */ /* avoid: */ @@ -7894,6 +8061,7 @@ set ter svg size 640, 480\nunset log y\n int modcovmax =1; int mobilavrange, mob; int iage=0; + int firstA1=0, firstA2=0; double sum=0., sumr=0.; double age; @@ -7902,8 +8070,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); @@ -7991,6 +8159,7 @@ set ter svg size 640, 480\nunset log y\n } /* age */ /* Thus we have agemingood and agemaxgood as well as goodr for raw (preobs) */ /* but they will change */ + firstA1=0;firstA2=0; for (age=fage-(mob-1)/2; age>=bage; age--){/* From oldest to youngest, filling up to the youngest */ sumnewm[cptcod]=0.; sumnewmr[cptcod]=0.; @@ -8023,11 +8192,19 @@ set ter svg size 640, 480\nunset log y\n sumr+=probs[(int)age][i][cptcod]; } if(fabs(sum - 1.) > 1.e-3) { /* bad */ - printf("Moving average A1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, (int)bage); + if(!firstA1){ + firstA1=1; + printf("Moving average A1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you increase bage=%d. Others in log file...\n",cptcod,sumr, (int)age, (int)bage); + } + fprintf(ficlog,"Moving average A1: For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one (%f) at any descending age! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, (int)bage); } /* end bad */ /* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */ if(fabs(sumr - 1.) > 1.e-3) { /* bad */ - printf("Moving average A2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, (int)bage); + if(!firstA2){ + firstA2=1; + printf("Moving average A2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase bage=%d. Others in log file...\n",cptcod,sumr, (int)age, (int)bage); + } + fprintf(ficlog,"Moving average A2: For this combination of covariate cptcod=%d, the raw prevalence doesn't sums to one (%f) even with smoothed values at young ages! age=%d, could you increase bage=%d\n",cptcod,sumr, (int)age, (int)bage); } /* end bad */ }/* age */ @@ -8227,11 +8404,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); @@ -8349,11 +8526,12 @@ set ter svg size 640, 480\nunset log y\n /* for (agec=bage; agec<=agemax-1; agec++){ /\* testing *\/ */ for (agec=bage; agec<=fage; agec++){ /* testing */ /* We compute bij at age agec over nhstepm, nhstepm decreases when agec increases because of agemax;*/ - nhstepm=(int) rint((agec-agelim)*YEARM/stepm); + nhstepm=(int) (agec-agelim) *YEARM/stepm;/* nhstepm=(int) rint((agec-agelim)*YEARM/stepm);*/ nhstepm = nhstepm/hstepm; p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); oldm=oldms;savm=savms; /* computes hbxij at age agec over 1 to nhstepm */ + /* printf("####prevbackforecast debug agec=%.2f nhstepm=%d\n",agec, nhstepm);fflush(stdout); */ hbxij(p3mat,nhstepm,agec,hstepm,p,prevacurrent,nlstate,stepm, k, nres); /* hpxij(p3mat,nhstepm,agec,hstepm,p, nlstate,stepm,oldm,savm, k,nres); */ /* Then we print p3mat for h corresponding to the right agec+h*stepms=yearp */ @@ -8401,7 +8579,7 @@ set ter svg size 640, 480\nunset log y\n /* Variance of prevalence limit: varprlim */ void varprlim(char fileresu[], int nresult, double ***prevacurrent, int mobilavproj, double bage, double fage, double **prlim, int *ncvyearp, double ftolpl, double p[], double **matcov, double *delti, int stepm, int cptcoveff){ - /*------- Variance of period (stable) prevalence------*/ + /*------- Variance of forward period (stable) prevalence------*/ char fileresvpl[FILENAMELENGTH]; FILE *ficresvpl; @@ -8412,11 +8590,11 @@ set ter svg size 640, 480\nunset log y\n strcpy(fileresvpl,"VPL_"); strcat(fileresvpl,fileresu); if((ficresvpl=fopen(fileresvpl,"w"))==NULL) { - printf("Problem with variance of period (stable) prevalence resultfile: %s\n", fileresvpl); + printf("Problem with variance of forward period (stable) prevalence resultfile: %s\n", fileresvpl); exit(0); } - printf("Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout); - fprintf(ficlog, "Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(ficlog); + printf("Computing Variance-covariance of forward period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout); + fprintf(ficlog, "Computing Variance-covariance of forward period (stable) prevalence: file '%s' ...", fileresvpl);fflush(ficlog); /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){ for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/ @@ -8453,8 +8631,8 @@ set ter svg size 640, 480\nunset log y\n } fclose(ficresvpl); - printf("done variance-covariance of period prevalence\n");fflush(stdout); - fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog); + printf("done variance-covariance of forward period prevalence\n");fflush(stdout); + fprintf(ficlog,"done variance-covariance of forward period prevalence\n");fflush(ficlog); } /* Variance of back prevalence: varbprlim */ @@ -9597,7 +9775,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative ( Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product \n\ Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\ Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model); - for(k=1;k<=cptcovt; k++){ Fixed[k]=0; Dummy[k]=0;} + for(k=-1;k<=cptcovt; k++){ Fixed[k]=0; Dummy[k]=0;} for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */ if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */ Fixed[k]= 0; @@ -9847,11 +10025,12 @@ Dummy[k] 0=dummy (0 1), 1 quantitative ( /* Searching for doublons in the model */ for(k1=1; k1<= cptcovt;k1++){ for(k2=1; k2 + + /* #include "syscompilerinfo.h"*/ /* command line Intel compiler 32bit windows, XP compatible:*/ /* /GS /W3 /Gy /Zc:wchar_t /Zi /O2 /Fd"Release\vc120.pdb" /D "WIN32" /D "NDEBUG" /D @@ -10059,6 +10240,8 @@ void syscompilerinfo(int logged) /ManifestFile:"Release\IMaCh.exe.intermediate.manifest" /OPT:ICF /NOLOGO /TLBID:1 */ + + #if defined __INTEL_COMPILER #if defined(__GNUC__) struct utsname sysInfo; /* For Intel on Linux and OS/X */ @@ -10075,8 +10258,6 @@ void syscompilerinfo(int logged) } #endif -#include - printf("Compiled with:");if(logged)fprintf(ficlog,"Compiled with:"); #if defined(__clang__) printf(" Clang/LLVM");if(logged)fprintf(ficlog," Clang/LLVM"); /* Clang/LLVM. ---------------------------------------------- */ @@ -10162,7 +10343,7 @@ void syscompilerinfo(int logged) #endif #endif - // void main() + // void main () // { #if defined(_MSC_VER) if (IsWow64()){ @@ -10183,7 +10364,7 @@ void syscompilerinfo(int logged) } int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp){ - /*--------------- Prevalence limit (period or stable prevalence) --------------*/ + /*--------------- Prevalence limit (forward period or forward stable prevalence) --------------*/ int i, j, k, i1, k4=0, nres=0 ; /* double ftolpl = 1.e-10; */ double age, agebase, agelim; @@ -10192,13 +10373,13 @@ int prevalence_limit(double *p, double * strcpy(filerespl,"PL_"); strcat(filerespl,fileresu); if((ficrespl=fopen(filerespl,"w"))==NULL) { - printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1; - fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1; + printf("Problem with forward period (stable) prevalence resultfile: %s\n", filerespl);return 1; + fprintf(ficlog,"Problem with forward period (stable) prevalence resultfile: %s\n", filerespl);return 1; } - printf("\nComputing period (stable) prevalence: result on file '%s' \n", filerespl); - fprintf(ficlog,"\nComputing period (stable) prevalence: result on file '%s' \n", filerespl); + printf("\nComputing forward period (stable) prevalence: result on file '%s' \n", filerespl); + fprintf(ficlog,"\nComputing forward period (stable) prevalence: result on file '%s' \n", filerespl); pstamp(ficrespl); - fprintf(ficrespl,"# Period (stable) prevalence. Precision given by ftolpl=%g \n", ftolpl); + fprintf(ficrespl,"# Forward period (stable) prevalence. Precision given by ftolpl=%g \n", ftolpl); fprintf(ficrespl,"#Age "); for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); fprintf(ficrespl,"\n"); @@ -10273,7 +10454,7 @@ int prevalence_limit(double *p, double * } int back_prevalence_limit(double *p, double **bprlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp, double dateprev1,double dateprev2, int firstpass, int lastpass, int mobilavproj){ - /*--------------- Back Prevalence limit (period or stable prevalence) --------------*/ + /*--------------- Back Prevalence limit (backward stable prevalence) --------------*/ /* Computes the back prevalence limit for any combination of covariate values * at any age between ageminpar and agemaxpar @@ -10288,13 +10469,13 @@ int back_prevalence_limit(double *p, dou strcpy(fileresplb,"PLB_"); strcat(fileresplb,fileresu); if((ficresplb=fopen(fileresplb,"w"))==NULL) { - printf("Problem with period (stable) back prevalence resultfile: %s\n", fileresplb);return 1; - fprintf(ficlog,"Problem with period (stable) back prevalence resultfile: %s\n", fileresplb);return 1; + printf("Problem with backward prevalence resultfile: %s\n", fileresplb);return 1; + fprintf(ficlog,"Problem with backward prevalence resultfile: %s\n", fileresplb);return 1; } - printf("Computing period (stable) back prevalence: result on file '%s' \n", fileresplb); - fprintf(ficlog,"Computing period (stable) back prevalence: result on file '%s' \n", fileresplb); + printf("Computing backward prevalence: result on file '%s' \n", fileresplb); + fprintf(ficlog,"Computing backward prevalence: result on file '%s' \n", fileresplb); pstamp(ficresplb); - fprintf(ficresplb,"# Period (stable) back prevalence. Precision given by ftolpl=%g \n", ftolpl); + fprintf(ficresplb,"# Backward prevalence. Precision given by ftolpl=%g \n", ftolpl); fprintf(ficresplb,"#Age "); for(i=1; i<=nlstate;i++) fprintf(ficresplb,"%d-%d ",i,i); fprintf(ficresplb,"\n"); @@ -10483,7 +10664,7 @@ int hPijx(double *p, int bage, int fage) /*if (stepm<=24) stepsize=2;*/ /* agelim=AGESUP; */ - ageminl=30; + ageminl=AGEINF; /* was 30 */ hstepm=stepsize*YEARM; /* Every year of age */ hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ @@ -10562,7 +10743,8 @@ int main(int argc, char *argv[]) double ssval; #endif int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav); - int i,j, k, n=MAXN,iter=0,m,size=100, cptcod; + int i,j, k, iter=0,m,size=100, cptcod; /* Suppressing because nobs */ + /* int i,j, k, n=MAXN,iter=0,m,size=100, cptcod; */ int ncvyear=0; /* Number of years needed for the period prevalence to converge */ int jj, ll, li, lj, lk; int numlinepar=0; /* Current linenumber of parameter file */ @@ -10572,7 +10754,9 @@ int main(int argc, char *argv[]) int vpopbased=0; int nres=0; int endishere=0; - + int noffset=0; + int ncurrv=0; /* Temporary variable */ + char ca[32], cb[32]; /* FILE *fichtm; *//* Html File */ /* FILE *ficgp;*/ /*Gnuplot File */ @@ -10595,7 +10779,7 @@ int main(int argc, char *argv[]) char pathr[MAXLINE], pathimach[MAXLINE]; char *tok, *val; /* pathtot */ - int firstobs=1, lastobs=10; + int firstobs=1, lastobs=10; /* nobs = lastobs-firstobs declared globally ;*/ int c, h , cpt, c2; int jl=0; int i1, j1, jk, stepsize=0; @@ -10603,7 +10787,7 @@ int main(int argc, char *argv[]) int *tab; int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */ - int backcast=0; + /* int backcast=0; */ /* defined as global for mlikeli and mle */ int mobilav=0,popforecast=0; int hstepm=0, nhstepm=0; int agemortsup; @@ -10627,8 +10811,9 @@ int main(int argc, char *argv[]) double *epj, vepp; - double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000; - double jback1=1,mback1=1,anback1=2000,jback2=1,mback2=1,anback2=2000; + double dateprev1, dateprev2; + double jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000, dateproj1=0, dateproj2=0; + double jback1=1,mback1=1,anback1=2000,jback2=1,mback2=1,anback2=2000, dateback1=0, dateback2=0; double **ximort; char *alph[]={"a","a","b","c","d","e"}, str[4]="1234"; @@ -10706,8 +10891,13 @@ int main(int argc, char *argv[]) if(pathr[0] == '\0') break; /* Dirty */ } } + else if (argc<=2){ + strcpy(pathtot,argv[1]); + } else{ strcpy(pathtot,argv[1]); + strcpy(z,argv[2]); + printf("\nargv[2]=%s z=%c\n",argv[2],z[0]); } /*if(getcwd(pathcd, MAXLINE)!= NULL)printf ("Error pathcd\n");*/ /*cygwin_split_path(pathtot,path,optionfile); @@ -10785,8 +10975,6 @@ int main(int argc, char *argv[]) exit(70); } - - strcpy(filereso,"o"); strcat(filereso,fileresu); if((ficparo=fopen(filereso,"w"))==NULL) { /* opened on subdirectory */ @@ -10795,17 +10983,52 @@ 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; - - /* First parameter line */ + /* Is it a BOM UTF-8 Windows file? */ + /* First parameter line */ while(fgets(line, MAXLINE, ficpar)) { + noffset=0; + if( line[0] == (char)0xEF && line[1] == (char)0xBB) /* EF BB BF */ + { + noffset=noffset+3; + printf("# File is an UTF8 Bom.\n"); // 0xBF + } + else if( line[0] == (char)0xFE && line[1] == (char)0xFF) + { + noffset=noffset+2; + printf("# File is an UTF16BE BOM file\n"); + } + else if( line[0] == 0 && line[1] == 0) + { + if( line[2] == (char)0xFE && line[3] == (char)0xFF){ + noffset=noffset+4; + printf("# File is an UTF16BE BOM file\n"); + } + } else{ + ;/*printf(" Not a BOM file\n");*/ + } + /* If line starts with a # it is a comment */ - if (line[0] == '#') { + if (line[noffset] == '#') { numlinepar++; fputs(line,stdout); fputs(line,ficparo); + fputs(line,ficres); fputs(line,ficlog); continue; }else @@ -10815,18 +11038,24 @@ int main(int argc, char *argv[]) title, datafile, &lastobs, &firstpass,&lastpass)) !=EOF){ if (num_filled != 5) { printf("Should be 5 parameters\n"); + fprintf(ficlog,"Should be 5 parameters\n"); } numlinepar++; printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", title, datafile, lastobs, firstpass,lastpass); + fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", title, datafile, lastobs, firstpass,lastpass); + fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", title, datafile, lastobs, firstpass,lastpass); + fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", title, datafile, lastobs, firstpass,lastpass); } /* Second parameter line */ while(fgets(line, MAXLINE, ficpar)) { - /* If line starts with a # it is a comment */ + /* while(fscanf(ficpar,"%[^\n]", line)) { */ + /* If line starts with a # it is a comment. Strangely fgets reads the EOL and fputs doesn't */ if (line[0] == '#') { numlinepar++; - fputs(line,stdout); - fputs(line,ficparo); - fputs(line,ficlog); + printf("%s",line); + fprintf(ficres,"%s",line); + fprintf(ficparo,"%s",line); + fprintf(ficlog,"%s",line); continue; }else break; @@ -10836,8 +11065,19 @@ int main(int argc, char *argv[]) if (num_filled != 11) { printf("Not 11 parameters, for example:ftol=1.e-8 stepm=12 ncovcol=2 nqv=1 ntv=2 nqtv=1 nlstate=2 ndeath=1 maxwav=3 mle=1 weight=1\n"); printf("but line=%s\n",line); + fprintf(ficlog,"Not 11 parameters, for example:ftol=1.e-8 stepm=12 ncovcol=2 nqv=1 ntv=2 nqtv=1 nlstate=2 ndeath=1 maxwav=3 mle=1 weight=1\n"); + fprintf(ficlog,"but line=%s\n",line); } - printf("ftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, mle, weightopt); + if( lastpass > maxwav){ + printf("Error (lastpass = %d) > (maxwav = %d)\n",lastpass, maxwav); + fprintf(ficlog,"Error (lastpass = %d) > (maxwav = %d)\n",lastpass, maxwav); + fflush(ficlog); + goto end; + } + printf("ftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, mle, weightopt); + fprintf(ficparo,"ftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, mle, weightopt); + fprintf(ficres,"ftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, 0, weightopt); + fprintf(ficlog,"ftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, mle, weightopt); } /* ftolpl=6*ftol*1.e5; /\* 6.e-3 make convergences in less than 80 loops for the prevalence limit *\/ */ /*ftolpl=6.e-4; *//* 6.e-3 make convergences in less than 80 loops for the prevalence limit */ @@ -10846,22 +11086,18 @@ int main(int argc, char *argv[]) /* If line starts with a # it is a comment */ if (line[0] == '#') { numlinepar++; - fputs(line,stdout); - fputs(line,ficparo); - fputs(line,ficlog); + printf("%s",line); + fprintf(ficres,"%s",line); + fprintf(ficparo,"%s",line); + fprintf(ficlog,"%s",line); 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; } @@ -10874,20 +11110,23 @@ int main(int argc, char *argv[]) } /* printf(" model=1+age%s modeltemp= %s, model=%s\n",model, modeltemp, model);fflush(stdout); */ printf("model=1+age+%s\n",model);fflush(stdout); + fprintf(ficparo,"model=1+age+%s\n",model);fflush(stdout); + fprintf(ficres,"model=1+age+%s\n",model);fflush(stdout); + fprintf(ficlog,"model=1+age+%s\n",model);fflush(stdout); } /* 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+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); */ - fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol, nqv, ntv, nqtv, nlstate,ndeath,maxwav, mle, weightopt,model); - fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol, nqv, ntv, nqtv, nlstate,ndeath,maxwav, mle, weightopt,model); + /* fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol, nqv, ntv, nqtv, nlstate,ndeath,maxwav, mle, weightopt,model); */ + /* fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol, nqv, ntv, nqtv, nlstate,ndeath,maxwav, mle, weightopt,model); */ 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); } } @@ -10907,10 +11146,10 @@ int main(int argc, char *argv[]) ungetc(c,ficpar); - covar=matrix(0,NCOVMAX,1,n); /**< used in readdata */ - if(nqv>=1)coqvar=matrix(1,nqv,1,n); /**< Fixed quantitative covariate */ - if(nqtv>=1)cotqvar=ma3x(1,maxwav,1,nqtv,1,n); /**< Time varying quantitative covariate */ - if(ntv+nqtv>=1)cotvar=ma3x(1,maxwav,1,ntv+nqtv,1,n); /**< Time varying covariate (dummy and quantitative)*/ + covar=matrix(0,NCOVMAX,firstobs,lastobs); /**< used in readdata */ + if(nqv>=1)coqvar=matrix(1,nqv,firstobs,lastobs); /**< Fixed quantitative covariate */ + if(nqtv>=1)cotqvar=ma3x(1,maxwav,1,nqtv,firstobs,lastobs); /**< Time varying quantitative covariate */ + if(ntv+nqtv>=1)cotvar=ma3x(1,maxwav,1,ntv+nqtv,firstobs,lastobs); /**< Time varying covariate (dummy and quantitative)*/ cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement plus one, indepently of n in Vn*/ /* v1+v2+v3+v2*v4+v5*age makes cptcovn = 5 v1+v2*age+v2*v3 makes cptcovn = 3 @@ -10973,6 +11212,15 @@ int main(int argc, char *argv[]) for(jj=1; jj <=nlstate+ndeath; jj++){ if(jj==i) continue; j++; + while((c=getc(ficpar))=='#' && c!= EOF){ + ungetc(c,ficpar); + fgets(line, MAXLINE, ficpar); + numlinepar++; + fputs(line,stdout); + fputs(line,ficparo); + fputs(line,ficlog); + } + ungetc(c,ficpar); fscanf(ficpar,"%1d%1d",&i1,&j1); if ((i1 != i) || (j1 != jj)){ printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n \ @@ -11109,30 +11357,29 @@ 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 */ - n= lastobs; - num=lvector(1,n); - moisnais=vector(1,n); - annais=vector(1,n); - moisdc=vector(1,n); - andc=vector(1,n); - weight=vector(1,n); - agedc=vector(1,n); - cod=ivector(1,n); - for(i=1;i<=n;i++){ + nobs=lastobs-firstobs+1; /* was = lastobs;*/ + /* num=lvector(1,n); */ + /* moisnais=vector(1,n); */ + /* annais=vector(1,n); */ + /* moisdc=vector(1,n); */ + /* andc=vector(1,n); */ + /* weight=vector(1,n); */ + /* agedc=vector(1,n); */ + /* cod=ivector(1,n); */ + /* for(i=1;i<=n;i++){ */ + num=lvector(firstobs,lastobs); + moisnais=vector(firstobs,lastobs); + annais=vector(firstobs,lastobs); + moisdc=vector(firstobs,lastobs); + andc=vector(firstobs,lastobs); + weight=vector(firstobs,lastobs); + agedc=vector(firstobs,lastobs); + cod=ivector(firstobs,lastobs); + for(i=firstobs;i<=lastobs;i++){ num[i]=0; moisnais[i]=0; annais[i]=0; @@ -11142,9 +11389,9 @@ Please run with mle=-1 to get a correct cod[i]=0; weight[i]=1.0; /* Equal weights, 1 by default */ } - mint=matrix(1,maxwav,1,n); - anint=matrix(1,maxwav,1,n); - s=imatrix(1,maxwav+1,1,n); /* s[i][j] health state for wave i and individual j */ + mint=matrix(1,maxwav,firstobs,lastobs); + anint=matrix(1,maxwav,firstobs,lastobs); + s=imatrix(1,maxwav+1,firstobs,lastobs); /* s[i][j] health state for wave i and individual j */ tab=ivector(1,NCOVMAX); ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ ncodemaxwundef=ivector(1,NCOVMAX); /* Number of code per covariate; if - 1 O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ @@ -11246,8 +11493,8 @@ Please run with mle=-1 to get a correct agegomp=(int)agemin; - free_vector(moisnais,1,n); - free_vector(annais,1,n); + free_vector(moisnais,firstobs,lastobs); + free_vector(annais,firstobs,lastobs); /* free_matrix(mint,1,maxwav,1,n); free_matrix(anint,1,maxwav,1,n);*/ /* free_vector(moisdc,1,n); */ @@ -11273,8 +11520,8 @@ Please run with mle=-1 to get a correct concatwav(wav, dh, bh, mw, s, agedc, agev, firstpass, lastpass, imx, nlstate, stepm); /* Concatenates waves */ - free_vector(moisdc,1,n); - free_vector(andc,1,n); + free_vector(moisdc,firstobs,lastobs); + free_vector(andc,firstobs,lastobs); /* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */ nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); @@ -11456,10 +11703,31 @@ Title=%s
      Datafile=%s Firstpass=%d La firstpass, lastpass, stepm, weightopt, model); fprintf(fichtm,"\n"); - fprintf(fichtm,"
      Total number of observations=%d
      \n\ + fprintf(fichtm,"

      Parameter line 2

      • Tolerance for the convergence of the likelihood: ftol=%g \n
      • Interval for the elementary matrix (in month): stepm=%d",\ + ftol, stepm); + fprintf(fichtm,"\n
      • Number of fixed dummy covariates: ncovcol=%d ", ncovcol); + ncurrv=1; + for(i=ncurrv; i <=ncovcol; i++) fprintf(fichtm,"V%d ", i); + fprintf(fichtm,"\n
      • Number of fixed quantitative variables: nqv=%d ", nqv); + ncurrv=i; + for(i=ncurrv; i <=ncurrv-1+nqv; i++) fprintf(fichtm,"V%d ", i); + fprintf(fichtm,"\n
      • Number of time varying (wave varying) dummy covariates: ntv=%d ", ntv); + ncurrv=i; + for(i=ncurrv; i <=ncurrv-1+ntv; i++) fprintf(fichtm,"V%d ", i); + fprintf(fichtm,"\n
      • Number of time varying quantitative covariates: nqtv=%d ", nqtv); + ncurrv=i; + for(i=ncurrv; i <=ncurrv-1+nqtv; i++) fprintf(fichtm,"V%d ", i); + fprintf(fichtm,"\n
      • Weights column \n
        Number of alive states: nlstate=%d
        Number of death states (not really implemented): ndeath=%d \n
      • Number of waves: maxwav=%d \n
      • Parameter for maximization (1), using parameter values (0), for design of parameters and variance-covariance matrix: mle=%d \n
      • Does the weight column be taken into account (1), or not (0): weight=%d
      \n", \ + nlstate, ndeath, maxwav, mle, weightopt); + + fprintf(fichtm,"

      Diagram of states %s_.svg

      \n\ +", subdirf2(optionfilefiname,"D_"),subdirf2(optionfilefiname,"D_"),subdirf2(optionfilefiname,"D_")); + + + fprintf(fichtm,"\n

      Some descriptive statistics

      \n
      Total number of observations=%d
      \n\ Youngest age at first (selected) pass %.2f, oldest age %.2f
      \n\ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
      \n",\ - imx,agemin,agemax,jmin,jmax,jmean); + imx,agemin,agemax,jmin,jmax,jmean); pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ @@ -11478,10 +11746,10 @@ Interval (in months) between two waves: for(j=1;j<=NDIM;j++) ximort[i][j]=0.; /* ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */ - cens=ivector(1,n); - ageexmed=vector(1,n); - agecens=vector(1,n); - dcwave=ivector(1,n); + cens=ivector(firstobs,lastobs); + ageexmed=vector(firstobs,lastobs); + agecens=vector(firstobs,lastobs); + dcwave=ivector(firstobs,lastobs); for (i=1; i<=imx; i++){ dcwave[i]=-1; @@ -11695,9 +11963,10 @@ Please run with mle=-1 to get a correct free_vector(lpop,1,AGESUP); free_vector(tpop,1,AGESUP); free_matrix(ximort,1,NDIM,1,NDIM); - free_ivector(cens,1,n); - free_vector(agecens,1,n); - free_ivector(dcwave,1,n); + free_ivector(dcwave,firstobs,lastobs); + free_vector(agecens,firstobs,lastobs); + free_vector(ageexmed,firstobs,lastobs); + free_ivector(cens,firstobs,lastobs); #ifdef GSL #endif } /* Endof if mle==-3 mortality only */ @@ -11731,7 +12000,7 @@ Please run with mle=-1 to get a correct printf("\n"); /*--------- results files --------------*/ - fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, weightopt,model); + /* fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nqv=%d ntv=%d nqtv=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nqv, ntv, nqtv, nlstate, ndeath, maxwav, weightopt,model); */ fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); @@ -12019,6 +12288,9 @@ Please run with mle=-1 to get a correct fprintf(ficlog,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); /* day and month of proj2 are not used but only year anproj2.*/ + dateproj1=anproj1+(mproj1-1)/12.+(jproj1-1)/365.; + dateproj2=anproj2+(mproj2-1)/12.+(jproj2-1)/365.; + } break; case 12: @@ -12034,6 +12306,8 @@ Please run with mle=-1 to get a correct fprintf(ficlog,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); fprintf(ficres,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); /* day and month of proj2 are not used but only year anproj2.*/ + dateback1=anback1+(mback1-1)/12.+(jback1-1)/365.; + dateback2=anback2+(mback2-1)/12.+(jback2-1)/365.; } break; case 13: @@ -12089,7 +12363,7 @@ Please run with mle=-1 to get a correct } printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt, \ model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,backcast, estepm, \ - jprev1,mprev1,anprev1,dateprev1,jprev2,mprev2,anprev2,dateprev2); + jprev1,mprev1,anprev1,dateprev1, dateproj1, dateback1,jprev2,mprev2,anprev2,dateprev2,dateproj2, dateback2); /*------------ free_vector -------------*/ /* chdir(path); */ @@ -12098,8 +12372,8 @@ Please run with mle=-1 to get a correct /* free_imatrix(dh,1,lastpass-firstpass+2,1,imx); */ /* free_imatrix(bh,1,lastpass-firstpass+2,1,imx); */ /* free_imatrix(mw,1,lastpass-firstpass+2,1,imx); */ - free_lvector(num,1,n); - free_vector(agedc,1,n); + free_lvector(num,firstobs,lastobs); + free_vector(agedc,firstobs,lastobs); /*free_matrix(covar,0,NCOVMAX,1,n);*/ /*free_matrix(covar,1,NCOVMAX,1,n);*/ fclose(ficparo); @@ -12356,13 +12630,13 @@ Please run with mle=-1 to get a correct if(vpopbased==1) fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav); else - fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n"); + fprintf(ficrest,"the age specific forward period (stable) prevalences in each health state \n"); fprintf(ficrest,"# Age popbased mobilav e.. (std) "); for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); fprintf(ficrest,"\n"); /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */ - printf("Computing age specific period (stable) prevalences in each health state \n"); - fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n"); + printf("Computing age specific forward period (stable) prevalences in each health state \n"); + fprintf(ficlog,"Computing age specific forward period (stable) prevalences in each health state \n"); for(age=bage; age <=fage ;age++){ prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k, nres); /*ZZ Is it the correct prevalim */ if (vpopbased==1) { @@ -12409,16 +12683,16 @@ Please run with mle=-1 to get a correct printf("done State-specific expectancies\n");fflush(stdout); fprintf(ficlog,"done State-specific expectancies\n");fflush(ficlog); - /* variance-covariance of period prevalence*/ + /* variance-covariance of forward period prevalence*/ varprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, prlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff); - free_vector(weight,1,n); + free_vector(weight,firstobs,lastobs); free_imatrix(Tvard,1,NCOVMAX,1,2); - free_imatrix(s,1,maxwav+1,1,n); - free_matrix(anint,1,maxwav,1,n); - free_matrix(mint,1,maxwav,1,n); - free_ivector(cod,1,n); + free_imatrix(s,1,maxwav+1,firstobs,lastobs); + free_matrix(anint,1,maxwav,firstobs,lastobs); + free_matrix(mint,1,maxwav,firstobs,lastobs); + free_ivector(cod,firstobs,lastobs); free_ivector(tab,1,NCOVMAX); fclose(ficresstdeij); fclose(ficrescveij); @@ -12438,10 +12712,10 @@ Please run with mle=-1 to get a correct free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath); free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath); - if(ntv+nqtv>=1)free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,1,n); - if(nqtv>=1)free_ma3x(cotqvar,1,maxwav,1,nqtv,1,n); - if(nqv>=1)free_matrix(coqvar,1,nqv,1,n); - free_matrix(covar,0,NCOVMAX,1,n); + if(ntv+nqtv>=1)free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,firstobs,lastobs); + if(nqtv>=1)free_ma3x(cotqvar,1,maxwav,1,nqtv,firstobs,lastobs); + if(nqv>=1)free_matrix(coqvar,1,nqv,firstobs,lastobs); + free_matrix(covar,0,NCOVMAX,firstobs,lastobs); free_matrix(matcov,1,npar,1,npar); free_matrix(hess,1,npar,1,npar); /*free_vector(delti,1,npar);*/ @@ -12524,6 +12798,8 @@ Please run with mle=-1 to get a correct fclose(ficlog); /*------ End -----------*/ + +/* Executes gnuplot */ printf("Before Current directory %s!\n",pathcd); #ifdef WIN32 @@ -12559,13 +12835,16 @@ Please run with mle=-1 to get a correct sprintf(plotcmd,"%s %s",pplotcmd, optionfilegnuplot); printf("Starting graphs with: '%s'\n",plotcmd);fflush(stdout); + strcpy(pplotcmd,plotcmd); if((outcmd=system(plotcmd)) != 0){ - printf("gnuplot command might not be in your path: '%s', err=%d\n", plotcmd, outcmd); + printf("Error in gnuplot, command might not be in your path: '%s', err=%d\n", plotcmd, outcmd); printf("\n Trying if gnuplot resides on the same directory that IMaCh\n"); sprintf(plotcmd,"%sgnuplot %s", pathimach, optionfilegnuplot); - if((outcmd=system(plotcmd)) != 0) + if((outcmd=system(plotcmd)) != 0){ printf("\n Still a problem with gnuplot command %s, err=%d\n", plotcmd, outcmd); + strcpy(plotcmd,pplotcmd); + } } printf(" Successful, please wait..."); while (z[0] != 'q') { @@ -12592,4 +12871,6 @@ end: printf("\nType q for exiting: "); fflush(stdout); scanf("%s",z); } + printf("End\n"); + exit(0); }