--- imach/src/imach.c 2017/04/26 16:22:11 1.265 +++ imach/src/imach.c 2017/05/13 07:26:12 1.266 @@ -1,6 +1,9 @@ -/* $Id: imach.c,v 1.265 2017/04/26 16:22:11 brouard Exp $ +/* $Id: imach.c,v 1.266 2017/05/13 07:26:12 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.266 2017/05/13 07:26:12 brouard + Summary: Version 0.99r13 (improvements and bugs fixed) + Revision 1.265 2017/04/26 16:22:11 brouard Summary: imach 0.99r13 Some bugs fixed @@ -997,12 +1000,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.265 2017/04/26 16:22:11 brouard Exp $ */ +/* $Id: imach.c,v 1.266 2017/05/13 07:26:12 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.265 $ $Date: 2017/04/26 16:22:11 $"; +char fullversion[]="$Revision: 1.266 $ $Date: 2017/05/13 07:26:12 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -2657,12 +2660,12 @@ Earliest age to start was %d-%d=%d, ncvl max=vector(1,nlstate); meandiff=vector(1,nlstate); - dnewm=ddnewms; doldm=ddoldms; dsavm=ddsavms; - oldm=oldms; savm=savms; - - /* Starting with matrix unity */ - for (ii=1;ii<=nlstate+ndeath;ii++) - for (j=1;j<=nlstate+ndeath;j++){ + dnewm=ddnewms; doldm=ddoldms; dsavm=ddsavms; + oldm=oldms; savm=savms; + + /* Starting with matrix unity */ + for (ii=1;ii<=nlstate+ndeath;ii++) + for (j=1;j<=nlstate+ndeath;j++){ oldm[ii][j]=(ii==j ? 1.0 : 0.0); } @@ -2736,8 +2739,27 @@ Earliest age to start was %d-%d=%d, ncvl /* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, ageminpar, agemaxpar, dnewm, doldm, dsavm,ij)); /\* Bug Valgrind *\/ */ /* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij)); /\* Bug Valgrind *\/ */ out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent,ij)); /* Bug Valgrind */ + /* if((int)age == 70){ */ + /* printf(" Backward prevalim age=%d agefin=%d \n", (int) age, (int) agefin); */ + /* for(i=1; i<=nlstate+ndeath; i++) { */ + /* printf("%d newm= ",i); */ + /* for(j=1;j<=nlstate+ndeath;j++) { */ + /* printf("%f ",newm[i][j]); */ + /* } */ + /* printf("oldm * "); */ + /* for(j=1;j<=nlstate+ndeath;j++) { */ + /* printf("%f ",oldm[i][j]); */ + /* } */ + /* printf(" pmmij "); */ + /* for(j=1;j<=nlstate+ndeath;j++) { */ + /* printf("%f ",pmmij[i][j]); */ + /* } */ + /* printf("\n"); */ + /* } */ + /* } */ savm=oldm; oldm=newm; + for(j=1; j<=nlstate; j++){ max[j]=0.; min[j]=1.; @@ -2788,7 +2810,7 @@ Oldest age to start was %d-%d=%d, ncvloo double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate ) { /* According to parameters values stored in x and the covariate's values stored in cov, - computes the probability to be observed in state j being in state i by appying the + computes the probability to be observed in state j (after stepm years) being in state i by appying the model to the ncovmodel covariates (including constant and age). lnpijopii=ln(pij/pii)= aij+bij*age+cij*v1+dij*v2+... = sum_nc=1^ncovmodel xij(nc)*cov[nc] and, according on how parameters are entered, the position of the coefficient xij(nc) of the @@ -2797,8 +2819,9 @@ double **pmij(double **ps, double *cov, j>=i nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel Computes ln(pij/pii) (lnpijopii), deduces pij/pii by exponentiation, sums on j different of i to get 1-pii/pii, deduces pii, and then all pij. - Outputs ps[i][j] the probability to be observed in j being in j according to + Outputs ps[i][j] or probability to be observed in j being in i according to the values of the covariates cov[nc] and corresponding parameter values x[nc+shiftij] + Sum on j ps[i][j] should equal to 1. */ double s1, lnpijopii; /*double t34;*/ @@ -2862,7 +2885,7 @@ double **pmij(double **ps, double *cov, /* for(i=1; i<= npar; i++) printf("%f ",x[i]); goto end;*/ - return ps; + return ps; /* Pointer is unchanged since its call */ } /*************** backward transition probabilities ***************/ @@ -2871,8 +2894,8 @@ double **pmij(double **ps, double *cov, /* double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, double ***dnewm, double **doldm, double **dsavm, int ij ) */ double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, int ij ) { - /* Computes the backward probability at age agefin and covariate ij - * and returns in **ps as well as **bmij. + /* Computes the backward probability at age agefin and covariate combination ij. In fact cov is already filled and x too. + * Call to pmij(cov and x), call to cross prevalence, sums and inverses, left multiply, and returns in **ps as well as **bmij. */ int i, ii, j,k; @@ -2889,43 +2912,56 @@ double **pmij(double **ps, double *cov, agefin=cov[2]; /* bmij *//* age is cov[2], ij is included in cov, but we need for - the observed prevalence (with this covariate ij) */ - dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate); + the observed prevalence (with this covariate ij) at beginning of transition */ + /* dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate); */ + pmmij=pmij(pmmij,cov,ncovmodel,x,nlstate); /*This is forward probability from agefin to agefin + stepm */ + /* outputs pmmij which is a stochastic matrix */ /* We do have the matrix Px in savm and we need pij */ for (j=1;j<=nlstate+ndeath;j++){ - sumnew=0.; /* w1 p11 + w2 p21 only on live states */ + sumnew=0.; /* w1 p11 + w2 p21 only on live states N1./N..*N11/N1. + N2./N..*N21/N2.=(N11+N21)/N..=N.1/N.. */ for (ii=1;ii<=nlstate;ii++){ - sumnew+=dsavm[ii][j]*prevacurrent[(int)agefin][ii][ij]; + /* sumnew+=dsavm[ii][j]*prevacurrent[(int)agefin][ii][ij]; */ + sumnew+=pmmij[ii][j]*prevacurrent[(int)agefin][ii][ij]; /* Yes prevalence at beginning of transition */ } /* sumnew is (N11+N21)/N..= N.1/N.. = sum on i of w_i pij */ - for (ii=1;ii<=nlstate+ndeath;ii++){ - if(sumnew >= 1.e-10){ + if(sumnew >= 1.e-10){ + for (ii=1;ii<=nlstate+ndeath;ii++){ /* if(agefin >= agemaxpar && agefin <= agemaxpar+stepm/YEARM){ */ /* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */ /* }else if(agefin >= agemaxpar+stepm/YEARM){ */ /* doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); */ /* }else */ doldm[ii][j]=(ii==j ? 1./sumnew : 0.0); - }else{ - ; - /* printf("ii=%d, i=%d, doldm=%lf dsavm=%lf, probs=%lf, sumnew=%lf,agefin=%d\n",ii,j,doldm[ii][j],dsavm[ii][j],prevacurrent[(int)agefin][ii][ij],sumnew, (int)agefin); */ - } - } /*End ii */ - } /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] */ - /* left Product of this diag matrix by dsavm=Px (newm=dsavm*doldm) */ - bbmij=matprod2(dnewm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, doldm); /* Bug Valgrind */ + } /*End ii */ + }else{ /* We put the identity matrix */ + for (ii=1;ii<=nlstate+ndeath;ii++){ + doldm[ii][j]=(ii==j ? 1. : 0.0); + } /*End ii */ + /* printf("Problem internal bmij A: sum_i w_i*p_ij=N.j/N.. <1.e-10 i=%d, j=%d, sumnew=%lf,agefin=%d\n",ii,j,sumnew, (int)agefin); */ + } + } /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] or identity*/ + /* left Product of this diag matrix by dsavm=Px (dnewm=dsavm*doldm) */ + /* bbmij=matprod2(dnewm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, doldm); /\* Bug Valgrind *\/ */ + bbmij=matprod2(dnewm, pmmij,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, doldm); /* Bug Valgrind */ /* dsavm=doldm; /\* dsavm is now diag [1/(w_1p1i+w_2 p2i)] but can be overwritten*\/ */ /* doldm=dnewm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */ /* dnewm=dsavm; /\* doldm is now Px * diag [1/(w_1p1i+w_2 p2i)] *\/ */ /* left Product of this matrix by diag matrix of prevalences (savm) */ for (j=1;j<=nlstate+ndeath;j++){ + sumnew=0.; for (ii=1;ii<=nlstate+ndeath;ii++){ + sumnew+=prevacurrent[(int)agefin][ii][ij]; dsavm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij] : 0.0); } - } /* End j, At the end oldm is diag[1/(w_1p1i+w_2 p2i)] */ - ps=matprod2(doldm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dnewm); /* Bug Valgrind */ + /* if(sumnew <0.9){ */ + /* printf("Problem internal bmij B: sum on i wi <0.9: j=%d, sum_i wi=%lf,agefin=%d\n",j,sumnew, (int)agefin); */ + /* } */ + } /* End j, At the end dsavm is diag[(w_i)] */ + /* What if dsavm doesn't sum ii to 1? */ + /* ps=matprod2(doldm, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dnewm); /\* Bug Valgrind *\/ */ + ps=matprod2(ps, dsavm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dnewm); /* Bug Valgrind */ /* newm or out is now diag[w_i] * Px * diag [1/(w_1p1i+w_2 p2i)] */ /* end bmij */ - return ps; + return ps; /*pointer is unchanged */ } /*************** transition probabilities ***************/ @@ -3151,7 +3187,7 @@ double ***hpxij(double ***po, int nhstep /* double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, int ij ) */ double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, int ij ) { - /* Computes the transition matrix starting at age 'age' over + /* For a combination of dummy covariate ij, computes the transition matrix starting at age 'age' over 'nhstepm*hstepm*stepm' months (i.e. until age (in years) age+nhstepm*hstepm*stepm/12) by multiplying nhstepm*hstepm matrices. @@ -3159,18 +3195,19 @@ double ***hbxij(double ***po, int nhstep (typically every 2 years instead of every month which is too big for the memory). Model is determined by parameters x and covariates have to be - included manually here. - + included manually here. Then we use a call to bmij(x and cov) + The addresss of po (p3mat allocated to the dimension of nhstepm) should be stored for output */ int i, j, d, h, k; - double **out, cov[NCOVMAX+1]; - double **newm; + double **out, cov[NCOVMAX+1], **bmij(); + double **newm, ***newmm; double agexact; double agebegin, ageend; double **oldm, **savm; - oldm=oldms;savm=savms; + newmm=po; /* To be saved */ + oldm=oldms;savm=savms; /* Global pointers */ /* Hstepm could be zero and should return the unit matrix */ for (i=1;i<=nlstate+ndeath;i++) for (j=1;j<=nlstate+ndeath;j++){ @@ -3183,14 +3220,18 @@ double ***hbxij(double ***po, int nhstep newm=savm; /* Covariates have to be included here again */ cov[1]=1.; - agexact=age-((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */ + 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) cov[3]= agexact*agexact; - for (k=1; k<=cptcovn;k++) - cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; - /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */ + for (k=1; k<=cptcovn;k++){ + /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; */ + /* /\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; *\/ */ + cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)]; + /* printf("hbxij Dummy agexact=%.0f combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov[%d]=%lf codtabm(%d,Tvar[%d])=%d \n",agexact,ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],2+nagesqr+TvarsDind[k],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */ + + } for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */ /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; @@ -3199,11 +3240,10 @@ double ***hbxij(double ***po, int nhstep cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */ - /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/ /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/ /* Careful transposed matrix */ - /* age is in cov[2] */ + /* age is in cov[2], prevacurrent at beginning of transition. */ /* out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ */ /* 1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); */ out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent,ij),\ @@ -4947,7 +4987,10 @@ void prevalence(double ***probs, double } else{ if(first==1){ first=0; - printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others in log file...\n",jk,i,j1,probs[i][jk][j1]); + 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]); } } } @@ -6185,7 +6228,7 @@ void varprob(char optionfilefiname[], do fprintf(fichtm,"\n
Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year-1\
: \
%s_%d%1d%1d-%1d%1d.svg, ",k1,l1,k2,l2,\
subdirf2(optionfilefiname,"VARPIJGR_"), j1,k1,l1,k2,l2, \
@@ -6413,16 +6456,16 @@ 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(lc2), \
- mu2,std,v21,sqrt(lc1),v22,sqrt(lc2));
+ mu1,std,v11,sqrt(lc1),v12,sqrt(fabs(lc2)), \
+ mu2,std,v21,sqrt(lc1),v22,sqrt(fabs(lc2))); /* For gnuplot only */
}else{
first=0;
fprintf(fichtmcov," %d (%.3f),",(int) age, c12);
fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);
fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);
fprintf(ficgp,"\nreplot %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(lc2), \
- mu2,std,v21,sqrt(lc1),v22,sqrt(lc2));
+ mu1,std,v11,sqrt(lc1),v12,sqrt(fabs(lc2)), \
+ mu2,std,v21,sqrt(lc1),v22,sqrt(fabs(lc2)));
}/* if first */
} /* age mod 5 */
} /* end loop age */
@@ -6713,7 +6756,7 @@ true period expectancies (those weighted
}
/******************* Gnuplot file **************/
-void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[]){
+ void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[], int offyear){
char dirfileres[132],optfileres[132];
char gplotcondition[132], gplotlabel[132];
@@ -6723,6 +6766,7 @@ void printinggnuplot(char fileresu[], ch
int vpopbased;
int ioffset; /* variable offset for columns */
int nres=0; /* Index of resultline */
+ int istart=1; /* For starting graphs in projections */
/* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */
/* printf("Problem with file %s",optionfilegnuplot); */
@@ -7253,12 +7297,16 @@ set ter svg size 640, 480\nunset log y\n
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 xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar);
- for (i=1; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */
+
+ /* for (i=1; i<= nlstate+1 ; i ++){ /\* nlstate +1 p11 p21 p.1 *\/ */
+ istart=nlstate+1; /* Could be one if by state, but nlstate+1 is w.i projection only */
+ /*istart=1;*/ /* Could be one if by state, but nlstate+1 is w.i projection only */
+ for (i=istart; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */
/*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
/*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
/*# yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
/*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
- if(i==1){
+ if(i==istart){
fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"F_"));
}else{
fprintf(ficgp,",\\\n '' ");
@@ -7270,10 +7318,15 @@ set ter svg size 640, 480\nunset log y\n
/*# V1 = 1 yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/
/*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */
fprintf(ficgp," u %d:(", ioffset);
- if(i==nlstate+1)
- fprintf(ficgp," $%d/(1.-$%d)) t 'pw.%d' with line ", \
+ if(i==nlstate+1){
+ fprintf(ficgp," $%d/(1.-$%d)):5 t 'pw.%d' with line lc variable ", \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
+ fprintf(ficgp,",\\\n '' ");
+ fprintf(ficgp," u %d:(",ioffset);
+ fprintf(ficgp," (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", \
+ offyear, \
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
- else
+ }else
fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,i,cpt );
}else{ /* more than 2 covariates */
@@ -7305,8 +7358,14 @@ set ter svg size 640, 480\nunset log y\n
/*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */
/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
if(i==nlstate+1){
- fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ", gplotcondition, \
+ fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0):5 t 'p.%d' with line lc variable", gplotcondition, \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
+ fprintf(ficgp,",\\\n '' ");
+ fprintf(ficgp," u %d:(",ioffset);
+ fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", gplotcondition, \
+ offyear, \
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
+/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0) && (($5-$6) == 1947) ? $10/(1.-$22) : 1/0):5 with labels center boxed not*/
}else{
fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset +1+(i-1)+(nlstate+1)*nlstate,i,cpt );
@@ -7526,10 +7585,11 @@ set ter svg size 640, 480\nunset log y\n
int mobilavrange, mob;
int iage=0;
- double sum=0.;
+ double sum=0., sumr=0.;
double age;
- double *sumnewp, *sumnewm;
- double *agemingood, *agemaxgood; /* Currently identical for all covariates */
+ double *sumnewp, *sumnewm, *sumnewmr;
+ double *agemingood, *agemaxgood;
+ double *agemingoodr, *agemaxgoodr;
/* modcovmax=2*cptcoveff;/\* Max number of modalities. We suppose */
@@ -7537,19 +7597,22 @@ set ter svg size 640, 480\nunset log y\n
sumnewp = vector(1,ncovcombmax);
sumnewm = vector(1,ncovcombmax);
+ sumnewmr = vector(1,ncovcombmax);
agemingood = vector(1,ncovcombmax);
+ agemingoodr = vector(1,ncovcombmax);
agemaxgood = vector(1,ncovcombmax);
+ agemaxgoodr = vector(1,ncovcombmax);
for (cptcod=1;cptcod<=ncovcombmax;cptcod++){
- sumnewm[cptcod]=0.;
+ sumnewm[cptcod]=0.; sumnewmr[cptcod]=0.;
sumnewp[cptcod]=0.;
- agemingood[cptcod]=0;
- agemaxgood[cptcod]=0;
+ agemingood[cptcod]=0, agemingoodr[cptcod]=0;
+ agemaxgood[cptcod]=0, agemaxgoodr[cptcod]=0;
}
if (cptcovn<1) ncovcombmax=1; /* At least 1 pass */
- if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){
- if(mobilav==1) mobilavrange=5; /* default */
+ if(mobilav==-1 || mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){
+ if(mobilav==1 || mobilav==-1) mobilavrange=5; /* default */
else mobilavrange=mobilav;
for (age=bage; age<=fage; age++)
for (i=1; i<=nlstate;i++)
@@ -7561,77 +7624,143 @@ set ter svg size 640, 480\nunset log y\n
*/
for (mob=3;mob <=mobilavrange;mob=mob+2){
for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){
- for (i=1; i<=nlstate;i++){
- for (cptcod=1;cptcod<=ncovcombmax;cptcod++){
+ for (cptcod=1;cptcod<=ncovcombmax;cptcod++){
+ sumnewm[cptcod]=0.;
+ for (i=1; i<=nlstate;i++){
mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod];
for (cpt=1;cpt<=(mob-1)/2;cpt++){
mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod];
mobaverage[(int)age][i][cptcod] +=probs[(int)age+cpt][i][cptcod];
}
mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/mob;
- }
- }
+ sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
+ } /* end i */
+ if(sumnewm[cptcod] >1.e-3) mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/sumnewm[cptcod]; /* Rescaling to sum one */
+ } /* end cptcod */
}/* end age */
}/* end mob */
- }else
+ }else{
+ printf("Error internal in movingaverage, mobilav=%d.\n",mobilav);
return -1;
- for (cptcod=1;cptcod<=ncovcombmax;cptcod++){
+ }
+
+ for (cptcod=1;cptcod<=ncovcombmax;cptcod++){ /* for each combination */
/* for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ */
if(invalidvarcomb[cptcod]){
printf("\nCombination (%d) ignored because no cases \n",cptcod);
continue;
}
- agemingood[cptcod]=fage-(mob-1)/2;
- for (age=fage-(mob-1)/2; age>=bage; age--){/* From oldest to youngest, finding the youngest wrong */
+ for (age=fage-(mob-1)/2; age>=bage+(mob-1)/2; age--){ /*looking for the youngest and oldest good age */
+ sumnewm[cptcod]=0.;
+ sumnewmr[cptcod]=0.;
+ for (i=1; i<=nlstate;i++){
+ sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
+ sumnewmr[cptcod]+=probs[(int)age][i][cptcod];
+ }
+ if(fabs(sumnewmr[cptcod] - 1.) <= 1.e-3) { /* good without smoothing */
+ agemingoodr[cptcod]=age;
+ }
+ if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */
+ agemingood[cptcod]=age;
+ }
+ } /* age */
+ for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ /*looking for the youngest and oldest good age */
sumnewm[cptcod]=0.;
+ sumnewmr[cptcod]=0.;
for (i=1; i<=nlstate;i++){
sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
+ sumnewmr[cptcod]+=probs[(int)age][i][cptcod];
+ }
+ if(fabs(sumnewmr[cptcod] - 1.) <= 1.e-3) { /* good without smoothing */
+ agemaxgoodr[cptcod]=age;
}
if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */
- agemingood[cptcod]=age;
- }else{ /* bad */
- for (i=1; i<=nlstate;i++){
- mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod];
- } /* i */
+ agemaxgood[cptcod]=age;
+ }
+ } /* age */
+ /* Thus we have agemingood and agemaxgood as well as goodr for raw (preobs) */
+ /* but they will change */
+ for (age=fage-(mob-1)/2; age>=bage; age--){/* From oldest to youngest, filling up to the youngest */
+ sumnewm[cptcod]=0.;
+ sumnewmr[cptcod]=0.;
+ for (i=1; i<=nlstate;i++){
+ sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
+ sumnewmr[cptcod]+=probs[(int)age][i][cptcod];
+ }
+ if(mobilav==-1){ /* Forcing raw ages if good else agemingood */
+ if(fabs(sumnewmr[cptcod] - 1.) <= 1.e-3) { /* good without smoothing */
+ agemaxgoodr[cptcod]=age; /* age min */
+ for (i=1; i<=nlstate;i++)
+ mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod];
+ }else{ /* bad we change the value with the values of good ages */
+ for (i=1; i<=nlstate;i++){
+ mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgoodr[cptcod]][i][cptcod];
+ } /* i */
+ } /* end bad */
+ }else{
+ if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */
+ agemaxgood[cptcod]=age;
+ }else{ /* bad we change the value with the values of good ages */
+ for (i=1; i<=nlstate;i++){
+ mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
+ } /* i */
+ } /* end bad */
+ }/* end else */
+ sum=0.;sumr=0.;
+ for (i=1; i<=nlstate;i++){
+ sum+=mobaverage[(int)age][i][cptcod];
+ 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, 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, bage);
} /* end bad */
}/* age */
- sum=0.;
- for (i=1; i<=nlstate;i++){
- sum+=mobaverage[(int)agemingood[cptcod]][i][cptcod];
- }
- if(fabs(sum - 1.) > 1.e-3) { /* bad */
- printf("For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one at any descending age!\n",cptcod);
- /* for (i=1; i<=nlstate;i++){ */
- /* mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */
- /* } /\* i *\/ */
- } /* end bad */
- /* else{ /\* We found some ages summing to one, we will smooth the oldest *\/ */
- /* From youngest, finding the oldest wrong */
- agemaxgood[cptcod]=bage+(mob-1)/2;
- for (age=bage+(mob-1)/2; age<=fage; age++){
+
+ for (age=bage+(mob-1)/2; age<=fage; age++){/* From youngest, finding the oldest wrong */
sumnewm[cptcod]=0.;
+ sumnewmr[cptcod]=0.;
for (i=1; i<=nlstate;i++){
sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
+ sumnewmr[cptcod]+=probs[(int)age][i][cptcod];
+ }
+ if(mobilav==-1){ /* Forcing raw ages if good else agemingood */
+ if(fabs(sumnewmr[cptcod] - 1.) <= 1.e-3) { /* good */
+ agemingoodr[cptcod]=age;
+ for (i=1; i<=nlstate;i++)
+ mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod];
+ }else{ /* bad we change the value with the values of good ages */
+ for (i=1; i<=nlstate;i++){
+ mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingoodr[cptcod]][i][cptcod];
+ } /* i */
+ } /* end bad */
+ }else{
+ if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */
+ agemingood[cptcod]=age;
+ }else{ /* bad */
+ for (i=1; i<=nlstate;i++){
+ mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod];
+ } /* i */
+ } /* end bad */
+ }/* end else */
+ sum=0.;sumr=0.;
+ for (i=1; i<=nlstate;i++){
+ sum+=mobaverage[(int)age][i][cptcod];
+ sumr+=mobaverage[(int)age][i][cptcod];
}
- if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */
- agemaxgood[cptcod]=age;
- }else{ /* bad */
- for (i=1; i<=nlstate;i++){
- mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
- } /* i */
+ if(fabs(sum - 1.) > 1.e-3) { /* bad */
+ printf("Moving average B1: 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 decrease fage=%d?\n",cptcod, sum, (int) age, fage);
+ } /* 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 B2: 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 fage=%d\n",cptcod,sumr, (int)age, fage);
} /* end bad */
}/* age */
- sum=0.;
- for (i=1; i<=nlstate;i++){
- sum+=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
- }
- if(fabs(sum - 1.) > 1.e-3) { /* bad */
- printf("For this combination of covariate cptcod=%d, we can't get a smoothed prevalence which sums to one at any ascending age!\n",cptcod);
- /* for (i=1; i<=nlstate;i++){ */
- /* mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */
- /* } /\* i *\/ */
- } /* end bad */
+
for (age=bage; age<=fage; age++){
/* printf("%d %d ", cptcod, (int)age); */
@@ -7646,28 +7775,32 @@ set ter svg size 640, 480\nunset log y\n
}
/* printf("\n"); */
/* } */
+
/* brutal averaging */
- for (i=1; i<=nlstate;i++){
- for (age=1; age<=bage; age++){
- mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod];
- /* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); */
- }
- for (age=fage; age<=AGESUP; age++){
- mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
- /* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); */
- }
- } /* end i status */
- for (i=nlstate+1; i<=nlstate+ndeath;i++){
- for (age=1; age<=AGESUP; age++){
- /*printf("i=%d, age=%d, cptcod=%d\n",i, (int)age, cptcod);*/
- mobaverage[(int)age][i][cptcod]=0.;
- }
- }
+ /* for (i=1; i<=nlstate;i++){ */
+ /* for (age=1; age<=bage; age++){ */
+ /* mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; */
+ /* /\* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); *\/ */
+ /* } */
+ /* for (age=fage; age<=AGESUP; age++){ */
+ /* mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod]; */
+ /* /\* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); *\/ */
+ /* } */
+ /* } /\* end i status *\/ */
+ /* for (i=nlstate+1; i<=nlstate+ndeath;i++){ */
+ /* for (age=1; age<=AGESUP; age++){ */
+ /* /\*printf("i=%d, age=%d, cptcod=%d\n",i, (int)age, cptcod);*\/ */
+ /* mobaverage[(int)age][i][cptcod]=0.; */
+ /* } */
+ /* } */
}/* end cptcod */
- free_vector(sumnewm,1, ncovcombmax);
- free_vector(sumnewp,1, ncovcombmax);
+ free_vector(agemaxgoodr,1, ncovcombmax);
free_vector(agemaxgood,1, ncovcombmax);
free_vector(agemingood,1, ncovcombmax);
+ free_vector(agemingoodr,1, ncovcombmax);
+ free_vector(sumnewmr,1, ncovcombmax);
+ free_vector(sumnewm,1, ncovcombmax);
+ free_vector(sumnewp,1, ncovcombmax);
return 0;
}/* End movingaverage */
@@ -7774,11 +7907,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)
+ /* if (mobilav>=1) */
ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][k];
- else {
- ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k];
- }
+ /* else { */ /* even if mobilav==-1 we use mobaverage */
+ /* ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */
+ /* } */
if (h*hstepm/YEARM*stepm== yearp) {
fprintf(ficresf," %.3f", p3mat[i][j][h]);
}
@@ -7790,6 +7923,8 @@ set ter svg size 640, 480\nunset log y\n
} /* end h */
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
} /* end agec */
+ /* diffyear=(int) anproj1+yearp-ageminpar-1; */
+ /*printf("Prevforecast %d+%d-%d=diffyear=%d\n",(int) anproj1, (int)yearp,(int)ageminpar,(int) anproj1-(int)ageminpar);*/
} /* end yearp */
} /* end k */
@@ -9771,6 +9906,8 @@ int back_prevalence_limit(double *p, dou
}else{
/* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k,nres);
+ /* printf("TOTOT\n"); */
+ /* exit(1); */
}
fprintf(ficresplb,"%.0f ",age );
for(j=1;j<=cptcoveff;j++)
@@ -9929,7 +10066,9 @@ int hPijx(double *p, int bage, int fage)
/* nhstepm=nhstepm*YEARM; aff par mois*/
- p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); /* We can't have it at an upper level because of nhstepm */
+ /* and memory limitations if stepm is small */
+
/* oldm=oldms;savm=savms; */
/* hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); */
hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k);
@@ -11492,7 +11631,7 @@ Please run with mle=-1 to get a correct
This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\
Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
}else{
- printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, backcast, pathc,p);
+ printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, backcast, pathc,p, (int)anproj1-(int)ageminpar);
}
printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt, \
model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,backcast, estepm, \
@@ -11554,8 +11693,14 @@ Please run with mle=-1 to get a correct
printf(" Error in movingaverage mobilav=%d\n",mobilav);
}
}
- /* /\* Prevalence for each covariates in probs[age][status][cov] *\/ */
- /* prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */
+ /* else if(mobilavproj==-1){ /\* Forcing raw observed prevalences *\/ */
+ /* for(i=1;i<=AGESUP;i++) */
+ /* for(j=1;j<=nlstate;j++) */
+ /* for(k=1;k<=ncovcombmax;k++) */
+ /* mobaverages[i][j][k]=probs[i][j][k]; */
+ /* /\* /\\* Prevalence for each covariates in probs[age][status][cov] *\\/ *\/ */
+ /* /\* prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); *\/ */
+ /* } */
else if (mobilavproj !=0) {
printf("Movingaveraging projected observed prevalence\n");
fprintf(ficlog,"Movingaveraging projected observed prevalence\n");