/* $Id$
$State$
$Log$
+ Revision 1.216 2015/12/18 17:32:11 brouard
+ Summary: 0.98r4 Warning and status=-2
+
+ Version 0.98r4 is now:
+ - displaying an error when status is -1, date of interview unknown and date of death known;
+ - permitting a status -2 when the vital status is unknown at a known date of right truncation.
+ Older changes concerning s=-2, dating from 2005 have been supersed.
+
Revision 1.215 2015/12/16 08:52:24 brouard
Summary: 0.98r4 working
double **oldm, **newm, **savm; /* Working pointers to matrices */
double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
/*FILE *fic ; */ /* Used in readdata only */
-FILE *ficpar, *ficparo,*ficres, *ficresp, *ficresphtm, *ficresphtmfr, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
+FILE *ficpar, *ficparo,*ficres, *ficresp, *ficresphtm, *ficresphtmfr, *ficrespl, *ficresplb,*ficrespij, *ficrespijb, *ficrest,*ficresf, *ficresfb,*ficrespop;
FILE *ficlog, *ficrespow;
int globpr=0; /* Global variable for printing or not */
double fretone; /* Only one call to likelihood */
FILE *ficresvpl;
char fileresvpl[FILENAMELENGTH];
char title[MAXLINE];
-char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH];
+char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH], fileresplb[FILENAMELENGTH];
char plotcmd[FILENAMELENGTH], pplotcmd[FILENAMELENGTH];
char tmpout[FILENAMELENGTH], tmpout2[FILENAMELENGTH];
char command[FILENAMELENGTH];
int outcmd=0;
-char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];
+char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filerespijb[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];
char fileresu[FILENAMELENGTH]; /* fileres without r in front */
char filelog[FILENAMELENGTH]; /* Log file */
char filerest[FILENAMELENGTH];
return prlim; /* should not reach here */
}
+
+ /**** Back Prevalence limit (stable or period prevalence) ****************/
+
+double **bprevalim(double **bprlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int *ncvyear, int ij)
+{
+ /* Computes the prevalence limit in each live state at age x by left multiplying the unit
+ matrix by transitions matrix until convergence is reached with precision ftolpl */
+ /* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1 = Wx-n Px-n ... Px-2 Px-1 I */
+ /* Wx is row vector: population in state 1, population in state 2, population dead */
+ /* or prevalence in state 1, prevalence in state 2, 0 */
+ /* newm is the matrix after multiplications, its rows are identical at a factor */
+ /* Initial matrix pimij */
+ /* {0.85204250825084937, 0.13044499163996345, 0.017512500109187184, */
+ /* 0.090851990222114765, 0.88271245433047185, 0.026435555447413338, */
+ /* 0, 0 , 1} */
+ /*
+ * and after some iteration: */
+ /* {0.45504275246439968, 0.42731458730878791, 0.11764266022681241, */
+ /* 0.45201005341706885, 0.42865420071559901, 0.11933574586733192, */
+ /* 0, 0 , 1} */
+ /* And prevalence by suppressing the deaths are close to identical rows in prlim: */
+ /* {0.51571254859325999, 0.4842874514067399, */
+ /* 0.51326036147820708, 0.48673963852179264} */
+ /* If we start from prlim again, prlim tends to a constant matrix */
+
+ int i, ii,j,k;
+ double *min, *max, *meandiff, maxmax,sumnew=0.;
+ /* double **matprod2(); */ /* test */
+ double **out, cov[NCOVMAX+1], **bmij();
+ double **newm;
+ double agefin, delaymax=200. ; /* 100 Max number of years to converge */
+ int ncvloop=0;
+
+ min=vector(1,nlstate);
+ max=vector(1,nlstate);
+ meandiff=vector(1,nlstate);
+
+ for (ii=1;ii<=nlstate+ndeath;ii++)
+ for (j=1;j<=nlstate+ndeath;j++){
+ oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+ }
+
+ cov[1]=1.;
+
+ /* Even if hstepm = 1, at least one multiplication by the unit matrix */
+ /* Start at agefin= age, computes the matrix of passage and loops decreasing agefin until convergence is reached */
+ for(agefin=age+stepm/YEARM; agefin<=age+delaymax; agefin=agefin+stepm/YEARM){
+ ncvloop++;
+ newm=savm;
+ /* Covariates have to be included here again */
+ cov[2]=agefin;
+ if(nagesqr==1)
+ cov[3]= agefin*agefin;;
+ for (k=1; k<=cptcovn;k++) {
+ /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
+ cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
+ /* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); */
+ }
+ /*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
+ /* for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]*cov[2]; */
+ for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2];
+ for (k=1; k<=cptcovprod;k++) /* Useless */
+ /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
+ cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
+
+ /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
+ /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
+ /*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/
+ /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
+ /* out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /\* Bug Valgrind *\/ */
+ out=matprod2(newm, oldm ,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate)); /* Bug Valgrind */
+
+ savm=oldm;
+ oldm=newm;
+
+ for(j=1; j<=nlstate; j++){
+ max[j]=0.;
+ min[j]=1.;
+ }
+ /* sumnew=0; */
+ /* for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k]; */
+ for(j=1; j<=nlstate; j++){
+ for(i=1;i<=nlstate;i++){
+ /* bprlim[i][j]= newm[i][j]/(1-sumnew); */
+ bprlim[i][j]= newm[i][j];
+ max[i]=FMAX(max[i],bprlim[i][j]);
+ min[i]=FMIN(min[i],bprlim[i][j]);
+ }
+ }
+
+ maxmax=0.;
+ for(i=1; i<=nlstate; i++){
+ meandiff[i]=(max[i]-min[i])/(max[i]+min[i])*2.; /* mean difference for each column */
+ maxmax=FMAX(maxmax,meandiff[i]);
+ /* printf("Back age= %d meandiff[%d]=%f, agefin=%d max[%d]=%f min[%d]=%f maxmax=%f\n", (int)age, i, meandiff[i],(int)agefin, i, max[i], i, min[i],maxmax); */
+ } /* j loop */
+ *ncvyear= -( (int)age- (int)agefin);
+ /* printf("Back maxmax=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear); */
+ if(maxmax < ftolpl){
+ printf("OK Back maxmax=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear);
+ free_vector(min,1,nlstate);
+ free_vector(max,1,nlstate);
+ free_vector(meandiff,1,nlstate);
+ return bprlim;
+ }
+ } /* age loop */
+ /* After some age loop it doesn't converge */
+ printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\
+Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);
+ /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */
+ free_vector(min,1,nlstate);
+ free_vector(max,1,nlstate);
+ free_vector(meandiff,1,nlstate);
+
+ return bprlim; /* should not reach here */
+}
+
/*************** transition probabilities ***************/
double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
return ps;
}
+/*************** transition probabilities ***************/
+
+double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
+{
+ /* According to parameters values stored in x and the covariate's values stored in cov,
+ computes the probability to be observed in state j being in state i by appying the
+ model to the ncovmodel covariates (including constant and age).
+ lnpijopii=ln(pij/pii)= aij+bij*age+cij*v1+dij*v2+... = sum_nc=1^ncovmodel xij(nc)*cov[nc]
+ and, according on how parameters are entered, the position of the coefficient xij(nc) of the
+ ncth covariate in the global vector x is given by the formula:
+ j<i nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel
+ j>=i nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel
+ Computes ln(pij/pii) (lnpijopii), deduces pij/pii by exponentiation,
+ sums on j different of i to get 1-pii/pii, deduces pii, and then all pij.
+ Outputs ps[i][j] the probability to be observed in j being in j according to
+ the values of the covariates cov[nc] and corresponding parameter values x[nc+shiftij]
+ */
+ double s1, lnpijopii;
+ /*double t34;*/
+ int i,j, nc, ii, jj;
+
+ for(i=1; i<= nlstate; i++){
+ for(j=1; j<i;j++){
+ for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
+ /*lnpijopii += param[i][j][nc]*cov[nc];*/
+ lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];
+/* printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
+ }
+ ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
+/* printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
+ }
+ for(j=i+1; j<=nlstate+ndeath;j++){
+ for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
+ /*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/
+ lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];
+/* printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */
+ }
+ ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
+ }
+ }
+
+ for(i=1; i<= nlstate; i++){
+ s1=0;
+ for(j=1; j<i; j++){
+ s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+ /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+ }
+ for(j=i+1; j<=nlstate+ndeath; j++){
+ s1+=exp(ps[i][j]); /* In fact sums pij/pii */
+ /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+ }
+ /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
+ ps[i][i]=1./(s1+1.);
+ /* Computing other pijs */
+ for(j=1; j<i; j++)
+ ps[i][j]= exp(ps[i][j])*ps[i][i];
+ for(j=i+1; j<=nlstate+ndeath; j++)
+ ps[i][j]= exp(ps[i][j])*ps[i][i];
+ /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
+ } /* end i */
+
+ for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
+ for(jj=1; jj<= nlstate+ndeath; jj++){
+ ps[ii][jj]=0;
+ ps[ii][ii]=1;
+ }
+ }
+ /* Added for backcast */
+ for(jj=1; jj<= nlstate; jj++){
+ s1=0.;
+ for(ii=1; ii<= nlstate; ii++){
+ s1+=ps[ii][jj];
+ }
+ for(ii=1; ii<= nlstate; ii++){
+ ps[ii][jj]=ps[ii][jj]/s1;
+ }
+ }
+
+ /* for(ii=1; ii<= nlstate+ndeath; ii++){ */
+ /* for(jj=1; jj<= nlstate+ndeath; jj++){ */
+ /* printf(" pmij ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */
+ /* } */
+ /* printf("\n "); */
+ /* } */
+ /* printf("\n ");printf("%lf ",cov[2]);*/
+ /*
+ for(i=1; i<= npar; i++) printf("%f ",x[i]);
+ goto end;*/
+ return ps;
+}
+
+
/**************** Product of 2 matrices ******************/
double **matprod2(double **out, double **in,int nrl, int nrh, int ncl, int nch, int ncolol, int ncoloh, double **b)
/*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath,
pmij(pmmij,cov,ncovmodel,x,nlstate));
+ /* if((int)age == 70){ */
+ /* printf(" Forward hpxij age=%d agexact=%f d=%d nhstepm=%d hstepm=%d\n", (int) age, agexact, d, nhstepm, hstepm); */
+ /* for(i=1; i<=nlstate+ndeath; i++) { */
+ /* printf("%d pmmij ",i); */
+ /* for(j=1;j<=nlstate+ndeath;j++) { */
+ /* printf("%f ",pmmij[i][j]); */
+ /* } */
+ /* printf(" oldm "); */
+ /* for(j=1;j<=nlstate+ndeath;j++) { */
+ /* printf("%f ",oldm[i][j]); */
+ /* } */
+ /* printf("\n"); */
+ /* } */
+ /* } */
savm=oldm;
oldm=newm;
}
return po;
}
+/************* Higher Back Matrix Product ***************/
+
+double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm, int ij )
+{
+ /* Computes the transition matrix starting at age 'age' over
+ 'nhstepm*hstepm*stepm' months (i.e. until
+ age (in years) age+nhstepm*hstepm*stepm/12) by multiplying
+ nhstepm*hstepm matrices.
+ Output is stored in matrix po[i][j][h] for h every 'hstepm' step
+ (typically every 2 years instead of every month which is too big
+ for the memory).
+ Model is determined by parameters x and covariates have to be
+ included manually here.
+
+ */
+
+ int i, j, d, h, k;
+ double **out, cov[NCOVMAX+1];
+ double **newm;
+ double agexact;
+ double agebegin, ageend;
+
+ /* Hstepm could be zero and should return the unit matrix */
+ for (i=1;i<=nlstate+ndeath;i++)
+ for (j=1;j<=nlstate+ndeath;j++){
+ oldm[i][j]=(i==j ? 1.0 : 0.0);
+ po[i][j][0]=(i==j ? 1.0 : 0.0);
+ }
+ /* Even if hstepm = 1, at least one multiplication by the unit matrix */
+ for(h=1; h <=nhstepm; h++){
+ for(d=1; d <=hstepm; d++){
+ newm=savm;
+ /* Covariates have to be included here again */
+ cov[1]=1.;
+ agexact=age-((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */
+ /* agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /\* age just before transition *\/ */
+ cov[2]=agexact;
+ if(nagesqr==1)
+ cov[3]= agexact*agexact;
+ for (k=1; k<=cptcovn;k++)
+ cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
+ /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
+ for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */
+ /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
+ cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
+ /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */
+ for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */
+ cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
+ /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
+
+
+ /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
+ /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
+ out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath,
+ oldm);
+ /* if((int)age == 70){ */
+ /* printf(" Backward hbxij age=%d agexact=%f d=%d nhstepm=%d hstepm=%d\n", (int) age, agexact, d, nhstepm, hstepm); */
+ /* for(i=1; i<=nlstate+ndeath; i++) { */
+ /* printf("%d pmmij ",i); */
+ /* for(j=1;j<=nlstate+ndeath;j++) { */
+ /* printf("%f ",pmmij[i][j]); */
+ /* } */
+ /* printf(" oldm "); */
+ /* for(j=1;j<=nlstate+ndeath;j++) { */
+ /* printf("%f ",oldm[i][j]); */
+ /* } */
+ /* printf("\n"); */
+ /* } */
+ /* } */
+ savm=oldm;
+ oldm=newm;
+ }
+ for(i=1; i<=nlstate+ndeath; i++)
+ for(j=1;j<=nlstate+ndeath;j++) {
+ po[i][j][h]=newm[i][j];
+ /*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/
+ }
+ /*printf("h=%d ",h);*/
+ } /* end h */
+/* printf("\n H=%d \n",h); */
+ return po;
+}
+
+
#ifdef NLOPT
double myfunc(unsigned n, const double *p1, double *grad, void *pd){
double fret;
s1=s[mw[mi][i]][i];
s2=s[mw[mi+1][i]][i];
- if(s2==-1){
- printf(" s1=%d, s2=%d i=%d \n", s1, s2, i);
- /* exit(1); */
- }
+ /* if(s2==-1){ */
+ /* printf(" s1=%d, s2=%d i=%d \n", s1, s2, i); */
+ /* /\* exit(1); *\/ */
+ /* } */
bbh=(double)bh[mi][i]/(double)stepm;
/* bias is positive if real duration
* is higher than the multiple of stepm and negative otherwise.
int i, mi, m;
/* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;
double sum=0., jmean=0.;*/
- int first, firstwo;
+ int first, firstwo, firsthree;
int j, k=0,jk, ju, jl;
double sum=0.;
first=0;
firstwo=0;
+ firsthree=0;
jmin=100000;
jmax=-1;
jmean=0.;
}
if(m >=lastpass){
if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){
- printf("Information! Unknown health status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
- fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
+ if(firsthree == 0){
+ printf("Information! Unknown health status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
+ fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
+ firsthree=1;
+ }else{
+ fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
+ }
mw[++mi][i]=m;
}
if(s[m][i]==-2){ /* Vital status is really unknown */
void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \
int lastpass, int stepm, int weightopt, char model[],\
int imx,int jmin, int jmax, double jmeanint,char rfileres[],\
- int popforecast, int prevfcast, int estepm , \
+ int popforecast, int prevfcast, int backcast, int estepm , \
double jprev1, double mprev1,double anprev1, double dateprev1, \
double jprev2, double mprev2,double anprev2, double dateprev2){
int jj1, k1, i1, cpt;
- Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ",
stepm,subdirf2(fileresu,"PIJ_"),subdirf2(fileresu,"PIJ_"));
fprintf(fichtm,"\
+ - Estimated back transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ",
+ stepm,subdirf2(fileresu,"PIJB_"),subdirf2(fileresu,"PIJB_"));
+ fprintf(fichtm,"\
- Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",
subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_"));
fprintf(fichtm,"\
+ - Period (stable) back prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",
+ subdirf2(fileresu,"PLB_"),subdirf2(fileresu,"PLB_"));
+ fprintf(fichtm,"\
- (a) Life expectancies by health status at initial age, e<sub>i.</sub> (b) health expectancies by health status at initial age, e<sub>ij</sub> . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \
<a href=\"%s\">%s</a> <br>\n",
estepm,subdirf2(fileresu,"E_"),subdirf2(fileresu,"E_"));
}
/* Period (stable) prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){
- fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \
+ fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \
<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1);
}
+ if(backcast==1){
+ /* Period (stable) back prevalence in each health state */
+ for(cpt=1; cpt<=nlstate;cpt++){
+ fprintf(fichtm,"<br>\n- Convergence to period (stable) back prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \
+<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1);
+ }
+ }
if(prevfcast==1){
/* Projection of prevalence up to period (stable) prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){
plot [%.f:%.f] ", ageminpar, agemaxpar);
k=3;
for (i=1; i<= nlstate ; i ++){
- if(i==1)
+ if(i==1){
fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
- else
+ }else{
fprintf(ficgp,", '' ");
+ }
l=(nlstate+ndeath)*(i-1)+1;
fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
for (j=2; j<= nlstate+ndeath ; j ++)
} /* end cpt state*/
} /* end covariate */
+ /* CV back preval stable (period) for each covariate */
+ 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#CV Back preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
+ for (k=1; k<=cptcoveff; k++){ /* For each covariate and each value */
+ lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
+ /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
+ /* 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,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1);
+ fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
+set ter svg size 640, 480\n\
+unset log y\n\
+plot [%.f:%.f] ", ageminpar, agemaxpar);
+ k=3; /* Offset */
+ for (i=1; i<= nlstate ; i ++){
+ if(i==1)
+ fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJB_"));
+ else
+ fprintf(ficgp,", '' ");
+ l=(nlstate+ndeath)*(i-1)+1;
+ fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); /* a vérifier */
+ for (j=2; j<= nlstate ; j ++)
+ fprintf(ficgp,"+$%d",k+l+j-1);
+ fprintf(ficgp,")) t \"prev(%d,%d)\" w l",i,cpt);
+ } /* nlstate */
+ fprintf(ficgp,"\nset out\n");
+ } /* end cpt state*/
+ } /* end covariate */
+
if(prevfcast==1){
/* Projection from cross-sectional to stable (period) for each covariate */
else
fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
}
- if(ng != 1){
- fprintf(ficgp,")/(1");
+ }else{
+ i=i-ncovmodel;
+ if(ng !=1 ) /* For logit formula of log p11 is more difficult to get */
+ fprintf(ficgp," (1.");
+ }
+
+ if(ng != 1){
+ fprintf(ficgp,")/(1");
- for(k1=1; k1 <=nlstate; k1++){
- if(nagesqr==0)
- fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
- else /* nagesqr =1 */
- fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr);
-
- ij=1;
- for(j=3; j <=ncovmodel-nagesqr; j++){
- if(ij <=cptcovage) { /* Bug valgrind */
- if((j-2)==Tage[ij]) { /* Bug valgrind */
- fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
- /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
- ij++;
- }
+ for(k1=1; k1 <=nlstate; k1++){
+ if(nagesqr==0)
+ fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
+ else /* nagesqr =1 */
+ fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr);
+
+ ij=1;
+ for(j=3; j <=ncovmodel-nagesqr; j++){
+ if(ij <=cptcovage) { /* Bug valgrind */
+ if((j-2)==Tage[ij]) { /* Bug valgrind */
+ fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
+ /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
+ ij++;
}
- else
- fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
}
- fprintf(ficgp,")");
+ else
+ fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
}
fprintf(ficgp,")");
- if(ng ==2)
- fprintf(ficgp," t \"p%d%d\" ", k2,k);
- else /* ng= 3 */
- fprintf(ficgp," t \"i%d%d\" ", k2,k);
- }else{ /* end ng <> 1 */
- fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k);
}
- if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
- i=i+ncovmodel;
+ fprintf(ficgp,")");
+ if(ng ==2)
+ fprintf(ficgp," t \"p%d%d\" ", k2,k);
+ else /* ng= 3 */
+ fprintf(ficgp," t \"i%d%d\" ", k2,k);
+ }else{ /* end ng <> 1 */
+ if( k !=k2) /* logit p11 is hard to draw */
+ fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k);
}
+ if ((k+k2)!= (nlstate*2+ndeath) && ng != 1)
+ fprintf(ficgp,",");
+ if (ng == 1 && k!=k2 && (k+k2)!= (nlstate*2+ndeath))
+ fprintf(ficgp,",");
+ i=i+ncovmodel;
} /* end k */
} /* end k2 */
fprintf(ficgp,"\n set out\n");
fprintf(ficresf," p%d%d",i,j);
fprintf(ficresf," p.%d",j);
}
- for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {
+ 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);
+ hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);
for (h=0; h<=nhstepm; h++){
if (h*hstepm/YEARM*stepm ==yearp) {
}
+/************** Back Forecasting ******************/
+void prevbackforecast(char fileres[], double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){
+ /* back1, year, month, day of starting backection
+ agemin, agemax range of age
+ dateprev1 dateprev2 range of dates during which prevalence is computed
+ anback2 year of en of backection (same day and month as back1).
+ */
+ int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1;
+ double agec; /* generic age */
+ double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
+ double *popeffectif,*popcount;
+ double ***p3mat;
+ double ***mobaverage;
+ char fileresfb[FILENAMELENGTH];
+
+ agelim=AGESUP;
+ /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people
+ in each health status at the date of interview (if between dateprev1 and dateprev2).
+ We still use firstpass and lastpass as another selection.
+ */
+ /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\ */
+ /* firstpass, lastpass, stepm, weightopt, model); */
+ prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
+
+ strcpy(fileresfb,"FB_");
+ strcat(fileresfb,fileresu);
+ if((ficresfb=fopen(fileresfb,"w"))==NULL) {
+ printf("Problem with back forecast resultfile: %s\n", fileresfb);
+ fprintf(ficlog,"Problem with back forecast resultfile: %s\n", fileresfb);
+ }
+ printf("Computing back forecasting: result on file '%s', please wait... \n", fileresfb);
+ fprintf(ficlog,"Computing back forecasting: result on file '%s', please wait... \n", fileresfb);
+
+ if (cptcoveff==0) ncodemax[cptcoveff]=1;
+
+ if (mobilav!=0) {
+ mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+ if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){
+ fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
+ printf(" Error in movingaverage mobilav=%d\n",mobilav);
+ }
+ }
+
+ stepsize=(int) (stepm+YEARM-1)/YEARM;
+ if (stepm<=12) stepsize=1;
+ if(estepm < stepm){
+ printf ("Problem %d lower than %d\n",estepm, stepm);
+ }
+ else hstepm=estepm;
+
+ hstepm=hstepm/stepm;
+ yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp and
+ fractional in yp1 */
+ anprojmean=yp;
+ yp2=modf((yp1*12),&yp);
+ mprojmean=yp;
+ yp1=modf((yp2*30.5),&yp);
+ jprojmean=yp;
+ if(jprojmean==0) jprojmean=1;
+ if(mprojmean==0) jprojmean=1;
+
+ i1=cptcoveff;
+ if (cptcovn < 1){i1=1;}
+
+ fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2);
+
+ fprintf(ficresfb,"#****** Routine prevbackforecast **\n");
+
+/* if (h==(int)(YEARM*yearp)){ */
+ for(cptcov=1, k=0;cptcov<=i1;cptcov++){
+ for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
+ k=k+1;
+ fprintf(ficresfb,"\n#****** hbijx=probability over h years, hp.jx is weighted by observed prev \n#");
+ for(j=1;j<=cptcoveff;j++) {
+ fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ }
+ fprintf(ficresfb," yearbproj age");
+ for(j=1; j<=nlstate+ndeath;j++){
+ for(i=1; i<=nlstate;i++)
+ fprintf(ficresfb," p%d%d",i,j);
+ fprintf(ficresfb," p.%d",j);
+ }
+ for (yearp=0; yearp>=(anback2-anback1);yearp -=stepsize) {
+ /* for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { */
+ fprintf(ficresfb,"\n");
+ fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp);
+ for (agec=fage; agec>=(ageminpar-1); agec--){
+ nhstepm=(int) rint((agelim-agec)*YEARM/stepm);
+ nhstepm = nhstepm/hstepm;
+ p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ oldm=oldms;savm=savms;
+ hbxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);
+
+ for (h=0; h<=nhstepm; h++){
+ if (h*hstepm/YEARM*stepm ==yearp) {
+ fprintf(ficresfb,"\n");
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficresfb,"%.f %.f ",anback1+yearp,agec+h*hstepm/YEARM*stepm);
+ }
+ for(j=1; j<=nlstate+ndeath;j++) {
+ ppij=0.;
+ for(i=1; i<=nlstate;i++) {
+ if (mobilav==1)
+ ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][cptcod];
+ else {
+ ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod];
+ }
+ if (h*hstepm/YEARM*stepm== yearp) {
+ fprintf(ficresfb," %.3f", p3mat[i][j][h]);
+ }
+ } /* end i */
+ if (h*hstepm/YEARM*stepm==yearp) {
+ fprintf(ficresfb," %.3f", ppij);
+ }
+ }/* end j */
+ } /* end h */
+ free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ } /* end agec */
+ } /* end yearp */
+ } /* end cptcod */
+ } /* end cptcov */
+
+ if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+
+ fclose(ficresfb);
+ printf("End of Computing Back forecasting \n");
+ fprintf(ficlog,"End of Computing Back forecasting\n");
+
+}
+
/************** Forecasting *****not tested NB*************/
void populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){
return 0;
}
+ int back_prevalence_limit(double *p, double **bprlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp){
+ /*--------------- Back Prevalence limit (period or stable prevalence) --------------*/
+ int i, j, k, i1 ;
+ /* double ftolpl = 1.e-10; */
+ double age, agebase, agelim;
+ double tot;
+
+ strcpy(fileresplb,"PLB_");
+ strcat(fileresplb,fileresu);
+ if((ficresplb=fopen(fileresplb,"w"))==NULL) {
+ printf("Problem with period (stable) back prevalence resultfile: %s\n", fileresplb);return 1;
+ fprintf(ficlog,"Problem with period (stable) back prevalence resultfile: %s\n", fileresplb);return 1;
+ }
+ printf("Computing period (stable) back prevalence: result on file '%s' \n", fileresplb);
+ fprintf(ficlog,"Computing period (stable) back prevalence: result on file '%s' \n", fileresplb);
+ pstamp(ficresplb);
+ fprintf(ficresplb,"# Period (stable) back prevalence. Precision given by ftolpl=%g \n", ftolpl);
+ fprintf(ficresplb,"#Age ");
+ for(i=1; i<=nlstate;i++) fprintf(ficresplb,"%d-%d ",i,i);
+ fprintf(ficresplb,"\n");
+
+ /* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */
+
+ agebase=ageminpar;
+ agelim=agemaxpar;
+
+ i1=pow(2,cptcoveff);
+ if (cptcovn < 1){i1=1;}
+
+ 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(ficresplb,"#******");
+ printf("#******");
+ fprintf(ficlog,"#******");
+ for(j=1;j<=cptcoveff;j++) {
+ fprintf(ficresplb," 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(ficresplb,"******\n");
+ printf("******\n");
+ fprintf(ficlog,"******\n");
+
+ fprintf(ficresplb,"#Age ");
+ for(j=1;j<=cptcoveff;j++) {
+ fprintf(ficresplb,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ }
+ for(i=1; i<=nlstate;i++) fprintf(ficresplb," %d-%d ",i,i);
+ fprintf(ficresplb,"Total Years_to_converge\n");
+
+ for (age=agebase; age<=agelim; age++){
+ /* for (age=agebase; age<=agebase; age++){ */
+ bprevalim(bprlim, nlstate, p, age, oldm, savm, 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)]);
+ tot=0.;
+ for(i=1; i<=nlstate;i++){
+ tot += bprlim[i][i];
+ fprintf(ficresplb," %.5f", bprlim[i][i]);
+ }
+ fprintf(ficresplb," %.3f %d\n", tot, *ncvyearp);
+ } /* Age */
+ /* was end of cptcod */
+ } /* cptcov */
+ return 0;
+}
+
int hPijx(double *p, int bage, int fage){
/*------------- h Pij x at various ages ------------*/
return 0;
}
+ int hBijx(double *p, int bage, int fage){
+ /*------------- h Bij x at various ages ------------*/
+
+ int stepsize;
+ int agelim;
+ int hstepm;
+ int nhstepm;
+ int h, i, i1, j, k;
+
+ double agedeb;
+ double ***p3mat;
+
+ strcpy(filerespijb,"PIJB_"); strcat(filerespijb,fileresu);
+ if((ficrespijb=fopen(filerespijb,"w"))==NULL) {
+ printf("Problem with Pij back resultfile: %s\n", filerespijb); return 1;
+ fprintf(ficlog,"Problem with Pij back resultfile: %s\n", filerespijb); return 1;
+ }
+ printf("Computing pij back: result on file '%s' \n", filerespijb);
+ fprintf(ficlog,"Computing pij back: result on file '%s' \n", filerespijb);
+
+ stepsize=(int) (stepm+YEARM-1)/YEARM;
+ /*if (stepm<=24) stepsize=2;*/
+
+ agelim=AGESUP;
+ hstepm=stepsize*YEARM; /* Every year of age */
+ hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */
+
+ /* hstepm=1; aff par mois*/
+ pstamp(ficrespijb);
+ fprintf(ficrespijb,"#****** h Pij x Back Probability to be in state i at age x-h being in j at x ");
+ i1= pow(2,cptcoveff);
+ /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
+ /* /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
+ /* k=k+1; */
+ for (k=1; k <= (int) pow(2,cptcoveff); k++){
+ fprintf(ficrespijb,"\n#****** ");
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficrespijb,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+ fprintf(ficrespijb,"******\n");
+
+ /* for (agedeb=fage; agedeb>=bage; agedeb--){ /\* If stepm=6 months *\/ */
+ for (agedeb=bage; agedeb<=fage; agedeb++){ /* If stepm=6 months */
+ nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
+ nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
+
+ /* nhstepm=nhstepm*YEARM; aff par mois*/
+
+ p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ oldm=oldms;savm=savms;
+ hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);
+ fprintf(ficrespijb,"# Cov Agex agex-h hpijx with i,j=");
+ for(i=1; i<=nlstate;i++)
+ for(j=1; j<=nlstate+ndeath;j++)
+ fprintf(ficrespijb," %1d-%1d",i,j);
+ fprintf(ficrespijb,"\n");
+ for (h=0; h<=nhstepm; h++){
+ /*agedebphstep = agedeb + h*hstepm/YEARM*stepm;*/
+ fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb - h*hstepm/YEARM*stepm );
+ /* fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm ); */
+ for(i=1; i<=nlstate;i++)
+ for(j=1; j<=nlstate+ndeath;j++)
+ fprintf(ficrespijb," %.5f", p3mat[i][j][h]);
+ fprintf(ficrespijb,"\n");
+ }
+ free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ fprintf(ficrespijb,"\n");
+ }
+ /*}*/
+ }
+ return 0;
+}
+
/***********************************************/
/**************** Main Program *****************/
int *tab;
int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */
+ int backcast=0;
int mobilav=0,popforecast=0;
int hstepm=0, nhstepm=0;
int agemortsup;
double bage=0, fage=110., age, agelim=0., agebase=0.;
double ftolpl=FTOL;
double **prlim;
+ double **bprlim;
double ***param; /* Matrix of parameters */
double *p;
double **matcov; /* Matrix of covariance */
double *epj, vepp;
double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;
+ double jback1=1,mback1=1,anback1=2000,jback2=1,mback2=1,anback2=2000;
+
double **ximort;
char *alph[]={"a","a","b","c","d","e"}, str[4]="1234";
int *dcwave;
fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
/* day and month of proj2 are not used but only year anproj2.*/
+ while((c=getc(ficpar))=='#' && c!= EOF){
+ ungetc(c,ficpar);
+ fgets(line, MAXLINE, ficpar);
+ fputs(line,stdout);
+ fputs(line,ficparo);
+ }
+ 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);
+ /* day and month of proj2 are not used but only year anproj2.*/
/* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint); */
printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, pathc,p);
printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt,\
- model,imx,jmin,jmax,jmean,rfileres,popforecast,prevfcast,estepm, \
+ model,imx,jmin,jmax,jmean,rfileres,popforecast,prevfcast,backcast, estepm, \
jprev1,mprev1,anprev1,dateprev1,jprev2,mprev2,anprev2,dateprev2);
/*------------ free_vector -------------*/
prevalence_limit(p, prlim, ageminpar, agemaxpar, ftolpl, &ncvyear);
fclose(ficrespl);
+ /*--------------- 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);
+ fclose(ficresplb);
+
+
#ifdef FREEEXIT2
#include "freeexit2.h"
#endif
hPijx(p, bage, fage);
fclose(ficrespij);
+ hBijx(p, bage, fage);
+ fclose(ficrespijb);
+
/*-------------- Variance of one-step probabilities---*/
k=1;
varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart);
if(prevfcast==1){
/* if(stepm ==1){*/
prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
- /* (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); */
- /* } */
}
+ if(backcast==1){
+ prevbackforecast(fileresu, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anback2, p, cptcoveff);
+ }
+ /* (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------------ */
git://henry.ined.fr / .git / commitdiff