/* $Id$
$State$
$Log$
+ Revision 1.218 2016/02/12 11:29:23 brouard
+ Summary: 0.99 Back projections
+
Revision 1.217 2015/12/23 17:18:31 brouard
Summary: Experimental backcast
undefined. Usually 3: -1, 0 and 1. */
double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;
double **pmmij, ***probs; /* Global pointer */
-double ***mobaverage; /* New global variable */
+double ***mobaverage, ***mobaverages; /* New global variable */
double *ageexmed,*agecens;
double dateintmean=0;
if (cptcovn<1) {j=1;ncodemax[1]=1;}
first=1;
- for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){
+ 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+3+AGEMARGE; iage++)
prop[i][iage]=0.0;
for (i=1; i<=imx; i++) { /* Each individual */
bool=1;
if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
- for (z1=1; z1<=cptcoveff; z1++)
+ for (z1=1; z1<=cptcoveff; z1++) /* For each covariate, look at the value for individual i and checks if it is equal to the corresponding value of this covariate according to current combination j1*/
if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)])
bool=0;
}
- if (bool==1) {
+ if (bool==1) { /* For this combination of covariates values, this individual fits */
/* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */
for(mi=1; mi<wav[i];mi++){
m=mw[mi][i];
for (j=1; j<=(cptcovs); j++) { /* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only */
for (k=-1; k < maxncov; k++) Ndum[k]=0;
for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the
- modality of this covariate Vj*/
+ modality of this covariate Vj*/
ij=(int)(covar[Tvar[j]][i]); /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i
- * If product of Vn*Vm, still boolean *:
- * If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables
- * 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0 */
+ * If product of Vn*Vm, still boolean *:
+ * If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables
+ * 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0 */
/* Finds for covariate j, n=Tvar[j] of Vn . ij is the
modality of the nth covariate of individual i. */
if (ij > modmaxcovj)
modmaxcovj=ij;
else if (ij < modmincovj)
- modmincovj=ij;
+ modmincovj=ij;
if ((ij < -1) && (ij > NCOVMAX)){
- printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX );
- exit(1);
+ printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX );
+ exit(1);
}else
Ndum[ij]++; /*counts and stores the occurence of this modality 0, 1, -1*/
/* If coded 1, 2, 3 , counts the number of 1 Ndum[1], number of 2, Ndum[2], etc */
printf("Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], k, Ndum[k]);
fprintf(ficlog, "Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], k, Ndum[k]);
if( Ndum[k] != 0 ){ /* Counts if nobody answered modality k ie empty modality, we skip it and reorder */
- if( k != -1){
- ncodemax[j]++; /* ncodemax[j]= Number of modalities of the j th
- covariate for which somebody answered excluding
- undefined. Usually 2: 0 and 1. */
- }
- ncodemaxwundef[j]++; /* ncodemax[j]= Number of modalities of the j th
- covariate for which somebody answered including
- undefined. Usually 3: -1, 0 and 1. */
+ if( k != -1){
+ ncodemax[j]++; /* ncodemax[j]= Number of modalities of the j th
+ covariate for which somebody answered excluding
+ undefined. Usually 2: 0 and 1. */
+ }
+ ncodemaxwundef[j]++; /* ncodemax[j]= Number of modalities of the j th
+ covariate for which somebody answered including
+ undefined. Usually 3: -1, 0 and 1. */
}
/* In fact ncodemax[j]=2 (dichotom. variables only) but it could be more for
- historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */
+ historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */
} /* Ndum[-1] number of undefined modalities */
-
+
/* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */
/* For covariate j, modalities could be 1, 2, 3, 4, 5, 6, 7.
If Ndum[1]=0, Ndum[2]=0, Ndum[3]= 635, Ndum[4]=0, Ndum[5]=0, Ndum[6]=27, Ndum[7]=125;
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*/
if (Ndum[i] == 0) { /* If nobody responded to this modality k */
- break;
- }
+ break;
+ }
ij++;
nbcode[Tvar[j]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/
cptcode = ij; /* New max modality for covar j */
/* } /\* 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*/
- for (k=-1; k< maxncov; k++) Ndum[k]=0;
+ for (k=-1; k< maxncov; k++) Ndum[k]=0;
for (i=1; i<=ncovmodel-2-nagesqr; i++) { /* -2, cste and age and eventually age*age */
- /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/
- ij=Tvar[i]; /* Tvar might be -1 if status was unknown */
- Ndum[ij]++; /* Might be supersed V1 + V1*age */
- }
-
- ij=0;
- for (i=0; i<= maxncov-1; i++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */
- /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/
- if((Ndum[i]!=0) && (i<=ncovcol)){
- ij++;
- /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/
- Tvaraff[ij]=i; /*For printing (unclear) */
- }else{
- /* Tvaraff[ij]=0; */
- }
- }
- /* ij--; */
- cptcoveff=ij; /*Number of total covariates*/
-
+ /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/
+ ij=Tvar[i]; /* Tvar might be -1 if status was unknown */
+ Ndum[ij]++; /* Might be supersed V1 + V1*age */
+ }
+
+ ij=0;
+ for (i=0; i<= maxncov-1; i++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */
+ /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/
+ if((Ndum[i]!=0) && (i<=ncovcol)){
+ ij++;
+ /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/
+ Tvaraff[ij]=i; /*For printing (unclear) */
+ }else{
+ /* Tvaraff[ij]=0; */
+ }
+ }
+ /* ij--; */
+ cptcoveff=ij; /*Number of total covariates*/
+
}
int lv=0, vlv=0, kl=0;
int ng=0;
int vpopbased;
+ int ioffset; /* variable offset for columns */
+
/* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */
/* printf("Problem with file %s",optionfilegnuplot); */
/* fprintf(ficlog,"Problem with file %s",optionfilegnuplot); */
fprintf(ficgp,"unset log;\n# plot weighted, mean weight should have point size of 0.5\n plot \"%s\"",subdirf(fileresilk));
fprintf(ficgp," u 2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable \\\n",i,1,i,1);
for (j=2; j<= nlstate+ndeath ; j ++) {
- fprintf(ficgp,",\\\n \"\" u 2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j);
+ fprintf(ficgp,",\\\n \"\" u 2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j);
}
fprintf(ficgp,";\nset out; unset ylabel;\n");
}
for (k1=1; k1<= m ; k1 ++) { /* For each combination of covariate */
/* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files ");
- 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 */
- /* 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[lv]][lv];
- fprintf(ficgp," V%d=%d ",k,vlv);
+ 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 */
+ /* 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]; /* vlv is the value of the covariate lv, 0 or 1 */
+ /* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */
+ fprintf(ficgp," V%d=%d ",k,vlv);
}
fprintf(ficgp,"\n#\n");
- fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1);
- fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);
- fprintf(ficgp,"set xlabel \"Age\" \n\
-set ylabel \"Probability\" \n\
-set ter svg size 640, 480\n\
+ fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1);
+ fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);
+ fprintf(ficgp,"set xlabel \"Age\" \n\
+set ylabel \"Probability\" \n \
+set ter svg size 640, 480\n \
plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1);
-
- for (i=1; i<= nlstate ; i ++) {
- 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+1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
- 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-1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
- 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,\"%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));
- if(backcast==1){
- 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,"\nset out \n");
+
+ for (i=1; i<= nlstate ; i ++) {
+ 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+1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
+ 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-1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);
+ 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,\"%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));
+ 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 */
+ 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 'Backward prevalence in state %d' with line ",kl+1, k,kl+1+1,nbcode[Tvaraff[k]][lv], \
+ 4+(cpt-1), cpt );
+ }else{
+ fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, k,kl+1+1,nbcode[Tvaraff[k]][lv]);
+ kl++;
+ }
+ } /* end covariate */
+ }
+ fprintf(ficgp,"\nset out \n");
} /* k1 */
} /* cpt */
/*2 eme*/
for (k1=1; k1<= m ; k1 ++) {
fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files ");
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 */
- /* 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[lv]][lv];
- fprintf(ficgp," V%d=%d ",k,vlv);
+ 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 */
+ /* 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];
+ fprintf(ficgp," V%d=%d ",k,vlv);
}
fprintf(ficgp,"\n#\n");
-
- fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1);
- for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
- if(vpopbased==0)
- fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
- else
- fprintf(ficgp,"\nreplot ");
- for (i=1; i<= nlstate+1 ; i ++) {
- k=2*i;
- fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1, vpopbased);
- for (j=1; j<= nlstate+1 ; j ++) {
- if (j==i) fprintf(ficgp," %%lf (%%lf)");
- else fprintf(ficgp," %%*lf (%%*lf)");
- }
- if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i);
- else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1);
- fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
- for (j=1; j<= nlstate+1 ; j ++) {
- if (j==i) fprintf(ficgp," %%lf (%%lf)");
- else fprintf(ficgp," %%*lf (%%*lf)");
- }
- fprintf(ficgp,"\" t\"\" w l lt 0,");
- fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4+$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
- for (j=1; j<= nlstate+1 ; j ++) {
- if (j==i) fprintf(ficgp," %%lf (%%lf)");
- else fprintf(ficgp," %%*lf (%%*lf)");
- }
- if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
- else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
- } /* state */
- } /* vpopbased */
- fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */
+
+ fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1);
+ for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
+ if(vpopbased==0)
+ fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
+ else
+ fprintf(ficgp,"\nreplot ");
+ for (i=1; i<= nlstate+1 ; i ++) {
+ k=2*i;
+ fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1, vpopbased);
+ for (j=1; j<= nlstate+1 ; j ++) {
+ if (j==i) fprintf(ficgp," %%lf (%%lf)");
+ else fprintf(ficgp," %%*lf (%%*lf)");
+ }
+ if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i);
+ else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1);
+ fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
+ for (j=1; j<= nlstate+1 ; j ++) {
+ if (j==i) fprintf(ficgp," %%lf (%%lf)");
+ else fprintf(ficgp," %%*lf (%%*lf)");
+ }
+ fprintf(ficgp,"\" t\"\" w l lt 0,");
+ fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4+$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),k1-1,k1-1,vpopbased);
+ for (j=1; j<= nlstate+1 ; j ++) {
+ if (j==i) fprintf(ficgp," %%lf (%%lf)");
+ else fprintf(ficgp," %%*lf (%%*lf)");
+ }
+ if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
+ else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
+ } /* state */
+ } /* vpopbased */
+ fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */
} /* k1 */
-
-
+
+
/*3eme*/
for (k1=1; k1<= m ; k1 ++) {
for (cpt=1; cpt<= nlstate ; cpt ++) {
fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files: cov=%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 */
- /* 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[lv]][lv];
- fprintf(ficgp," V%d=%d ",k,vlv);
+ 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 */
+ /* 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];
+ fprintf(ficgp," V%d=%d ",k,vlv);
}
fprintf(ficgp,"\n#\n");
-
+
/* k=2+nlstate*(2*cpt-2); */
k=2+(nlstate+1)*(cpt-1);
fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1);
fprintf(ficgp,"set ter svg size 640, 480\n\
plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileresu,"E_"),k1-1,k1-1,k,cpt);
/*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
- for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
- fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
- fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
- for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
- fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
-
+ for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
+ fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
+ fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d+2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
+ for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
+ fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
+
*/
for (i=1; i< nlstate ; i ++) {
- fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+i,cpt,i+1);
- /* fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/
-
+ fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+i,cpt,i+1);
+ /* fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d%d\" w l",subdirf2(fileres,"e"),k1-1,k1-1,k+2*i,cpt,i+1);*/
+
}
fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+nlstate,cpt);
}
/* 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[lv]][lv];
+ vlv= nbcode[Tvaraff[k]][lv];
fprintf(ficgp," V%d=%d ",k,vlv);
}
fprintf(ficgp,"\n#\n");
/* 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[lv]][lv];
+ vlv= nbcode[Tvaraff[k]][lv];
fprintf(ficgp," V%d=%d ",k,vlv);
}
fprintf(ficgp,"\n#\n");
/* 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[lv]][lv];
+ vlv= nbcode[Tvaraff[k]][lv];
fprintf(ficgp," V%d=%d ",k,vlv);
}
fprintf(ficgp,"\n#\n");
/* 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[lv]][lv];
+ vlv= nbcode[Tvaraff[k]][lv];
fprintf(ficgp," V%d=%d ",k,vlv);
}
fprintf(ficgp,"\n#\n");
for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
- fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (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 */
- /* 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[lv]][lv];
- fprintf(ficgp," V%d=%d ",k,vlv);
- }
- fprintf(ficgp,"\n#\n");
-
- fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n ");
- fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1);
- fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
-set ter svg size 640, 480\n\
-unset log y\n\
+ fprintf(ficgp,"\n#\n#\n#Projection of prevalence to stable (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 */
+ /* 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];
+ fprintf(ficgp," V%d=%d ",k,vlv);
+ }
+ fprintf(ficgp,"\n#\n");
+
+ fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n ");
+ fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1);
+ fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
+set ter svg size 640, 480\n \
+unset log y\n \
plot [%.f:%.f] ", ageminpar, agemaxpar);
- for (i=1; 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){
- fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"F_"));
- }else{
- fprintf(ficgp,",\\\n '' ");
- }
- if(cptcoveff ==0){ /* No covariate */
- fprintf(ficgp," u 2:("); /* Age is in 2 */
- /*# 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 */
- if(i==nlstate+1)
- fprintf(ficgp," $%d/(1.-$%d)) t 'p.%d' with line ", \
- 2+(cpt-1)*(nlstate+1)+1+(i-1), 2+1+(i-1)+(nlstate+1)*nlstate,cpt );
- else
- fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \
- 2+(cpt-1)*(nlstate+1)+1+(i-1), 2+1+(i-1)+(nlstate+1)*nlstate,i,cpt );
- }else{
- fprintf(ficgp,"u 6:(("); /* Age is in 6 */
- /*# 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 */
- kl=0;
- for (k=1; k<=cptcoveff; k++){ /* For each 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[lv]][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)
- if(i==nlstate+1)
- fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv], \
- 6+(cpt-1)*(nlstate+1)+1+(i-1), 6+1+(i-1)+(nlstate+1)*nlstate,cpt );
- else
- fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv], \
- 6+(cpt-1)*(nlstate+1)+1+(i-1), 6+1+(i-1)+(nlstate+1)*nlstate,i,cpt );
- else{
- fprintf(ficgp,"$%d==%d && $%d==%d && ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv]);
- kl++;
- }
- } /* end covariate */
- } /* end if covariate */
- } /* nlstate */
- fprintf(ficgp,"\nset out\n");
- } /* end cpt state*/
- } /* end covariate */
- } /* End if prevfcast */
-
-
- /* proba elementaires */
- fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n");
+ for (i=1; 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){
+ fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"F_"));
+ }else{
+ fprintf(ficgp,",\\\n '' ");
+ }
+ if(cptcoveff ==0){ /* No covariate */
+ ioffset=2; /* Age is in 2 */
+ /*# 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 */
+ /*# 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 ", \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
+ 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 */
+ if(cptcoveff ==1){
+ ioffset=4; /* Age is in 4 */
+ }else{
+ ioffset=6; /* Age is in 6 */
+ /*# 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 */
+ }
+ fprintf(ficgp," u %d:((",ioffset);
+ kl=0;
+ for (k=1; k<=cptcoveff; k++){ /* For each covariate */
+ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to combination k1 and covariate k */
+ /* 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){
+ if(i==nlstate+1){
+ if(cptcoveff ==1){
+ fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ",kl, k,kl+1,nbcode[Tvaraff[k]][lv], \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
+ }else{
+ fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ",kl, k,kl+1,nbcode[Tvaraff[k]][lv], \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
+ }
+ }else{
+ if(cptcoveff ==1){
+ fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ",kl, k,kl+1,nbcode[Tvaraff[k]][lv], \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset +1+(i-1)+(nlstate+1)*nlstate,i,cpt );
+ }else{
+ fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ",kl, k,kl+1,nbcode[Tvaraff[k]][lv], \
+ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset +1+(i-1)+(nlstate+1)*nlstate,i,cpt );
+ }
+ }
+ }else{ /* k < cptcoveff */
+ fprintf(ficgp,"$%d==%d && $%d==%d && ",kl, k,kl+1,nbcode[Tvaraff[k]][lv]);
+ kl++;
+ }
+ } /* end covariate */
+ } /* end if covariate */
+ } /* nlstate */
+ fprintf(ficgp,"\nset out\n");
+ } /* end cpt state*/
+ } /* end covariate */
+ } /* End if prevfcast */
+
+
+ /* proba elementaires */
+ fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n");
for(i=1,jk=1; i <=nlstate; i++){
fprintf(ficgp,"# initial state %d\n",i);
for(k=1; k <=(nlstate+ndeath); k++){
/*************** Moving average **************/
- /* int movingaverage(double ***probs, double bage, double fage, double ***mobaverage, int mobilav, double bageout, double fageout){ */
+/* int movingaverage(double ***probs, double bage, double fage, double ***mobaverage, int mobilav, double bageout, double fageout){ */
int movingaverage(double ***probs, double bage, double fage, double ***mobaverage, int mobilav){
int i, cpt, cptcod;
int modcovmax =1;
int mobilavrange, mob;
- double age;
int iage=0;
+
+ double sum=0.;
+ double age;
double *sumnewp, *sumnewm;
double *agemingood, *agemaxgood; /* Currently identical for all covariates */
+
+ modcovmax=2*cptcoveff;/* Max number of modalities. We suppose
+ a covariate has 2 modalities, should be equal to ncovcombmax */
+
sumnewp = vector(1,modcovmax);
sumnewm = vector(1,modcovmax);
agemingood = vector(1,modcovmax);
agemaxgood = vector(1,modcovmax);
-
-
- modcovmax=2*cptcoveff;/* Max number of modalities. We suppose
- a covariate has 2 modalities, should be equal to ncovcombmax */
+
+ for (cptcod=1;cptcod<=modcovmax;cptcod++){
+ sumnewm[cptcod]=0.;
+ sumnewp[cptcod]=0.;
+ agemingood[cptcod]=0;
+ agemaxgood[cptcod]=0;
+ }
if (cptcovn<1) modcovmax=1; /* At least 1 pass */
if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){
else mobilavrange=mobilav;
for (age=bage; age<=fage; age++)
for (i=1; i<=nlstate;i++)
- for (cptcod=1;cptcod<=modcovmax;cptcod++)
- mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod];
+ for (cptcod=1;cptcod<=modcovmax;cptcod++)
+ mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod];
/* We keep the original values on the extreme ages bage, fage and for
fage+1 and bage-1 we use a 3 terms moving average; for fage+2 bage+2
we use a 5 terms etc. until the borders are no more concerned.
*/
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<=modcovmax;cptcod++){
- 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;
- }
- }
+ for (i=1; i<=nlstate;i++){
+ for (cptcod=1;cptcod<=modcovmax;cptcod++){
+ 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;
+ }
+ }
}/* end age */
}/* end mob */
}else
return -1;
for (cptcod=1;cptcod<=modcovmax;cptcod++){
/* for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ */
- agemingood[cptcod]=fage+(mob-1)/2;
+ agemingood[cptcod]=fage-(mob-1)/2;
for (age=fage-(mob-1)/2; age>=bage; age--){/* From oldest to youngest, finding the youngest wrong */
sumnewm[cptcod]=0.;
for (i=1; i<=nlstate;i++){
sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
}
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 */
- /* From youngest, finding the oldest wrong */
- agemaxgood[cptcod]=bage+(mob-1)/2;
- for (age=bage+(mob-1)/2; age<=fage; age++){
- sumnewm[cptcod]=0.;
- for (i=1; i<=nlstate;i++){
- sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
- }
- if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */
- agemaxgood[cptcod]=age;
+ agemingood[cptcod]=age;
}else{ /* bad */
- for (i=1; i<=nlstate;i++){
- mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod];
- } /* i */
+ for (i=1; i<=nlstate;i++){
+ mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod];
+ } /* i */
} /* end bad */
}/* age */
- for (age=bage; age<=fage; age++){
- printf("%d %d ", cptcod, (int)age);
- sumnewp[cptcod]=0.;
- sumnewm[cptcod]=0.;
- for (i=1; i<=nlstate;i++){
- sumnewp[cptcod]+=probs[(int)age][i][cptcod];
- sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
- printf("%.4f %.4f ",probs[(int)age][i][cptcod], mobaverage[(int)age][i][cptcod]);
- }
- printf("%.4f %.4f \n",sumnewp[cptcod], sumnewm[cptcod]);
+ sum=0.;
+ for (i=1; i<=nlstate;i++){
+ sum+=mobaverage[(int)agemingood[cptcod]][i][cptcod];
}
- printf("\n");
+ 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++){
+ sumnewm[cptcod]=0.;
+ for (i=1; i<=nlstate;i++){
+ sumnewm[cptcod]+=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 */
+ } /* 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);
+ sumnewp[cptcod]=0.;
+ sumnewm[cptcod]=0.;
+ for (i=1; i<=nlstate;i++){
+ sumnewp[cptcod]+=probs[(int)age][i][cptcod];
+ sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod];
+ /* printf("%.4f %.4f ",probs[(int)age][i][cptcod], mobaverage[(int)age][i][cptcod]); */
+ }
+ /* printf("%.4f %.4f \n",sumnewp[cptcod], sumnewm[cptcod]); */
+ }
+ /* 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]);
+ 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]);
+ 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.;
+ /*printf("i=%d, age=%d, cptcod=%d\n",i, (int)age, cptcod);*/
+ mobaverage[(int)age][i][cptcod]=0.;
}
}
}/* end cptcod */
k=k+1;
fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#");
for(j=1;j<=cptcoveff;j++) {
- fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
}
fprintf(ficresf," yearproj age");
for(j=1; j<=nlstate+ndeath;j++){
- for(i=1; i<=nlstate;i++)
+ for(i=1; i<=nlstate;i++)
fprintf(ficresf," p%d%d",i,j);
- fprintf(ficresf," p.%d",j);
+ fprintf(ficresf," wp.%d",j);
}
for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {
- fprintf(ficresf,"\n");
- fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+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;
- hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);
-
- for (h=0; h<=nhstepm; h++){
- if (h*hstepm/YEARM*stepm ==yearp) {
+ fprintf(ficresf,"\n");
+ fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+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;
+ hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);
+
+ for (h=0; h<=nhstepm; h++){
+ if (h*hstepm/YEARM*stepm ==yearp) {
fprintf(ficresf,"\n");
for(j=1;j<=cptcoveff;j++)
fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- fprintf(ficresf,"%.f %.f ",anproj1+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(ficresf," %.3f", p3mat[i][j][h]);
- }
- } /* end i */
- if (h*hstepm/YEARM*stepm==yearp) {
- fprintf(ficresf," %.3f", ppij);
- }
- }/* end j */
- } /* end h */
- free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
- } /* end agec */
+ fprintf(ficresf,"%.f %.f ",anproj1+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(ficresf," %.3f", p3mat[i][j][h]);
+ }
+ } /* end i */
+ if (h*hstepm/YEARM*stepm==yearp) {
+ fprintf(ficresf," %.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 */
-
+
fclose(ficresf);
printf("End of Computing forecasting \n");
fprintf(ficlog,"End of Computing forecasting\n");
#endif
- }
+}
- int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp){
+int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp){
/*--------------- Prevalence limit (period or stable prevalence) --------------*/
int i, j, k, i1 ;
/* double ftolpl = 1.e-10; */
for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
fprintf(ficrespl,"\n");
- /* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */
+ /* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */
- agebase=ageminpar;
- agelim=agemaxpar;
+ agebase=ageminpar;
+ agelim=agemaxpar;
- i1=pow(2,cptcoveff);
- if (cptcovn < 1){i1=1;}
+ i1=pow(2,cptcoveff);
+ if (cptcovn < 1){i1=1;}
- for(cptcov=1,k=0;cptcov<=i1;cptcov++){
+ for(cptcov=1,k=0;cptcov<=i1;cptcov++){
/* for(cptcov=1,k=0;cptcov<=1;cptcov++){ */
- //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
- k=k+1;
- /* to clean */
- //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
- fprintf(ficrespl,"#******");
- printf("#******");
- fprintf(ficlog,"#******");
- for(j=1;j<=cptcoveff;j++) {
- fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- }
- fprintf(ficrespl,"******\n");
- printf("******\n");
- fprintf(ficlog,"******\n");
+ //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
+ k=k+1;
+ /* to clean */
+ //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
+ fprintf(ficrespl,"#******");
+ printf("#******");
+ fprintf(ficlog,"#******");
+ for(j=1;j<=cptcoveff;j++) {
+ fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ }
+ fprintf(ficrespl,"******\n");
+ printf("******\n");
+ fprintf(ficlog,"******\n");
- fprintf(ficrespl,"#Age ");
- for(j=1;j<=cptcoveff;j++) {
- fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- }
- for(i=1; i<=nlstate;i++) fprintf(ficrespl," %d-%d ",i,i);
- fprintf(ficrespl,"Total Years_to_converge\n");
+ fprintf(ficrespl,"#Age ");
+ for(j=1;j<=cptcoveff;j++) {
+ fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ }
+ for(i=1; i<=nlstate;i++) fprintf(ficrespl," %d-%d ",i,i);
+ fprintf(ficrespl,"Total Years_to_converge\n");
- for (age=agebase; age<=agelim; age++){
- /* for (age=agebase; age<=agebase; age++){ */
- prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k);
- fprintf(ficrespl,"%.0f ",age );
- for(j=1;j<=cptcoveff;j++)
- fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- tot=0.;
- for(i=1; i<=nlstate;i++){
- tot += prlim[i][i];
- fprintf(ficrespl," %.5f", prlim[i][i]);
- }
- fprintf(ficrespl," %.3f %d\n", tot, *ncvyearp);
- } /* Age */
- /* was end of cptcod */
- } /* cptcov */
- return 0;
+ for (age=agebase; age<=agelim; age++){
+ /* for (age=agebase; age<=agebase; age++){ */
+ prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k);
+ fprintf(ficrespl,"%.0f ",age );
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ tot=0.;
+ for(i=1; i<=nlstate;i++){
+ tot += prlim[i][i];
+ fprintf(ficrespl," %.5f", prlim[i][i]);
+ }
+ fprintf(ficrespl," %.3f %d\n", tot, *ncvyearp);
+ } /* Age */
+ /* was end of cptcod */
+ } /* cptcov */
+ return 0;
}
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){
if(mobilavproj > 0){
/* bprevalim(bprlim, mobaverage, nlstate, p, age, ageminpar, agemaxpar, oldm, savm, doldm, dsavm, ftolpl, ncvyearp, k); */
/* bprevalim(bprlim, mobaverage, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
- bprevalim(bprlim, mobaverage, nlstate, p, age, ftolpl, ncvyearp, k);
+ bprevalim(bprlim, mobaverage, nlstate, p, age, ftolpl, ncvyearp, k);
}else if (mobilavproj == 0){
- printf("There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
- fprintf(ficlog,"There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
- exit(1);
+ printf("There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
+ fprintf(ficlog,"There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
+ exit(1);
}else{
- /* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
- bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k);
+ /* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
+ bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k);
}
fprintf(ficresplb,"%.0f ",age );
for(j=1;j<=cptcoveff;j++)
- fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
tot=0.;
for(i=1; i<=nlstate;i++){
- tot += bprlim[i][i];
- fprintf(ficresplb," %.5f", bprlim[i][i]);
+ tot += bprlim[i][i];
+ fprintf(ficresplb," %.5f", bprlim[i][i]);
}
fprintf(ficresplb," %.3f %d\n", tot, *ncvyearp);
} /* Age */
double agedeb=0.;
double ageminpar=AGEOVERFLOW,agemin=AGEOVERFLOW, agemaxpar=-AGEOVERFLOW, agemax=-AGEOVERFLOW;
- double ageminout=-AGEOVERFLOW,agemaxout=AGEOVERFLOW; /* Smaller Age range redefined after movingaverage */
+ double ageminout=-AGEOVERFLOW,agemaxout=AGEOVERFLOW; /* Smaller Age range redefined after movingaverage */
double fret;
double dum=0.; /* Dummy variable */
* bbbbbbbb
* 76543210
* h-1 00000101 (6-1=5)
- *(h-1)>>(k-1)= 00000001 >> (2-1) = 1 right shift
+ *(h-1)>>(k-1)= 00000010 >> (2-1) = 1 right shift
* &
* 1 00000001 (1)
- * 00000001 = 1 & ((h-1) >> (k-1))
- * +1= 00000010 =2
+ * 00000000 = 1 & ((h-1) >> (k-1))
+ * +1= 00000001 =1
*
* h=14, k=3 => h'=h-1=13, k'=k-1=2
* h' 1101 =2^3+2^2+0x2^1+2^0
if((num_filled=sscanf(line,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%lf\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm, &ftolpl)) !=EOF){
if (num_filled != 6) {
- printf("Not 6 parameters in line, for example:agemin=60 agemax=95 bage=55 fage=95 estepm=24 ftolpl=6e-4\n");
- printf("but line=%s\n",line);
+ printf("Error: Not 6 parameters in line, for example:agemin=60 agemax=95 bage=55 fage=95 estepm=24 ftolpl=6e-4\n, your line=%s . Probably you are running an older format.\n",line);
+ fprintf(ficlog,"Error: Not 6 parameters in line, for example:agemin=60 agemax=95 bage=55 fage=95 estepm=24 ftolpl=6e-4\n, your line=%s . Probably you are running an older format.\n",line);
goto end;
}
printf("agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%lf\n",ageminpar,agemaxpar, bage, fage, estepm, ftolpl);
ungetc(c,ficpar);
fscanf(ficpar,"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);
- fscanf(ficparo,"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);
- fscanf(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);
- fscanf(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);
+ fprintf(ficparo,"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(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.*/
hPijx(p, bage, fage);
fclose(ficrespij);
- ncovcombmax= pow(2,cptcoveff);
- /*-------------- Variance of one-step probabilities---*/
+ ncovcombmax= pow(2,cptcoveff);
+ /*-------------- Variance of one-step probabilities---*/
k=1;
varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart);
- /* Prevalence for each covariates in probs[age][status][cov] */
+ /* Prevalence for each covariates in probs[age][status][cov] */
probs= ma3x(1,AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
for(i=1;i<=AGESUP;i++)
- for(j=1;j<=nlstate;j++)
+ for(j=1;j<=nlstate+ndeath;j++) /* ndeath is useless but a necessity to be compared with mobaverages */
for(k=1;k<=ncovcombmax;k++)
probs[i][j][k]=0.;
- prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
- if (mobilav!=0 ||mobilavproj !=0 ) {
- mobaverage= ma3x(1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
- if (mobilav!=0) {
+ prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
+ if (mobilav!=0 ||mobilavproj !=0 ) {
+ mobaverages= ma3x(1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
+ for(i=1;i<=AGESUP;i++)
+ for(j=1;j<=nlstate;j++)
+ for(k=1;k<=ncovcombmax;k++)
+ mobaverages[i][j][k]=0.;
+ mobaverage=mobaverages;
+ if (mobilav!=0) {
if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilav)!=0){
fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
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 !=0) {
+ }
+ /* /\* 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) {
if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilavproj)!=0){
fprintf(ficlog," Error in movingaverage mobilavproj=%d\n",mobilavproj);
printf(" Error in movingaverage mobilavproj=%d\n",mobilavproj);
}
- }
- }/* end if moving average */
-
+ }
+ }/* end if moving average */
+
/*---------- Forecasting ------------------*/
/*if((stepm == 1) && (strcmp(model,".")==0)){*/
if(prevfcast==1){
prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
}
if(backcast==1){
- ddnewms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);
- ddoldms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);
- ddsavms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);
-
- /*--------------- Back Prevalence limit (period or stable prevalence) --------------*/
- /*#include "prevlim.h"*/ /* Use ficresplb, ficlog */
- bprlim=matrix(1,nlstate,1,nlstate);
- back_prevalence_limit(p, bprlim, ageminpar, agemaxpar, ftolpl, &ncvyear, dateprev1, dateprev2, firstpass, lastpass, mobilavproj);
- fclose(ficresplb);
-
- hBijx(p, bage, fage, mobaverage);
- 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); */
- 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);
- }
- /* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
- /* if (mobilavproj!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
-
- /* (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);*/
- /* } */
- /* else{ */
- /* erreur=108; */
- /* printf("Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */
- /* fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */
- /* } */
+ ddnewms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);
+ ddoldms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);
+ ddsavms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);
+
+ /*--------------- Back Prevalence limit (period or stable prevalence) --------------*/
+
+ bprlim=matrix(1,nlstate,1,nlstate);
+ back_prevalence_limit(p, bprlim, ageminpar, agemaxpar, ftolpl, &ncvyear, dateprev1, dateprev2, firstpass, lastpass, mobilavproj);
+ fclose(ficresplb);
+
+ hBijx(p, bage, fage, mobaverage);
+ 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); */
+ 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);
+ }
/* ------ Other prevalence ratios------------ */
- /* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */
-
- /* prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */
- /* printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d, mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\
- ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass);
- */
free_ivector(wav,1,imx);
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);
- /* if (mobilav!=0) { */
- /* mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
- /* if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){ */
- /* fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); */
- /* printf(" Error in movingaverage mobilav=%d\n",mobilav); */
- /* } */
- /* } */
-
-
/*---------- Health expectancies, no variances ------------*/
strcpy(filerese,"E_");
}
printf("Computing Health Expectancies: result on file '%s' ...", filerese);fflush(stdout);
fprintf(ficlog,"Computing Health Expectancies: result on file '%s' ...", filerese);fflush(ficlog);
- /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
- for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
-
+
for (k=1; k <= (int) pow(2,cptcoveff); k++){
- fprintf(ficreseij,"\n#****** ");
- for(j=1;j<=cptcoveff;j++) {
- fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- }
- fprintf(ficreseij,"******\n");
-
- eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
- oldm=oldms;savm=savms;
- evsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart);
+ fprintf(ficreseij,"\n#****** ");
+ for(j=1;j<=cptcoveff;j++) {
+ fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ }
+ fprintf(ficreseij,"******\n");
- free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
- /*}*/
+ eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
+ oldm=oldms;savm=savms;
+ evsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart);
+
+ free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
}
fclose(ficreseij);
printf("done evsij\n");fflush(stdout);
for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
- oldm=oldms;savm=savms; /* ZZ Segmentation fault */
- cptcod= 0; /* To be deleted */
- printf("varevsij %d \n",vpopbased);
- fprintf(ficlog, "varevsij %d \n",vpopbased);
- varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */
- fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are ");
- 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,"# 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"); */
- epj=vector(1,nlstate+1);
- 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");
- for(age=bage; age <=fage ;age++){
- prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k); /*ZZ Is it the correct prevalim */
- if (vpopbased==1) {
- if(mobilav ==0){
- for(i=1; i<=nlstate;i++)
- prlim[i][i]=probs[(int)age][i][k];
- }else{ /* mobilav */
- for(i=1; i<=nlstate;i++)
- prlim[i][i]=mobaverage[(int)age][i][k];
- }
- }
-
- fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav);
- /* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */
- /* printf(" age %4.0f ",age); */
- for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
- for(i=1, epj[j]=0.;i <=nlstate;i++) {
- epj[j] += prlim[i][i]*eij[i][j][(int)age];
- /*ZZZ printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
- /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */
- }
- epj[nlstate+1] +=epj[j];
- }
- /* printf(" age %4.0f \n",age); */
-
- for(i=1, vepp=0.;i <=nlstate;i++)
- for(j=1;j <=nlstate;j++)
- vepp += vareij[i][j][(int)age];
- fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
- for(j=1;j <=nlstate;j++){
- fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
- }
- fprintf(ficrest,"\n");
- }
+ oldm=oldms;savm=savms; /* ZZ Segmentation fault */
+ cptcod= 0; /* To be deleted */
+ printf("varevsij %d \n",vpopbased);
+ fprintf(ficlog, "varevsij %d \n",vpopbased);
+ varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */
+ fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are ");
+ 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,"# 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"); */
+ epj=vector(1,nlstate+1);
+ 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");
+ for(age=bage; age <=fage ;age++){
+ prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k); /*ZZ Is it the correct prevalim */
+ if (vpopbased==1) {
+ if(mobilav ==0){
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=probs[(int)age][i][k];
+ }else{ /* mobilav */
+ for(i=1; i<=nlstate;i++)
+ prlim[i][i]=mobaverage[(int)age][i][k];
+ }
+ }
+
+ fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav);
+ /* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */
+ /* printf(" age %4.0f ",age); */
+ for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
+ for(i=1, epj[j]=0.;i <=nlstate;i++) {
+ epj[j] += prlim[i][i]*eij[i][j][(int)age];
+ /*ZZZ printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
+ /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */
+ }
+ epj[nlstate+1] +=epj[j];
+ }
+ /* printf(" age %4.0f \n",age); */
+
+ for(i=1, vepp=0.;i <=nlstate;i++)
+ for(j=1;j <=nlstate;j++)
+ vepp += vareij[i][j][(int)age];
+ fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
+ for(j=1;j <=nlstate;j++){
+ fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
+ }
+ fprintf(ficrest,"\n");
+ }
} /* End vpopbased */
free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog);
/*---------- End : free ----------------*/
- if (mobilav!=0 ||mobilavproj !=0) free_ma3x(mobaverage,1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); /* We need to have a squared matrix with prevalence of the dead! */
+ if (mobilav!=0 ||mobilavproj !=0)
+ free_ma3x(mobaverages,1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); /* We need to have a squared matrix with prevalence of the dead! */
free_ma3x(probs,1,AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
} /* mle==-3 arrives here for freeing */
/* endfree:*/
git://henry.ined.fr / .git / commitdiff