");
fprintf(ficresphtmfr,"Age | ");
- for(jk=-1; jk <=nlstate+ndeath; jk++){
+ for(s2=-1; s2 <=nlstate+ndeath; s2++){
for(m=-1; m <=nlstate+ndeath; m++){
- if(jk!=0 && m!=0)
- fprintf(ficresphtmfr,"%d%d | ",jk,m);
+ if(s2!=0 && m!=0)
+ fprintf(ficresphtmfr,"%d%d | ",s2,m);
}
}
fprintf(ficresphtmfr, "\n");
@@ -4577,97 +4661,112 @@ Title=%s
Datafile=%s Firstpass=%d La
fprintf(ficresphtmfr,"%d | ",iage);
fprintf(ficlog,"Age %d", iage);
}
- for(jk=1; jk <=nlstate ; jk++){
- for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
- pp[jk] += freq[jk][m][iage];
+ for(s1=1; s1 <=nlstate ; s1++){
+ for(m=-1, pp[s1]=0; m <=nlstate+ndeath ; m++)
+ pp[s1] += freq[s1][m][iage];
}
- for(jk=1; jk <=nlstate ; jk++){
+ for(s1=1; s1 <=nlstate ; s1++){
for(m=-1, pos=0; m <=0 ; m++)
- pos += freq[jk][m][iage];
- if(pp[jk]>=1.e-10){
+ pos += freq[s1][m][iage];
+ if(pp[s1]>=1.e-10){
if(first==1){
- printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
+ printf(" %d.=%.0f loss[%d]=%.1f%%",s1,pp[s1],s1,100*pos/pp[s1]);
}
- fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
+ fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",s1,pp[s1],s1,100*pos/pp[s1]);
}else{
if(first==1)
- printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
- fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
+ printf(" %d.=%.0f loss[%d]=NaNQ%%",s1,pp[s1],s1);
+ fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",s1,pp[s1],s1);
}
}
- for(jk=1; jk <=nlstate ; jk++){
- /* posprop[jk]=0; */
- for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */
- pp[jk] += freq[jk][m][iage];
- } /* pp[jk] is the total number of transitions starting from state jk and any ending status until this age */
-
- for(jk=1,pos=0, pospropta=0.; jk <=nlstate ; jk++){
- pos += pp[jk]; /* pos is the total number of transitions until this age */
- posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state
- from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
- pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state
- from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
+ for(s1=1; s1 <=nlstate ; s1++){
+ /* posprop[s1]=0; */
+ for(m=0, pp[s1]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */
+ pp[s1] += freq[s1][m][iage];
+ } /* pp[s1] is the total number of transitions starting from state s1 and any ending status until this age */
+
+ for(s1=1,pos=0, pospropta=0.; s1 <=nlstate ; s1++){
+ pos += pp[s1]; /* pos is the total number of transitions until this age */
+ posprop[s1] += prop[s1][iage]; /* prop is the number of transitions from a live state
+ from s1 at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
+ pospropta += prop[s1][iage]; /* prop is the number of transitions from a live state
+ from s1 at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
+ }
+
+ /* Writing ficresp */
+ if(cptcoveff==0 && nj==1){ /* no covariate and first pass */
+ if( iage <= iagemax){
+ fprintf(ficresp," %d",iage);
+ }
+ }else if( nj==2){
+ if( iage <= iagemax){
+ fprintf(ficresp," %d",iage);
+ for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, " %d %d",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
+ }
}
- for(jk=1; jk <=nlstate ; jk++){
+ for(s1=1; s1 <=nlstate ; s1++){
if(pos>=1.e-5){
if(first==1)
- printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
- fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
+ printf(" %d.=%.0f prev[%d]=%.1f%%",s1,pp[s1],s1,100*pp[s1]/pos);
+ fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",s1,pp[s1],s1,100*pp[s1]/pos);
}else{
if(first==1)
- printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
- fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
+ printf(" %d.=%.0f prev[%d]=NaNQ%%",s1,pp[s1],s1);
+ fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",s1,pp[s1],s1);
}
if( iage <= iagemax){
if(pos>=1.e-5){
- fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
- /* fprintf(ficresp, "%d %d %d %.5f %.0f %.0f",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)],iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); */
-
- fprintf(ficresphtm,"%d | %.5f | %.0f | %.0f | ",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
- /*probs[iage][jk][j1]= pp[jk]/pos;*/
- /*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/
- }
- else{
- fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta);
- fprintf(ficresphtm,"%d | NaNq | %.0f | %.0f | ",iage, prop[jk][iage],pospropta);
+ if(cptcoveff==0 && nj==1){ /* no covariate and first pass */
+ fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta);
+ }else if( nj==2){
+ fprintf(ficresp," %.5f %.0f %.0f",prop[s1][iage]/pospropta, prop[s1][iage],pospropta);
+ }
+ fprintf(ficresphtm,"%d | %.5f | %.0f | %.0f | ",iage,prop[s1][iage]/pospropta, prop[s1][iage],pospropta);
+ /*probs[iage][s1][j1]= pp[s1]/pos;*/
+ /*printf("\niage=%d s1=%d j1=%d %.5f %.0f %.0f %f",iage,s1,j1,pp[s1]/pos, pp[s1],pos,probs[iage][s1][j1]);*/
+ } else{
+ if((cptcoveff==0 && nj==1)|| nj==2 ) fprintf(ficresp," NaNq %.0f %.0f",prop[s1][iage],pospropta);
+ fprintf(ficresphtm,"%d | NaNq | %.0f | %.0f | ",iage, prop[s1][iage],pospropta);
}
}
- pospropt[jk] +=posprop[jk];
- } /* end loop jk */
+ pospropt[s1] +=posprop[s1];
+ } /* end loop s1 */
/* pospropt=0.; */
- for(jk=-1; jk <=nlstate+ndeath; jk++){
+ for(s1=-1; s1 <=nlstate+ndeath; s1++){
for(m=-1; m <=nlstate+ndeath; m++){
- if(freq[jk][m][iage] !=0 ) { /* minimizing output */
+ if(freq[s1][m][iage] !=0 ) { /* minimizing output */
if(first==1){
- printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]);
+ printf(" %d%d=%.0f",s1,m,freq[s1][m][iage]);
}
- /* printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]); */
- fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]);
+ /* printf(" %d%d=%.0f",s1,m,freq[s1][m][iage]); */
+ fprintf(ficlog," %d%d=%.0f",s1,m,freq[s1][m][iage]);
}
- if(jk!=0 && m!=0)
- fprintf(ficresphtmfr,"%.0f | ",freq[jk][m][iage]);
+ if(s1!=0 && m!=0)
+ fprintf(ficresphtmfr,"%.0f | ",freq[s1][m][iage]);
}
- } /* end loop jk */
+ } /* end loop s1 */
posproptt=0.;
- for(jk=1; jk <=nlstate; jk++){
- posproptt += pospropt[jk];
+ for(s1=1; s1 <=nlstate; s1++){
+ posproptt += pospropt[s1];
}
fprintf(ficresphtmfr,"
\n ");
- if(iage <= iagemax){
- fprintf(ficresp,"\n");
- fprintf(ficresphtm,"\n");
+ fprintf(ficresphtm,"\n");
+ if((cptcoveff==0 && nj==1)|| nj==2 ) {
+ if(iage <= iagemax)
+ fprintf(ficresp,"\n");
}
if(first==1)
printf("Others in log...\n");
fprintf(ficlog,"\n");
} /* end loop age iage */
+
fprintf(ficresphtm,"Tot | ");
- for(jk=1; jk <=nlstate ; jk++){
+ for(s1=1; s1 <=nlstate ; s1++){
if(posproptt < 1.e-5){
- fprintf(ficresphtm,"Nanq | %.0f | %.0f | ",pospropt[jk],posproptt);
+ fprintf(ficresphtm,"Nanq | %.0f | %.0f | ",pospropt[s1],posproptt);
}else{
- fprintf(ficresphtm,"%.5f | %.0f | %.0f | ",pospropt[jk]/posproptt,pospropt[jk],posproptt);
+ fprintf(ficresphtm,"%.5f | %.0f | %.0f | ",pospropt[s1]/posproptt,pospropt[s1],posproptt);
}
}
fprintf(ficresphtm,"
\n");
@@ -4687,66 +4786,66 @@ Title=%s
Datafile=%s Firstpass=%d La
fprintf(ficlog,"\n");
if(j!=0){
printf("#Freqsummary: Starting values for combination j1=%d:\n", j1);
- for(i=1,jk=1; i <=nlstate; i++){
+ for(i=1,s1=1; i <=nlstate; i++){
for(k=1; k <=(nlstate+ndeath); k++){
if (k != i) {
- for(jj=1; jj <=ncovmodel; jj++){ /* For counting jk */
+ for(jj=1; jj <=ncovmodel; jj++){ /* For counting s1 */
if(jj==1){ /* Constant case (in fact cste + age) */
if(j1==1){ /* All dummy covariates to zero */
freq[i][k][iagemax+4]=freq[i][k][iagemax+3]; /* Stores case 0 0 0 */
freq[i][i][iagemax+4]=freq[i][i][iagemax+3]; /* Stores case 0 0 0 */
printf("%d%d ",i,k);
fprintf(ficlog,"%d%d ",i,k);
- printf("%12.7f ln(%.0f/%.0f)= %f, OR=%f sd=%f \n",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]),freq[i][k][iagemax+3]/freq[i][i][iagemax+3], sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3]));
- fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f \n",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
- pstart[jk]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]);
+ printf("%12.7f ln(%.0f/%.0f)= %f, OR=%f sd=%f \n",p[s1],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]),freq[i][k][iagemax+3]/freq[i][i][iagemax+3], sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3]));
+ fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f \n",p[s1],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
+ pstart[s1]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]);
}
}else if((j1==1) && (jj==2 || nagesqr==1)){ /* age or age*age parameter without covariate V4*age (to be done later) */
for(iage=iagemin; iage <= iagemax+3; iage++){
x[iage]= (double)iage;
y[iage]= log(freq[i][k][iage]/freq[i][i][iage]);
- /* printf("i=%d, k=%d, jk=%d, j1=%d, jj=%d, y[%d]=%f\n",i,k,jk,j1,jj, iage, y[iage]); */
+ /* printf("i=%d, k=%d, s1=%d, j1=%d, jj=%d, y[%d]=%f\n",i,k,s1,j1,jj, iage, y[iage]); */
}
linreg(iagemin,iagemax,&no,x,y,&a,&b,&r, &sa, &sb ); /* y= a+b*x with standard errors */
- pstart[jk]=b;
- pstart[jk-1]=a;
+ pstart[s1]=b;
+ pstart[s1-1]=a;
}else if( j1!=1 && (j1==2 || (log(j1-1.)/log(2.)-(int)(log(j1-1.)/log(2.))) <0.010) && ( TvarsDind[(int)(log(j1-1.)/log(2.))+1]+2+nagesqr == jj) && Dummy[jj-2-nagesqr]==0){ /* We want only if the position, jj, in model corresponds to unique covariate equal to 1 in j1 combination */
printf("j1=%d, jj=%d, (int)(log(j1-1.)/log(2.))+1=%d, TvarsDind[(int)(log(j1-1.)/log(2.))+1]=%d\n",j1, jj,(int)(log(j1-1.)/log(2.))+1,TvarsDind[(int)(log(j1-1.)/log(2.))+1]);
printf("j1=%d, jj=%d, (log(j1-1.)/log(2.))+1=%f, TvarsDind[(int)(log(j1-1.)/log(2.))+1]=%d\n",j1, jj,(log(j1-1.)/log(2.))+1,TvarsDind[(int)(log(j1-1.)/log(2.))+1]);
- pstart[jk]= log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4]));
+ pstart[s1]= log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4]));
printf("%d%d ",i,k);
fprintf(ficlog,"%d%d ",i,k);
- printf("jk=%d,i=%d,k=%d,p[%d]=%12.7f ln((%.0f/%.0f)/(%.0f/%.0f))= %f, OR=%f sd=%f \n",jk,i,k,jk,p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3],freq[i][k][iagemax+4],freq[i][i][iagemax+4], log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4])),(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4]), sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3]+1/freq[i][k][iagemax+4]+1/freq[i][i][iagemax+4]));
+ printf("s1=%d,i=%d,k=%d,p[%d]=%12.7f ln((%.0f/%.0f)/(%.0f/%.0f))= %f, OR=%f sd=%f \n",s1,i,k,s1,p[s1],freq[i][k][iagemax+3],freq[i][i][iagemax+3],freq[i][k][iagemax+4],freq[i][i][iagemax+4], log((freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4])),(freq[i][k][iagemax+3]/freq[i][i][iagemax+3])/(freq[i][k][iagemax+4]/freq[i][i][iagemax+4]), sqrt(1/freq[i][k][iagemax+3]+1/freq[i][i][iagemax+3]+1/freq[i][k][iagemax+4]+1/freq[i][i][iagemax+4]));
}else{ /* Other cases, like quantitative fixed or varying covariates */
;
}
/* printf("%12.7f )", param[i][jj][k]); */
/* fprintf(ficlog,"%12.7f )", param[i][jj][k]); */
- jk++;
+ s1++;
} /* end jj */
} /* end k!= i */
} /* end k */
- } /* end i, jk */
+ } /* end i, s1 */
} /* end j !=0 */
} /* end selected combination of covariate j1 */
if(j==0){ /* We can estimate starting values from the occurences in each case */
printf("#Freqsummary: Starting values for the constants:\n");
fprintf(ficlog,"\n");
- for(i=1,jk=1; i <=nlstate; i++){
+ for(i=1,s1=1; i <=nlstate; i++){
for(k=1; k <=(nlstate+ndeath); k++){
if (k != i) {
printf("%d%d ",i,k);
fprintf(ficlog,"%d%d ",i,k);
for(jj=1; jj <=ncovmodel; jj++){
- pstart[jk]=p[jk]; /* Setting pstart to p values by default */
+ pstart[s1]=p[s1]; /* Setting pstart to p values by default */
if(jj==1){ /* Age has to be done */
- pstart[jk]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]);
- printf("%12.7f ln(%.0f/%.0f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
- fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f ",p[jk],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
+ pstart[s1]= log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]);
+ printf("%12.7f ln(%.0f/%.0f)= %12.7f ",p[s1],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
+ fprintf(ficlog,"%12.7f ln(%.0f/%.0f)= %12.7f ",p[s1],freq[i][k][iagemax+3],freq[i][i][iagemax+3], log(freq[i][k][iagemax+3]/freq[i][i][iagemax+3]));
}
/* printf("%12.7f )", param[i][jj][k]); */
/* fprintf(ficlog,"%12.7f )", param[i][jj][k]); */
- jk++;
+ s1++;
}
printf("\n");
fprintf(ficlog,"\n");
@@ -4755,17 +4854,17 @@ Title=%s
Datafile=%s Firstpass=%d La
}
printf("#Freqsummary\n");
fprintf(ficlog,"\n");
- for(jk=-1; jk <=nlstate+ndeath; jk++){
- for(m=-1; m <=nlstate+ndeath; m++){
- /* param[i]|j][k]= freq[jk][m][iagemax+3] */
- printf(" %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]);
- fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]);
- /* if(freq[jk][m][iage] !=0 ) { /\* minimizing output *\/ */
- /* printf(" %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); */
- /* fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iagemax+3]); */
+ for(s1=-1; s1 <=nlstate+ndeath; s1++){
+ for(s2=-1; s2 <=nlstate+ndeath; s2++){
+ /* param[i]|j][k]= freq[s1][s2][iagemax+3] */
+ printf(" %d%d=%.0f",s1,s2,freq[s1][s2][iagemax+3]);
+ fprintf(ficlog," %d%d=%.0f",s1,s2,freq[s1][s2][iagemax+3]);
+ /* if(freq[s1][s2][iage] !=0 ) { /\* minimizing output *\/ */
+ /* printf(" %d%d=%.0f",s1,s2,freq[s1][s2][iagemax+3]); */
+ /* fprintf(ficlog," %d%d=%.0f",s1,s2,freq[s1][s2][iagemax+3]); */
/* } */
}
- } /* end loop jk */
+ } /* end loop s1 */
printf("\n");
fprintf(ficlog,"\n");
@@ -4899,7 +4998,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]);
}
}
}
@@ -6137,7 +6239,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 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.\n",optionfilehtmcov);
+ fprintf(fichtm,"\n 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(fichtmcov,"Current page is file %s
\n\nMatrix 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.\
@@ -6354,7 +6456,7 @@ To be simple, these graphs help to under
fprintf(ficgp,"\nset parametric;unset label");
fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2);
fprintf(ficgp,"\nset ter svg size 640, 480");
- fprintf(fichtmcov,"\n
Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year-1\
+ fprintf(fichtmcov,"\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, \
@@ -6365,16 +6467,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 */
@@ -6665,7 +6767,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];
@@ -6675,6 +6777,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); */
@@ -6773,7 +6876,34 @@ void printinggnuplot(char fileresu[], ch
if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
else fprintf(ficgp," %%*lf (%%*lf)");
}
- fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1));
+ /* fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1)); */
+
+ 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 );
+ }else{
+ kl=0;
+ for (k=1; k<=cptcoveff; k++){ /* For each combination of covariate */
+ 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 */
+ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */
+ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
+ vlv= nbcode[Tvaraff[k]][lv];
+ kl++;
+ /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */
+ /*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */
+ /*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(k==cptcoveff){
+ fprintf(ficgp,"$%d==%d && $%d==%d)? $%d : 1/0) t 'Observed prevalence in state %d' w l lt 2",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \
+ 2+cptcoveff*2+3*(cpt-1), cpt ); /* 4 or 6 ?*/
+ }else{
+ fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]);
+ kl++;
+ }
+ } /* end covariate */
+ } /* end if no covariate */
+
if(backcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */
/* fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt); */
fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1, nres in 2 to be fixed */
@@ -7178,12 +7308,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 '' ");
@@ -7195,10 +7329,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 */
@@ -7230,8 +7369,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 );
@@ -7451,10 +7596,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 */
@@ -7462,19 +7608,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++)
@@ -7486,77 +7635,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;
- }else{ /* bad */
- for (i=1; i<=nlstate;i++){
- mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod];
- } /* i */
- } /* 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++){
+ 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 */
agemaxgood[cptcod]=age;
- }else{ /* bad */
- for (i=1; i<=nlstate;i++){
- mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
- } /* i */
+ }
+ } /* 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)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+(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(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 */
+
for (age=bage; age<=fage; age++){
/* printf("%d %d ", cptcod, (int)age); */
@@ -7571,40 +7786,44 @@ 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 */
/************** Forecasting ******************/
- void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){
+void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){
/* proj1, year, month, day of starting projection
agemin, agemax range of age
dateprev1 dateprev2 range of dates during which prevalence is computed
anproj2 year of en of projection (same day and month as proj1).
*/
- int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
+ int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
double agec; /* generic age */
double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
double *popeffectif,*popcount;
@@ -7699,11 +7918,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]);
}
@@ -7715,6 +7934,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 */
@@ -7725,134 +7946,143 @@ set ter svg size 640, 480\nunset log y\n
}
/* /\************** Back Forecasting ******************\/ */
-/* void prevbackforecast(char fileres[], double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){ */
-/* /\* back1, year, month, day of starting backection */
-/* agemin, agemax range of age */
-/* dateprev1 dateprev2 range of dates during which prevalence is computed */
-/* anback2 year of en of backection (same day and month as back1). */
-/* *\/ */
-/* int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1; */
-/* double agec; /\* generic age *\/ */
-/* double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; */
-/* double *popeffectif,*popcount; */
-/* double ***p3mat; */
-/* /\* double ***mobaverage; *\/ */
-/* char fileresfb[FILENAMELENGTH]; */
-
-/* agelim=AGESUP; */
-/* /\* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people */
-/* in each health status at the date of interview (if between dateprev1 and dateprev2). */
-/* We still use firstpass and lastpass as another selection. */
-/* *\/ */
-/* /\* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\ *\/ */
-/* /\* firstpass, lastpass, stepm, weightopt, model); *\/ */
-/* prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */
-
-/* strcpy(fileresfb,"FB_"); */
-/* strcat(fileresfb,fileresu); */
-/* if((ficresfb=fopen(fileresfb,"w"))==NULL) { */
-/* printf("Problem with back forecast resultfile: %s\n", fileresfb); */
-/* fprintf(ficlog,"Problem with back forecast resultfile: %s\n", fileresfb); */
-/* } */
-/* printf("Computing back forecasting: result on file '%s', please wait... \n", fileresfb); */
-/* fprintf(ficlog,"Computing back forecasting: result on file '%s', please wait... \n", fileresfb); */
-
-/* if (cptcoveff==0) ncodemax[cptcoveff]=1; */
-
-/* /\* if (mobilav!=0) { *\/ */
-/* /\* mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); *\/ */
-/* /\* if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){ *\/ */
-/* /\* fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); *\/ */
-/* /\* printf(" Error in movingaverage mobilav=%d\n",mobilav); *\/ */
-/* /\* } *\/ */
-/* /\* } *\/ */
-
-/* stepsize=(int) (stepm+YEARM-1)/YEARM; */
-/* if (stepm<=12) stepsize=1; */
-/* if(estepm < stepm){ */
-/* printf ("Problem %d lower than %d\n",estepm, stepm); */
-/* } */
-/* else hstepm=estepm; */
-
-/* hstepm=hstepm/stepm; */
-/* yp1=modf(dateintmean,&yp);/\* extracts integral of datemean in yp and */
-/* fractional in yp1 *\/ */
-/* anprojmean=yp; */
-/* yp2=modf((yp1*12),&yp); */
-/* mprojmean=yp; */
-/* yp1=modf((yp2*30.5),&yp); */
-/* jprojmean=yp; */
-/* if(jprojmean==0) jprojmean=1; */
-/* if(mprojmean==0) jprojmean=1; */
-
-/* i1=cptcoveff; */
-/* if (cptcovn < 1){i1=1;} */
+void prevbackforecast(char fileres[], double ***prevacurrent, double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){
+ /* back1, year, month, day of starting backection
+ agemin, agemax range of age
+ dateprev1 dateprev2 range of dates during which prevalence is computed
+ anback2 year of en of backection (same day and month as back1).
+ */
+ int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
+ double agec; /* generic age */
+ double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
+ double *popeffectif,*popcount;
+ double ***p3mat;
+ /* double ***mobaverage; */
+ char fileresfb[FILENAMELENGTH];
+
+ agelim=AGESUP;
+ /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people
+ in each health status at the date of interview (if between dateprev1 and dateprev2).
+ We still use firstpass and lastpass as another selection.
+ */
+ /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\ */
+ /* firstpass, lastpass, stepm, weightopt, model); */
+
+ /*Do we need to compute prevalence again?*/
+
+ /* prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */
-/* fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2); */
+ strcpy(fileresfb,"FB_");
+ strcat(fileresfb,fileresu);
+ if((ficresfb=fopen(fileresfb,"w"))==NULL) {
+ printf("Problem with back forecast resultfile: %s\n", fileresfb);
+ fprintf(ficlog,"Problem with back forecast resultfile: %s\n", fileresfb);
+ }
+ printf("\nComputing back forecasting: result on file '%s', please wait... \n", fileresfb);
+ fprintf(ficlog,"\nComputing back forecasting: result on file '%s', please wait... \n", fileresfb);
-/* fprintf(ficresfb,"#****** Routine prevbackforecast **\n"); */
-
-/* /\* if (h==(int)(YEARM*yearp)){ *\/ */
-/* for(cptcov=1, k=0;cptcov<=i1;cptcov++){ */
-/* for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ */
-/* k=k+1; */
-/* fprintf(ficresfb,"\n#****** hbijx=probability over h years, hp.jx is weighted by observed prev \n#"); */
-/* for(j=1;j<=cptcoveff;j++) { */
-/* fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */
-/* } */
-/* fprintf(ficresfb," yearbproj age"); */
-/* for(j=1; j<=nlstate+ndeath;j++){ */
-/* for(i=1; i<=nlstate;i++) */
-/* fprintf(ficresfb," p%d%d",i,j); */
-/* fprintf(ficresfb," p.%d",j); */
-/* } */
-/* for (yearp=0; yearp>=(anback2-anback1);yearp -=stepsize) { */
-/* /\* for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { *\/ */
-/* fprintf(ficresfb,"\n"); */
-/* fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp); */
-/* for (agec=fage; agec>=(ageminpar-1); agec--){ */
-/* nhstepm=(int) rint((agelim-agec)*YEARM/stepm); */
-/* nhstepm = nhstepm/hstepm; */
-/* p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); */
-/* oldm=oldms;savm=savms; */
-/* hbxij(p3mat,nhstepm,agec,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm,oldm,savm, dnewm, doldm, dsavm, k); */
-/* for (h=0; h<=nhstepm; h++){ */
-/* if (h*hstepm/YEARM*stepm ==yearp) { */
-/* fprintf(ficresfb,"\n"); */
-/* for(j=1;j<=cptcoveff;j++) */
-/* fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */
-/* fprintf(ficresfb,"%.f %.f ",anback1+yearp,agec+h*hstepm/YEARM*stepm); */
-/* } */
-/* for(j=1; j<=nlstate+ndeath;j++) { */
-/* ppij=0.; */
-/* for(i=1; i<=nlstate;i++) { */
-/* if (mobilav==1) */
-/* ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][cptcod]; */
-/* else { */
-/* ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod]; */
-/* } */
-/* if (h*hstepm/YEARM*stepm== yearp) { */
-/* fprintf(ficresfb," %.3f", p3mat[i][j][h]); */
-/* } */
-/* } /\* end i *\/ */
-/* if (h*hstepm/YEARM*stepm==yearp) { */
-/* fprintf(ficresfb," %.3f", ppij); */
-/* } */
-/* }/\* end j *\/ */
-/* } /\* end h *\/ */
-/* free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); */
-/* } /\* end agec *\/ */
-/* } /\* end yearp *\/ */
-/* } /\* end cptcod *\/ */
-/* } /\* end cptcov *\/ */
-
-/* /\* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); *\/ */
-
-/* fclose(ficresfb); */
-/* printf("End of Computing Back forecasting \n"); */
-/* fprintf(ficlog,"End of Computing Back forecasting\n"); */
+ if (cptcoveff==0) ncodemax[cptcoveff]=1;
+
+
+ stepsize=(int) (stepm+YEARM-1)/YEARM;
+ if (stepm<=12) stepsize=1;
+ if(estepm < stepm){
+ printf ("Problem %d lower than %d\n",estepm, stepm);
+ }
+ else hstepm=estepm;
+
+ hstepm=hstepm/stepm;
+ yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp and
+ fractional in yp1 */
+ anprojmean=yp;
+ yp2=modf((yp1*12),&yp);
+ mprojmean=yp;
+ yp1=modf((yp2*30.5),&yp);
+ jprojmean=yp;
+ if(jprojmean==0) jprojmean=1;
+ if(mprojmean==0) jprojmean=1;
+
+ i1=pow(2,cptcoveff);
+ if (cptcovn < 1){i1=1;}
+
+ fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2);
+
+ fprintf(ficresfb,"#****** Routine prevbackforecast **\n");
+
+ /* if (h==(int)(YEARM*yearp)){ */
+ /* for(cptcov=1, k=0;cptcov<=i1;cptcov++){ */
+ /* for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ */
+ /* k=k+1; */
+ for(nres=1; nres <= nresult; nres++) /* For each resultline */
+ for(k=1; k<=i1;k++){
+ if(i1 != 1 && TKresult[nres]!= k)
+ continue;
+ if(invalidvarcomb[k]){
+ printf("\nCombination (%d) projection ignored because no cases \n",k);
+ continue;
+ }
+ fprintf(ficresfb,"\n#****** hbijx=probability over h years, hp.jx is weighted by observed prev \n#");
+ for(j=1;j<=cptcoveff;j++) {
+ fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ }
+ for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+ fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+ }
+ fprintf(ficresfb," yearbproj age");
+ for(j=1; j<=nlstate+ndeath;j++){
+ for(i=1; i<=nlstate;i++)
+ fprintf(ficresfb," p%d%d",i,j);
+ fprintf(ficresfb," p.%d",j);
+ }
+ for (yearp=0; yearp>=(anback2-anback1);yearp -=stepsize) {
+ /* for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { */
+ fprintf(ficresfb,"\n");
+ fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp);
+ /* for (agec=fage; agec>=(ageminpar-1); agec--){ */
+ /* nhstepm=(int) rint((agelim-agec)*YEARM/stepm); */
+ for (agec=fage; agec>=fage-20; agec--){ /* testing up to 10 */
+ nhstepm=(int) rint((agelim-agec)*YEARM/stepm);
+ nhstepm = nhstepm/hstepm;
+ p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ oldm=oldms;savm=savms;
+ hbxij(p3mat,nhstepm,agec,hstepm,p,prevacurrent,nlstate,stepm, k, nres);
+
+ for (h=0; h<=nhstepm; h++){
+ if (h*hstepm/YEARM*stepm ==yearp) {
+ fprintf(ficresfb,"\n");
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficresfb,"%.f %.f ",anback1+yearp,agec+h*hstepm/YEARM*stepm);
+ }
+ for(j=1; j<=nlstate+ndeath;j++) {
+ ppij=0.;
+ for(i=1; i<=nlstate;i++) {
+ /* 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]; */
+ /* } */
+ if (h*hstepm/YEARM*stepm== yearp) {
+ fprintf(ficresfb," %.3f", p3mat[i][j][h]);
+ }
+ } /* end i */
+ if (h*hstepm/YEARM*stepm==yearp) {
+ fprintf(ficresfb," %.3f", ppij);
+ }
+ }/* end j */
+ } /* end h */
+ free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ } /* end agec */
+ } /* end yearp */
+ } /* end k */
+
+ /* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
+
+ fclose(ficresfb);
+ printf("End of Computing Back forecasting \n");
+ fprintf(ficlog,"End of Computing Back forecasting\n");
-/* } */
+}
/************** Forecasting *****not tested NB*************/
/* void populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2s, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){ */
@@ -9696,6 +9926,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++)
@@ -9854,10 +10086,12 @@ 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);
+ hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k, nres);
/* hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm, dnewm, doldm, dsavm, k); */
fprintf(ficrespijb,"# Cov Agex agex-h hbijx with i,j=");
for(i=1; i<=nlstate;i++)
@@ -11417,7 +11651,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, \
@@ -11479,8 +11713,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");
@@ -11512,8 +11752,8 @@ Please run with mle=-1 to get a correct
fclose(ficrespijb);
free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */
- /* prevbackforecast(fileresu, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj,
- bage, fage, firstpass, lastpass, anback2, p, cptcoveff); */
+ prevbackforecast(fileresu, mobaverage, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj,
+ bage, fage, firstpass, lastpass, anback2, p, cptcoveff);
free_matrix(ddnewms, 1, nlstate+ndeath, 1, nlstate+ndeath);
free_matrix(ddsavms, 1, nlstate+ndeath, 1, nlstate+ndeath);
free_matrix(ddoldms, 1, nlstate+ndeath, 1, nlstate+ndeath);