|
version 1.218, 2016/02/12 11:29:23
|
version 1.219, 2016/02/15 00:48:12
|
|
Line 1
|
Line 1
|
| /* $Id$ |
/* $Id$ |
| $State$ |
$State$ |
| $Log$ |
$Log$ |
| |
Revision 1.219 2016/02/15 00:48:12 brouard |
| |
*** empty log message *** |
| |
|
| Revision 1.218 2016/02/12 11:29:23 brouard |
Revision 1.218 2016/02/12 11:29:23 brouard |
| Summary: 0.99 Back projections |
Summary: 0.99 Back projections |
| |
|
|
Line 976 int *ncodemaxwundef; /* ncodemax[j]= Nu
|
Line 979 int *ncodemaxwundef; /* ncodemax[j]= Nu
|
| undefined. Usually 3: -1, 0 and 1. */ |
undefined. Usually 3: -1, 0 and 1. */ |
| double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint; |
double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint; |
| double **pmmij, ***probs; /* Global pointer */ |
double **pmmij, ***probs; /* Global pointer */ |
| double ***mobaverage; /* New global variable */ |
double ***mobaverage, ***mobaverages; /* New global variable */ |
| double *ageexmed,*agecens; |
double *ageexmed,*agecens; |
| double dateintmean=0; |
double dateintmean=0; |
| |
|
|
Line 3969 void prevalence(double ***probs, double
|
Line 3972 void prevalence(double ***probs, double
|
| if (cptcovn<1) {j=1;ncodemax[1]=1;} |
if (cptcovn<1) {j=1;ncodemax[1]=1;} |
| |
|
| first=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 (i=1; i<=nlstate; i++) |
| for(iage=iagemin-AGEMARGE; iage <= iagemax+3+AGEMARGE; iage++) |
for(iage=iagemin-AGEMARGE; iage <= iagemax+3+AGEMARGE; iage++) |
| prop[i][iage]=0.0; |
prop[i][iage]=0.0; |
|
Line 3977 void prevalence(double ***probs, double
|
Line 3980 void prevalence(double ***probs, double
|
| for (i=1; i<=imx; i++) { /* Each individual */ |
for (i=1; i<=imx; i++) { /* Each individual */ |
| bool=1; |
bool=1; |
| if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ |
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)]) |
if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) |
| bool=0; |
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(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */ |
| for(mi=1; mi<wav[i];mi++){ |
for(mi=1; mi<wav[i];mi++){ |
| m=mw[mi][i]; |
m=mw[mi][i]; |
|
Line 4242 void tricode(int *Tvar, int **nbcode, in
|
Line 4245 void tricode(int *Tvar, int **nbcode, in
|
| for (j=1; j<=(cptcovs); j++) { /* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only */ |
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 (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 |
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 |
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 product of Vn*Vm, still boolean *: |
| * If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables |
* 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 */ |
* 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 |
/* Finds for covariate j, n=Tvar[j] of Vn . ij is the |
| modality of the nth covariate of individual i. */ |
modality of the nth covariate of individual i. */ |
| if (ij > modmaxcovj) |
if (ij > modmaxcovj) |
| modmaxcovj=ij; |
modmaxcovj=ij; |
| else if (ij < modmincovj) |
else if (ij < modmincovj) |
| modmincovj=ij; |
modmincovj=ij; |
| if ((ij < -1) && (ij > NCOVMAX)){ |
if ((ij < -1) && (ij > NCOVMAX)){ |
| printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX ); |
printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX ); |
| exit(1); |
exit(1); |
| }else |
}else |
| Ndum[ij]++; /*counts and stores the occurence of this modality 0, 1, -1*/ |
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 */ |
/* If coded 1, 2, 3 , counts the number of 1 Ndum[1], number of 2, Ndum[2], etc */ |
|
Line 4273 void tricode(int *Tvar, int **nbcode, in
|
Line 4276 void tricode(int *Tvar, int **nbcode, in
|
| printf("Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], k, Ndum[k]); |
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]); |
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( Ndum[k] != 0 ){ /* Counts if nobody answered modality k ie empty modality, we skip it and reorder */ |
| if( k != -1){ |
if( k != -1){ |
| ncodemax[j]++; /* ncodemax[j]= Number of modalities of the j th |
ncodemax[j]++; /* ncodemax[j]= Number of modalities of the j th |
| covariate for which somebody answered excluding |
covariate for which somebody answered excluding |
| undefined. Usually 2: 0 and 1. */ |
undefined. Usually 2: 0 and 1. */ |
| } |
} |
| ncodemaxwundef[j]++; /* ncodemax[j]= Number of modalities of the j th |
ncodemaxwundef[j]++; /* ncodemax[j]= Number of modalities of the j th |
| covariate for which somebody answered including |
covariate for which somebody answered including |
| undefined. Usually 3: -1, 0 and 1. */ |
undefined. Usually 3: -1, 0 and 1. */ |
| } |
} |
| /* In fact ncodemax[j]=2 (dichotom. variables only) but it could be more for |
/* 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 */ |
} /* Ndum[-1] number of undefined modalities */ |
| |
|
| /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */ |
/* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */ |
| /* For covariate j, modalities could be 1, 2, 3, 4, 5, 6, 7. |
/* 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; |
If Ndum[1]=0, Ndum[2]=0, Ndum[3]= 635, Ndum[4]=0, Ndum[5]=0, Ndum[6]=27, Ndum[7]=125; |
|
Line 4302 void tricode(int *Tvar, int **nbcode, in
|
Line 4305 void tricode(int *Tvar, int **nbcode, in
|
| ij=0; /* ij is similar to i but can jump over null modalities */ |
ij=0; /* ij is similar to i but can jump over null modalities */ |
| for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/ |
for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/ |
| if (Ndum[i] == 0) { /* If nobody responded to this modality k */ |
if (Ndum[i] == 0) { /* If nobody responded to this modality k */ |
| break; |
break; |
| } |
} |
| ij++; |
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.*/ |
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 */ |
cptcode = ij; /* New max modality for covar j */ |
|
Line 4324 void tricode(int *Tvar, int **nbcode, in
|
Line 4327 void tricode(int *Tvar, int **nbcode, in
|
| /* } /\* end of loop on modality k *\/ */ |
/* } /\* end of loop on modality k *\/ */ |
| } /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/ |
} /* end of loop on model-covariate 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 */ |
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.*/ |
/* 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 */ |
ij=Tvar[i]; /* Tvar might be -1 if status was unknown */ |
| Ndum[ij]++; /* Might be supersed V1 + V1*age */ |
Ndum[ij]++; /* Might be supersed V1 + V1*age */ |
| } |
} |
| |
|
| ij=0; |
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) */ |
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]);*/ |
/*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/ |
| if((Ndum[i]!=0) && (i<=ncovcol)){ |
if((Ndum[i]!=0) && (i<=ncovcol)){ |
| ij++; |
ij++; |
| /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/ |
/*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/ |
| Tvaraff[ij]=i; /*For printing (unclear) */ |
Tvaraff[ij]=i; /*For printing (unclear) */ |
| }else{ |
}else{ |
| /* Tvaraff[ij]=0; */ |
/* Tvaraff[ij]=0; */ |
| } |
} |
| } |
} |
| /* ij--; */ |
/* ij--; */ |
| cptcoveff=ij; /*Number of total covariates*/ |
cptcoveff=ij; /*Number of total covariates*/ |
| |
|
| } |
} |
| |
|
| |
|
|
Line 5632 true period expectancies (those weighted
|
Line 5635 true period expectancies (those weighted
|
| int lv=0, vlv=0, kl=0; |
int lv=0, vlv=0, kl=0; |
| int ng=0; |
int ng=0; |
| int vpopbased; |
int vpopbased; |
| |
int ioffset; /* variable offset for columns */ |
| |
|
| /* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */ |
/* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */ |
| /* printf("Problem with file %s",optionfilegnuplot); */ |
/* printf("Problem with file %s",optionfilegnuplot); */ |
| /* fprintf(ficlog,"Problem with file %s",optionfilegnuplot); */ |
/* fprintf(ficlog,"Problem with file %s",optionfilegnuplot); */ |
|
Line 5661 true period expectancies (those weighted
|
Line 5666 true period expectancies (those weighted
|
| fprintf(ficgp,"unset log;\n# plot weighted, mean weight should have point size of 0.5\n plot \"%s\"",subdirf(fileresilk)); |
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); |
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 ++) { |
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"); |
fprintf(ficgp,";\nset out; unset ylabel;\n"); |
| } |
} |
|
Line 5678 true period expectancies (those weighted
|
Line 5683 true period expectancies (those weighted
|
| for (k1=1; k1<= m ; k1 ++) { /* For each combination of covariate */ |
for (k1=1; k1<= m ; k1 ++) { /* For each combination of covariate */ |
| /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */ |
/* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */ |
| fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files "); |
fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files "); |
| for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
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 covariate number corresponding to k1 combination */ |
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,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(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 */ |
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| vlv= nbcode[Tvaraff[lv]][lv]; |
vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */ |
| fprintf(ficgp," V%d=%d ",k,vlv); |
/* 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,"\n#\n"); |
| |
|
| fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1); |
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,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1); |
| fprintf(ficgp,"set xlabel \"Age\" \n\ |
fprintf(ficgp,"set xlabel \"Age\" \n\ |
| set ylabel \"Probability\" \n\ |
set ylabel \"Probability\" \n \ |
| set ter svg size 640, 480\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); |
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 ++) { |
for (i=1; i<= nlstate ; i ++) { |
| if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
| else 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); |
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 ++) { |
for (i=1; i<= nlstate ; i ++) { |
| if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
| else 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); |
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 ++) { |
for (i=1; i<= nlstate ; i ++) { |
| if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
| else 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)); |
| if(backcast==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\" 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 */ |
| fprintf(ficgp,"\nset out \n"); |
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 */ |
} /* k1 */ |
| } /* cpt */ |
} /* cpt */ |
| /*2 eme*/ |
/*2 eme*/ |
| for (k1=1; k1<= m ; k1 ++) { |
for (k1=1; k1<= m ; k1 ++) { |
| fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files "); |
fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files "); |
| for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */ |
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 */ |
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,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(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 */ |
/* 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," V%d=%d ",k,vlv); |
| } |
} |
| fprintf(ficgp,"\n#\n"); |
fprintf(ficgp,"\n#\n"); |
| |
|
| fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1); |
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*/ |
for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ |
| if(vpopbased==0) |
if(vpopbased==0) |
| fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage); |
fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage); |
| else |
else |
| fprintf(ficgp,"\nreplot "); |
fprintf(ficgp,"\nreplot "); |
| for (i=1; i<= nlstate+1 ; i ++) { |
for (i=1; i<= nlstate+1 ; i ++) { |
| k=2*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); |
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 ++) { |
for (j=1; j<= nlstate+1 ; j ++) { |
| if (j==i) fprintf(ficgp," %%lf (%%lf)"); |
if (j==i) fprintf(ficgp," %%lf (%%lf)"); |
| else fprintf(ficgp," %%*lf (%%*lf)"); |
else fprintf(ficgp," %%*lf (%%*lf)"); |
| } |
} |
| if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i); |
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); |
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); |
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 ++) { |
for (j=1; j<= nlstate+1 ; j ++) { |
| if (j==i) fprintf(ficgp," %%lf (%%lf)"); |
if (j==i) fprintf(ficgp," %%lf (%%lf)"); |
| else fprintf(ficgp," %%*lf (%%*lf)"); |
else fprintf(ficgp," %%*lf (%%*lf)"); |
| } |
} |
| fprintf(ficgp,"\" t\"\" w l lt 0,"); |
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); |
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 ++) { |
for (j=1; j<= nlstate+1 ; j ++) { |
| if (j==i) fprintf(ficgp," %%lf (%%lf)"); |
if (j==i) fprintf(ficgp," %%lf (%%lf)"); |
| else fprintf(ficgp," %%*lf (%%*lf)"); |
else fprintf(ficgp," %%*lf (%%*lf)"); |
| } |
} |
| if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0"); |
if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0"); |
| else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n"); |
else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n"); |
| } /* state */ |
} /* state */ |
| } /* vpopbased */ |
} /* vpopbased */ |
| fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */ |
fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */ |
| } /* k1 */ |
} /* k1 */ |
| |
|
| |
|
| /*3eme*/ |
/*3eme*/ |
| for (k1=1; k1<= m ; k1 ++) { |
for (k1=1; k1<= m ; k1 ++) { |
| for (cpt=1; cpt<= nlstate ; cpt ++) { |
for (cpt=1; cpt<= nlstate ; cpt ++) { |
| fprintf(ficgp,"\n# 3d: Life expectancy with EXP_ files: cov=%d state=%d",k1, 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 */ |
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 */ |
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,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(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 */ |
/* 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," V%d=%d ",k,vlv); |
| } |
} |
| fprintf(ficgp,"\n#\n"); |
fprintf(ficgp,"\n#\n"); |
| |
|
| /* k=2+nlstate*(2*cpt-2); */ |
/* k=2+nlstate*(2*cpt-2); */ |
| k=2+(nlstate+1)*(cpt-1); |
k=2+(nlstate+1)*(cpt-1); |
| fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1); |
fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1); |
| fprintf(ficgp,"set ter svg size 640, 480\n\ |
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); |
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); |
/*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) "); |
for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) "); |
| fprintf(ficgp,"\" t \"e%d1\" w l",cpt); |
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); |
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) "); |
for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) "); |
| fprintf(ficgp,"\" t \"e%d1\" w l",cpt); |
fprintf(ficgp,"\" t \"e%d1\" w l",cpt); |
| |
|
| */ |
*/ |
| for (i=1; i< nlstate ; i ++) { |
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(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(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); |
fprintf(ficgp," ,\"%s\" every :::%d::%d u 1:%d t \"e%d.\" w l",subdirf2(fileresu,"E_"),k1-1,k1-1,k+nlstate,cpt); |
| } |
} |
|
Line 5808 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 5835 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
| /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
/* 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(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 */ |
/* 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," V%d=%d ",k,vlv); |
| } |
} |
| fprintf(ficgp,"\n#\n"); |
fprintf(ficgp,"\n#\n"); |
|
Line 5844 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 5871 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
/* 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(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 */ |
/* 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," V%d=%d ",k,vlv); |
| } |
} |
| fprintf(ficgp,"\n#\n"); |
fprintf(ficgp,"\n#\n"); |
|
Line 5888 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 5915 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
/* 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(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 */ |
/* 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," V%d=%d ",k,vlv); |
| } |
} |
| fprintf(ficgp,"\n#\n"); |
fprintf(ficgp,"\n#\n"); |
|
Line 5923 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 5950 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
/* 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(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 */ |
/* 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," V%d=%d ",k,vlv); |
| } |
} |
| fprintf(ficgp,"\n#\n"); |
fprintf(ficgp,"\n#\n"); |
|
Line 5959 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 5986 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| |
|
| for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */ |
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 */ |
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); |
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 */ |
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 */ |
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,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(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 */ |
/* 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," V%d=%d ",k,vlv); |
| } |
} |
| fprintf(ficgp,"\n#\n"); |
fprintf(ficgp,"\n#\n"); |
| |
|
| fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\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,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1); |
| fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\ |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\ |
| set ter svg size 640, 480\n\ |
set ter svg size 640, 480\n \ |
| unset log y\n\ |
unset log y\n \ |
| plot [%.f:%.f] ", ageminpar, agemaxpar); |
plot [%.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 */ |
| /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ |
/*# 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 */ |
/*# 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*/ |
/*# 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 */ |
/*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ |
| if(i==1){ |
if(i==1){ |
| fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"F_")); |
fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"F_")); |
| }else{ |
}else{ |
| fprintf(ficgp,",\\\n '' "); |
fprintf(ficgp,",\\\n '' "); |
| } |
} |
| if(cptcoveff ==0){ /* No covariate */ |
if(cptcoveff ==0){ /* No covariate */ |
| fprintf(ficgp," u 2:("); /* Age is in 2 */ |
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*/ |
/*# 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 */ |
/*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */ |
| if(i==nlstate+1) |
/*# V1 = 1 yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/ |
| fprintf(ficgp," $%d/(1.-$%d)) t 'p.%d' with line ", \ |
/*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */ |
| 2+(cpt-1)*(nlstate+1)+1+(i-1), 2+1+(i-1)+(nlstate+1)*nlstate,cpt ); |
fprintf(ficgp," u %d:(", ioffset); |
| else |
if(i==nlstate+1) |
| fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \ |
fprintf(ficgp," $%d/(1.-$%d)) t 'pw.%d' with line ", \ |
| 2+(cpt-1)*(nlstate+1)+1+(i-1), 2+1+(i-1)+(nlstate+1)*nlstate,i,cpt ); |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt ); |
| }else{ |
else |
| fprintf(ficgp,"u 6:(("); /* Age is in 6 */ |
fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \ |
| /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,i,cpt ); |
| /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ |
}else{ /* more than 2 covariates */ |
| kl=0; |
if(cptcoveff ==1){ |
| for (k=1; k<=cptcoveff; k++){ /* For each covariate */ |
ioffset=4; /* Age is in 4 */ |
| lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */ |
}else{ |
| /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
ioffset=6; /* Age is in 6 */ |
| /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
/*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ |
| /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
/*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ |
| vlv= nbcode[Tvaraff[lv]][lv]; |
} |
| kl++; |
fprintf(ficgp," u %d:((",ioffset); |
| /* 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 *\/ */ |
kl=0; |
| /*6+(cpt-1)*(nlstate+1)+1+(i-1)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */ |
for (k=1; k<=cptcoveff; k++){ /* For each covariate */ |
| /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */ |
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to combination k1 and covariate k */ |
| /* '' 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*/ |
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ |
| if(k==cptcoveff) |
/* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ |
| if(i==nlstate+1) |
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ |
| fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv], \ |
vlv= nbcode[Tvaraff[k]][lv]; |
| 6+(cpt-1)*(nlstate+1)+1+(i-1), 6+1+(i-1)+(nlstate+1)*nlstate,cpt ); |
kl++; |
| else |
/* 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 *\/ */ |
| 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)+(nlstate+1)*nlstate; 6+(1-1)*(2+1)+1+(1-1) +(2+1)*2=13 */ |
| 6+(cpt-1)*(nlstate+1)+1+(i-1), 6+1+(i-1)+(nlstate+1)*nlstate,i,cpt ); |
/*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */ |
| else{ |
/* '' 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*/ |
| fprintf(ficgp,"$%d==%d && $%d==%d && ",kl, k,kl+1,nbcode[Tvaraff[lv]][lv]); |
if(k==cptcoveff){ |
| kl++; |
if(i==nlstate+1){ |
| } |
if(cptcoveff ==1){ |
| } /* end covariate */ |
fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ",kl, k,kl+1,nbcode[Tvaraff[k]][lv], \ |
| } /* end if covariate */ |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt ); |
| } /* nlstate */ |
}else{ |
| fprintf(ficgp,"\nset out\n"); |
fprintf(ficgp,"$%d==%d && $%d==%d)? $%d/(1.-$%d) : 1/0) t 'p.%d' with line ",kl, k,kl+1,nbcode[Tvaraff[k]][lv], \ |
| } /* end cpt state*/ |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt ); |
| } /* end covariate */ |
} |
| } /* End if prevfcast */ |
}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], \ |
| /* proba elementaires */ |
ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset +1+(i-1)+(nlstate+1)*nlstate,i,cpt ); |
| fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n"); |
}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++){ |
for(i=1,jk=1; i <=nlstate; i++){ |
| fprintf(ficgp,"# initial state %d\n",i); |
fprintf(ficgp,"# initial state %d\n",i); |
| for(k=1; k <=(nlstate+ndeath); k++){ |
for(k=1; k <=(nlstate+ndeath); k++){ |
|
Line 6175 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 6221 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| |
|
| |
|
| /*************** Moving average **************/ |
/*************** 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 movingaverage(double ***probs, double bage, double fage, double ***mobaverage, int mobilav){ |
| |
|
| int i, cpt, cptcod; |
int i, cpt, cptcod; |
| int modcovmax =1; |
int modcovmax =1; |
| int mobilavrange, mob; |
int mobilavrange, mob; |
| double age; |
|
| int iage=0; |
int iage=0; |
| |
|
| |
double sum=0.; |
| |
double age; |
| double *sumnewp, *sumnewm; |
double *sumnewp, *sumnewm; |
| double *agemingood, *agemaxgood; /* Currently identical for all covariates */ |
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); |
sumnewp = vector(1,modcovmax); |
| sumnewm = vector(1,modcovmax); |
sumnewm = vector(1,modcovmax); |
| agemingood = vector(1,modcovmax); |
agemingood = vector(1,modcovmax); |
| agemaxgood = vector(1,modcovmax); |
agemaxgood = vector(1,modcovmax); |
| |
|
| |
for (cptcod=1;cptcod<=modcovmax;cptcod++){ |
| modcovmax=2*cptcoveff;/* Max number of modalities. We suppose |
sumnewm[cptcod]=0.; |
| a covariate has 2 modalities, should be equal to ncovcombmax */ |
sumnewp[cptcod]=0.; |
| |
agemingood[cptcod]=0; |
| |
agemaxgood[cptcod]=0; |
| |
} |
| if (cptcovn<1) modcovmax=1; /* At least 1 pass */ |
if (cptcovn<1) modcovmax=1; /* At least 1 pass */ |
| |
|
| if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){ |
if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){ |
|
Line 6201 int movingaverage(double ***probs, doubl
|
Line 6256 int movingaverage(double ***probs, doubl
|
| else mobilavrange=mobilav; |
else mobilavrange=mobilav; |
| for (age=bage; age<=fage; age++) |
for (age=bage; age<=fage; age++) |
| for (i=1; i<=nlstate;i++) |
for (i=1; i<=nlstate;i++) |
| for (cptcod=1;cptcod<=modcovmax;cptcod++) |
for (cptcod=1;cptcod<=modcovmax;cptcod++) |
| mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod]; |
mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod]; |
| /* We keep the original values on the extreme ages bage, fage and for |
/* 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 |
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. |
we use a 5 terms etc. until the borders are no more concerned. |
| */ |
*/ |
| for (mob=3;mob <=mobilavrange;mob=mob+2){ |
for (mob=3;mob <=mobilavrange;mob=mob+2){ |
| for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ |
for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ |
| for (i=1; i<=nlstate;i++){ |
for (i=1; i<=nlstate;i++){ |
| for (cptcod=1;cptcod<=modcovmax;cptcod++){ |
for (cptcod=1;cptcod<=modcovmax;cptcod++){ |
| mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod]; |
mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod]; |
| for (cpt=1;cpt<=(mob-1)/2;cpt++){ |
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] +=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; |
mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/mob; |
| } |
} |
| } |
} |
| }/* end age */ |
}/* end age */ |
| }/* end mob */ |
}/* end mob */ |
| }else |
}else |
| return -1; |
return -1; |
| for (cptcod=1;cptcod<=modcovmax;cptcod++){ |
for (cptcod=1;cptcod<=modcovmax;cptcod++){ |
| /* for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ */ |
/* 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 */ |
for (age=fage-(mob-1)/2; age>=bage; age--){/* From oldest to youngest, finding the youngest wrong */ |
| sumnewm[cptcod]=0.; |
sumnewm[cptcod]=0.; |
| for (i=1; i<=nlstate;i++){ |
for (i=1; i<=nlstate;i++){ |
| sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod]; |
sumnewm[cptcod]+=mobaverage[(int)age][i][cptcod]; |
| } |
} |
| if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */ |
if(fabs(sumnewm[cptcod] - 1.) <= 1.e-3) { /* good */ |
| agemingood[cptcod]=age; |
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; |
|
| }else{ /* bad */ |
}else{ /* bad */ |
| for (i=1; i<=nlstate;i++){ |
for (i=1; i<=nlstate;i++){ |
| mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][i][cptcod]; |
mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][i][cptcod]; |
| } /* i */ |
} /* i */ |
| } /* end bad */ |
} /* end bad */ |
| }/* age */ |
}/* age */ |
| for (age=bage; age<=fage; age++){ |
sum=0.; |
| printf("%d %d ", cptcod, (int)age); |
for (i=1; i<=nlstate;i++){ |
| sumnewp[cptcod]=0.; |
sum+=mobaverage[(int)agemingood[cptcod]][i][cptcod]; |
| 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"); |
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 */ |
/* brutal averaging */ |
| for (i=1; i<=nlstate;i++){ |
for (i=1; i<=nlstate;i++){ |
| for (age=1; age<=bage; age++){ |
for (age=1; age<=bage; age++){ |
| mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemingood[cptcod]][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]); |
/* 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++){ |
for (age=fage; age<=AGESUP; age++){ |
| mobaverage[(int)age][i][cptcod]=mobaverage[(int)agemaxgood[cptcod]][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]); |
/* printf("age=%d i=%d cptcod=%d mobaverage=%.4f \n",(int)age,i, cptcod, mobaverage[(int)age][i][cptcod]); */ |
| } |
} |
| } /* end i status */ |
} /* end i status */ |
| for (i=nlstate+1; i<=nlstate+ndeath;i++){ |
for (i=nlstate+1; i<=nlstate+ndeath;i++){ |
| for (age=1; age<=AGESUP; age++){ |
for (age=1; age<=AGESUP; age++){ |
| /*printf("i=%d, age=%d, cptcod=%d\n",i, (int)age, cptcod);*/ |
/*printf("i=%d, age=%d, cptcod=%d\n",i, (int)age, cptcod);*/ |
| mobaverage[(int)age][i][cptcod]=0.; |
mobaverage[(int)age][i][cptcod]=0.; |
| } |
} |
| } |
} |
| }/* end cptcod */ |
}/* end cptcod */ |
|
Line 6358 void prevforecast(char fileres[], double
|
Line 6436 void prevforecast(char fileres[], double
|
| k=k+1; |
k=k+1; |
| fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#"); |
fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#"); |
| for(j=1;j<=cptcoveff;j++) { |
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"); |
fprintf(ficresf," yearproj age"); |
| for(j=1; j<=nlstate+ndeath;j++){ |
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%d",i,j); |
| fprintf(ficresf," p.%d",j); |
fprintf(ficresf," wp.%d",j); |
| } |
} |
| for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { |
for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { |
| fprintf(ficresf,"\n"); |
fprintf(ficresf,"\n"); |
| fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp); |
fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp); |
| for (agec=fage; agec>=(ageminpar-1); agec--){ |
for (agec=fage; agec>=(ageminpar-1); agec--){ |
| nhstepm=(int) rint((agelim-agec)*YEARM/stepm); |
nhstepm=(int) rint((agelim-agec)*YEARM/stepm); |
| nhstepm = nhstepm/hstepm; |
nhstepm = nhstepm/hstepm; |
| p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
| oldm=oldms;savm=savms; |
oldm=oldms;savm=savms; |
| hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k); |
hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k); |
| |
|
| for (h=0; h<=nhstepm; h++){ |
for (h=0; h<=nhstepm; h++){ |
| if (h*hstepm/YEARM*stepm ==yearp) { |
if (h*hstepm/YEARM*stepm ==yearp) { |
| fprintf(ficresf,"\n"); |
fprintf(ficresf,"\n"); |
| for(j=1;j<=cptcoveff;j++) |
for(j=1;j<=cptcoveff;j++) |
| fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
| fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm); |
fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm); |
| } |
} |
| for(j=1; j<=nlstate+ndeath;j++) { |
for(j=1; j<=nlstate+ndeath;j++) { |
| ppij=0.; |
ppij=0.; |
| for(i=1; i<=nlstate;i++) { |
for(i=1; i<=nlstate;i++) { |
| if (mobilav==1) |
if (mobilav==1) |
| ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][cptcod]; |
ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][cptcod]; |
| else { |
else { |
| ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod]; |
ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod]; |
| } |
} |
| if (h*hstepm/YEARM*stepm== yearp) { |
if (h*hstepm/YEARM*stepm== yearp) { |
| fprintf(ficresf," %.3f", p3mat[i][j][h]); |
fprintf(ficresf," %.3f", p3mat[i][j][h]); |
| } |
} |
| } /* end i */ |
} /* end i */ |
| if (h*hstepm/YEARM*stepm==yearp) { |
if (h*hstepm/YEARM*stepm==yearp) { |
| fprintf(ficresf," %.3f", ppij); |
fprintf(ficresf," %.3f", ppij); |
| } |
} |
| }/* end j */ |
}/* end j */ |
| } /* end h */ |
} /* end h */ |
| free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
| } /* end agec */ |
} /* end agec */ |
| } /* end yearp */ |
} /* end yearp */ |
| } /* end cptcod */ |
} /* end cptcod */ |
| } /* end cptcov */ |
} /* end cptcov */ |
| |
|
| fclose(ficresf); |
fclose(ficresf); |
| printf("End of Computing forecasting \n"); |
printf("End of Computing forecasting \n"); |
| fprintf(ficlog,"End of Computing forecasting\n"); |
fprintf(ficlog,"End of Computing forecasting\n"); |
|
Line 7656 void syscompilerinfo(int logged)
|
Line 7734 void syscompilerinfo(int logged)
|
| #endif |
#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) --------------*/ |
/*--------------- Prevalence limit (period or stable prevalence) --------------*/ |
| int i, j, k, i1 ; |
int i, j, k, i1 ; |
| /* double ftolpl = 1.e-10; */ |
/* double ftolpl = 1.e-10; */ |
|
Line 7679 void syscompilerinfo(int logged)
|
Line 7757 void syscompilerinfo(int logged)
|
| for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); |
for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); |
| fprintf(ficrespl,"\n"); |
fprintf(ficrespl,"\n"); |
| |
|
| /* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */ |
/* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */ |
| |
|
| agebase=ageminpar; |
agebase=ageminpar; |
| agelim=agemaxpar; |
agelim=agemaxpar; |
| |
|
| i1=pow(2,cptcoveff); |
i1=pow(2,cptcoveff); |
| if (cptcovn < 1){i1=1;} |
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(cptcov=1,k=0;cptcov<=1;cptcov++){ */ |
| //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
//for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
| k=k+1; |
k=k+1; |
| /* to clean */ |
/* to clean */ |
| //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov)); |
//printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov)); |
| fprintf(ficrespl,"#******"); |
fprintf(ficrespl,"#******"); |
| printf("#******"); |
printf("#******"); |
| fprintf(ficlog,"#******"); |
fprintf(ficlog,"#******"); |
| for(j=1;j<=cptcoveff;j++) { |
for(j=1;j<=cptcoveff;j++) { |
| fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,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)]); |
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(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
| } |
} |
| fprintf(ficrespl,"******\n"); |
fprintf(ficrespl,"******\n"); |
| printf("******\n"); |
printf("******\n"); |
| fprintf(ficlog,"******\n"); |
fprintf(ficlog,"******\n"); |
| |
|
| fprintf(ficrespl,"#Age "); |
fprintf(ficrespl,"#Age "); |
| for(j=1;j<=cptcoveff;j++) { |
for(j=1;j<=cptcoveff;j++) { |
| fprintf(ficrespl,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,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); |
for(i=1; i<=nlstate;i++) fprintf(ficrespl," %d-%d ",i,i); |
| fprintf(ficrespl,"Total Years_to_converge\n"); |
fprintf(ficrespl,"Total Years_to_converge\n"); |
| |
|
| for (age=agebase; age<=agelim; age++){ |
for (age=agebase; age<=agelim; age++){ |
| /* for (age=agebase; age<=agebase; age++){ */ |
/* for (age=agebase; age<=agebase; age++){ */ |
| prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k); |
prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k); |
| fprintf(ficrespl,"%.0f ",age ); |
fprintf(ficrespl,"%.0f ",age ); |
| for(j=1;j<=cptcoveff;j++) |
for(j=1;j<=cptcoveff;j++) |
| fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
| tot=0.; |
tot=0.; |
| for(i=1; i<=nlstate;i++){ |
for(i=1; i<=nlstate;i++){ |
| tot += prlim[i][i]; |
tot += prlim[i][i]; |
| fprintf(ficrespl," %.5f", prlim[i][i]); |
fprintf(ficrespl," %.5f", prlim[i][i]); |
| } |
} |
| fprintf(ficrespl," %.3f %d\n", tot, *ncvyearp); |
fprintf(ficrespl," %.3f %d\n", tot, *ncvyearp); |
| } /* Age */ |
} /* Age */ |
| /* was end of cptcod */ |
/* was end of cptcod */ |
| } /* cptcov */ |
} /* cptcov */ |
| return 0; |
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){ |
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){ |
|
Line 7798 int back_prevalence_limit(double *p, dou
|
Line 7876 int back_prevalence_limit(double *p, dou
|
| if(mobilavproj > 0){ |
if(mobilavproj > 0){ |
| /* bprevalim(bprlim, mobaverage, nlstate, p, age, ageminpar, agemaxpar, oldm, savm, doldm, dsavm, ftolpl, ncvyearp, k); */ |
/* 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, 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){ |
}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); |
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); |
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); |
exit(1); |
| }else{ |
}else{ |
| /* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, 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); |
bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k); |
| } |
} |
| fprintf(ficresplb,"%.0f ",age ); |
fprintf(ficresplb,"%.0f ",age ); |
| for(j=1;j<=cptcoveff;j++) |
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.; |
tot=0.; |
| for(i=1; i<=nlstate;i++){ |
for(i=1; i<=nlstate;i++){ |
| tot += bprlim[i][i]; |
tot += bprlim[i][i]; |
| fprintf(ficresplb," %.5f", bprlim[i][i]); |
fprintf(ficresplb," %.5f", bprlim[i][i]); |
| } |
} |
| fprintf(ficresplb," %.3f %d\n", tot, *ncvyearp); |
fprintf(ficresplb," %.3f %d\n", tot, *ncvyearp); |
| } /* Age */ |
} /* Age */ |
|
Line 8004 int main(int argc, char *argv[])
|
Line 8082 int main(int argc, char *argv[])
|
| double agedeb=0.; |
double agedeb=0.; |
| |
|
| double ageminpar=AGEOVERFLOW,agemin=AGEOVERFLOW, agemaxpar=-AGEOVERFLOW, agemax=-AGEOVERFLOW; |
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 fret; |
| double dum=0.; /* Dummy variable */ |
double dum=0.; /* Dummy variable */ |
|
Line 8702 Please run with mle=-1 to get a correct
|
Line 8780 Please run with mle=-1 to get a correct
|
| * bbbbbbbb |
* bbbbbbbb |
| * 76543210 |
* 76543210 |
| * h-1 00000101 (6-1=5) |
* 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) |
* 1 00000001 (1) |
| * 00000001 = 1 & ((h-1) >> (k-1)) |
* 00000000 = 1 & ((h-1) >> (k-1)) |
| * +1= 00000010 =2 |
* +1= 00000001 =1 |
| * |
* |
| * h=14, k=3 => h'=h-1=13, k'=k-1=2 |
* h=14, k=3 => h'=h-1=13, k'=k-1=2 |
| * h' 1101 =2^3+2^2+0x2^1+2^0 |
* h' 1101 =2^3+2^2+0x2^1+2^0 |
|
Line 9284 Please run with mle=-1 to get a correct
|
Line 9362 Please run with mle=-1 to get a correct
|
| 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=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) { |
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("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); |
| printf("but line=%s\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; |
goto end; |
| } |
} |
| printf("agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%lf\n",ageminpar,agemaxpar, bage, fage, estepm, ftolpl); |
printf("agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d ftolpl=%lf\n",ageminpar,agemaxpar, bage, fage, estepm, ftolpl); |
|
Line 9360 Please run with mle=-1 to get a correct
|
Line 9438 Please run with mle=-1 to get a correct
|
| ungetc(c,ficpar); |
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(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); |
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); |
| 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); |
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); |
| 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(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.*/ |
/* day and month of proj2 are not used but only year anproj2.*/ |
| |
|
| |
|
|
Line 9413 Please run with mle=-1 to get a correct
|
Line 9491 Please run with mle=-1 to get a correct
|
| hPijx(p, bage, fage); |
hPijx(p, bage, fage); |
| fclose(ficrespij); |
fclose(ficrespij); |
| |
|
| ncovcombmax= pow(2,cptcoveff); |
ncovcombmax= pow(2,cptcoveff); |
| /*-------------- Variance of one-step probabilities---*/ |
/*-------------- Variance of one-step probabilities---*/ |
| k=1; |
k=1; |
| varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart); |
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); |
probs= ma3x(1,AGESUP,1,nlstate+ndeath, 1,ncovcombmax); |
| for(i=1;i<=AGESUP;i++) |
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++) |
for(k=1;k<=ncovcombmax;k++) |
| probs[i][j][k]=0.; |
probs[i][j][k]=0.; |
| prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); |
prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); |
| if (mobilav!=0 ||mobilavproj !=0 ) { |
if (mobilav!=0 ||mobilavproj !=0 ) { |
| mobaverage= ma3x(1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); |
mobaverages= ma3x(1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); |
| if (mobilav!=0) { |
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){ |
if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilav)!=0){ |
| fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); |
fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); |
| printf(" 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 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); */ |
/* prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */ |
| else if (mobilavproj !=0) { |
else if (mobilavproj !=0) { |
| if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilavproj)!=0){ |
if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilavproj)!=0){ |
| fprintf(ficlog," Error in movingaverage mobilavproj=%d\n",mobilavproj); |
fprintf(ficlog," Error in movingaverage mobilavproj=%d\n",mobilavproj); |
| printf(" Error in movingaverage mobilavproj=%d\n",mobilavproj); |
printf(" Error in movingaverage mobilavproj=%d\n",mobilavproj); |
| } |
} |
| } |
} |
| }/* end if moving average */ |
}/* end if moving average */ |
| |
|
| /*---------- Forecasting ------------------*/ |
/*---------- Forecasting ------------------*/ |
| /*if((stepm == 1) && (strcmp(model,".")==0)){*/ |
/*if((stepm == 1) && (strcmp(model,".")==0)){*/ |
| if(prevfcast==1){ |
if(prevfcast==1){ |
|
Line 9450 Please run with mle=-1 to get a correct
|
Line 9533 Please run with mle=-1 to get a correct
|
| prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff); |
prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff); |
| } |
} |
| if(backcast==1){ |
if(backcast==1){ |
| ddnewms=matrix(1,nlstate+ndeath,1,nlstate+ndeath); |
ddnewms=matrix(1,nlstate+ndeath,1,nlstate+ndeath); |
| ddoldms=matrix(1,nlstate+ndeath,1,nlstate+ndeath); |
ddoldms=matrix(1,nlstate+ndeath,1,nlstate+ndeath); |
| ddsavms=matrix(1,nlstate+ndeath,1,nlstate+ndeath); |
ddsavms=matrix(1,nlstate+ndeath,1,nlstate+ndeath); |
| |
|
| /*--------------- Back Prevalence limit (period or stable prevalence) --------------*/ |
/*--------------- Back Prevalence limit (period or stable prevalence) --------------*/ |
| /*#include "prevlim.h"*/ /* Use ficresplb, ficlog */ |
|
| bprlim=matrix(1,nlstate,1,nlstate); |
bprlim=matrix(1,nlstate,1,nlstate); |
| back_prevalence_limit(p, bprlim, ageminpar, agemaxpar, ftolpl, &ncvyear, dateprev1, dateprev2, firstpass, lastpass, mobilavproj); |
back_prevalence_limit(p, bprlim, ageminpar, agemaxpar, ftolpl, &ncvyear, dateprev1, dateprev2, firstpass, lastpass, mobilavproj); |
| fclose(ficresplb); |
fclose(ficresplb); |
| |
|
| hBijx(p, bage, fage, mobaverage); |
hBijx(p, bage, fage, mobaverage); |
| fclose(ficrespijb); |
fclose(ficrespijb); |
| free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */ |
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, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj, |
| free_matrix(ddnewms, 1, nlstate+ndeath, 1, nlstate+ndeath); |
bage, fage, firstpass, lastpass, anback2, p, cptcoveff); */ |
| free_matrix(ddsavms, 1, nlstate+ndeath, 1, nlstate+ndeath); |
free_matrix(ddnewms, 1, nlstate+ndeath, 1, nlstate+ndeath); |
| free_matrix(ddoldms, 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); */ |
|
| /* } */ |
|
| |
|
| |
|
| /* ------ Other prevalence ratios------------ */ |
/* ------ 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_ivector(wav,1,imx); |
| free_imatrix(dh,1,lastpass-firstpass+2,1,imx); |
free_imatrix(dh,1,lastpass-firstpass+2,1,imx); |
| free_imatrix(bh,1,lastpass-firstpass+2,1,imx); |
free_imatrix(bh,1,lastpass-firstpass+2,1,imx); |
| free_imatrix(mw,1,lastpass-firstpass+2,1,imx); |
free_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 ------------*/ |
/*---------- Health expectancies, no variances ------------*/ |
| |
|
| strcpy(filerese,"E_"); |
strcpy(filerese,"E_"); |
|
Line 9514 Please run with mle=-1 to get a correct
|
Line 9573 Please run with mle=-1 to get a correct
|
| } |
} |
| printf("Computing Health Expectancies: result on file '%s' ...", filerese);fflush(stdout); |
printf("Computing Health Expectancies: result on file '%s' ...", filerese);fflush(stdout); |
| fprintf(ficlog,"Computing Health Expectancies: result on file '%s' ...", filerese);fflush(ficlog); |
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++){ |
for (k=1; k <= (int) pow(2,cptcoveff); k++){ |
| fprintf(ficreseij,"\n#****** "); |
fprintf(ficreseij,"\n#****** "); |
| for(j=1;j<=cptcoveff;j++) { |
for(j=1;j<=cptcoveff;j++) { |
| fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
| } |
} |
| fprintf(ficreseij,"******\n"); |
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); |
|
| |
|
| 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); |
fclose(ficreseij); |
| printf("done evsij\n");fflush(stdout); |
printf("done evsij\n");fflush(stdout); |
|
Line 9615 Please run with mle=-1 to get a correct
|
Line 9671 Please run with mle=-1 to get a correct
|
| |
|
| |
|
| for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ |
for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ |
| oldm=oldms;savm=savms; /* ZZ Segmentation fault */ |
oldm=oldms;savm=savms; /* ZZ Segmentation fault */ |
| cptcod= 0; /* To be deleted */ |
cptcod= 0; /* To be deleted */ |
| printf("varevsij %d \n",vpopbased); |
printf("varevsij %d \n",vpopbased); |
| fprintf(ficlog, "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 */ |
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 "); |
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) |
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); |
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 |
else |
| fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n"); |
fprintf(ficrest,"the age specific period (stable) prevalences in each health state \n"); |
| fprintf(ficrest,"# Age popbased mobilav e.. (std) "); |
fprintf(ficrest,"# Age popbased mobilav e.. (std) "); |
| for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); |
for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); |
| fprintf(ficrest,"\n"); |
fprintf(ficrest,"\n"); |
| /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\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); |
epj=vector(1,nlstate+1); |
| printf("Computing age specific period (stable) prevalences in each health state \n"); |
printf("Computing age specific period (stable) prevalences in each health state \n"); |
| fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n"); |
fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n"); |
| for(age=bage; age <=fage ;age++){ |
for(age=bage; age <=fage ;age++){ |
| prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k); /*ZZ Is it the correct prevalim */ |
prevalim(prlim, nlstate, p, age, oldm, savm, ftolpl, &ncvyear, k); /*ZZ Is it the correct prevalim */ |
| if (vpopbased==1) { |
if (vpopbased==1) { |
| if(mobilav ==0){ |
if(mobilav ==0){ |
| for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
| prlim[i][i]=probs[(int)age][i][k]; |
prlim[i][i]=probs[(int)age][i][k]; |
| }else{ /* mobilav */ |
}else{ /* mobilav */ |
| for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
| prlim[i][i]=mobaverage[(int)age][i][k]; |
prlim[i][i]=mobaverage[(int)age][i][k]; |
| } |
} |
| } |
} |
| |
|
| fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav); |
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 */ |
/* 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); */ |
/* printf(" age %4.0f ",age); */ |
| for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){ |
for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){ |
| for(i=1, epj[j]=0.;i <=nlstate;i++) { |
for(i=1, epj[j]=0.;i <=nlstate;i++) { |
| epj[j] += prlim[i][i]*eij[i][j][(int)age]; |
epj[j] += prlim[i][i]*eij[i][j][(int)age]; |
| /*ZZZ printf("%lf %lf ", 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]); */ |
/* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */ |
| } |
} |
| epj[nlstate+1] +=epj[j]; |
epj[nlstate+1] +=epj[j]; |
| } |
} |
| /* printf(" age %4.0f \n",age); */ |
/* printf(" age %4.0f \n",age); */ |
| |
|
| for(i=1, vepp=0.;i <=nlstate;i++) |
for(i=1, vepp=0.;i <=nlstate;i++) |
| for(j=1;j <=nlstate;j++) |
for(j=1;j <=nlstate;j++) |
| vepp += vareij[i][j][(int)age]; |
vepp += vareij[i][j][(int)age]; |
| fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp)); |
fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp)); |
| for(j=1;j <=nlstate;j++){ |
for(j=1;j <=nlstate;j++){ |
| fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age])); |
fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age])); |
| } |
} |
| fprintf(ficrest,"\n"); |
fprintf(ficrest,"\n"); |
| } |
} |
| } /* End vpopbased */ |
} /* End vpopbased */ |
| free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage); |
free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage); |
| free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage); |
free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage); |
|
Line 9722 Please run with mle=-1 to get a correct
|
Line 9778 Please run with mle=-1 to get a correct
|
| fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog); |
fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog); |
| |
|
| /*---------- End : free ----------------*/ |
/*---------- 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); |
free_ma3x(probs,1,AGESUP,1,nlstate+ndeath, 1,ncovcombmax); |
| } /* mle==-3 arrives here for freeing */ |
} /* mle==-3 arrives here for freeing */ |
| /* endfree:*/ |
/* endfree:*/ |
![[BACK]](/icons/small/back.gif){kind=icon}
Return to imach.c CVS log ![[TXT]](/icons/small/text.gif){kind=icon}
![[DIR]](/icons/small/dir.gif){kind=icon}
Up to [Local Repository] / imach / src