--- imach/src/imach.c 2019/05/09 15:17:34 1.293
+++ imach/src/imach.c 2021/03/31 13:11:57 1.308
@@ -1,6 +1,69 @@
-/* $Id: imach.c,v 1.293 2019/05/09 15:17:34 brouard Exp $
+/* $Id: imach.c,v 1.308 2021/03/31 13:11:57 brouard Exp $
$State: Exp $
$Log: imach.c,v $
+ Revision 1.308 2021/03/31 13:11:57 brouard
+ Summary: Version 0.99r23
+
+
+ * imach.c (Module): Still bugs in the result loop. Thank to Holly Benett
+
+ Revision 1.307 2021/03/08 18:11:32 brouard
+ Summary: 0.99r22 fixed bug on result:
+
+ Revision 1.306 2021/02/20 15:44:02 brouard
+ Summary: Version 0.99r21
+
+ * imach.c (Module): Fix bug on quitting after result lines!
+ (Module): Version 0.99r21
+
+ Revision 1.305 2021/02/20 15:28:30 brouard
+ * imach.c (Module): Fix bug on quitting after result lines!
+
+ Revision 1.304 2021/02/12 11:34:20 brouard
+ * imach.c (Module): The use of a Windows BOM (huge) file is now an error
+
+ Revision 1.303 2021/02/11 19:50:15 brouard
+ * (Module): imach.c Someone entered 'results:' instead of 'result:'. Now it is an error which is printed.
+
+ Revision 1.302 2020/02/22 21:00:05 brouard
+ * (Module): imach.c Update mle=-3 (for computing Life expectancy
+ and life table from the data without any state)
+
+ Revision 1.301 2019/06/04 13:51:20 brouard
+ Summary: Error in 'r'parameter file backcast yearsbproj instead of yearsfproj
+
+ Revision 1.300 2019/05/22 19:09:45 brouard
+ Summary: version 0.99r19 of May 2019
+
+ Revision 1.299 2019/05/22 18:37:08 brouard
+ Summary: Cleaned 0.99r19
+
+ Revision 1.298 2019/05/22 18:19:56 brouard
+ *** empty log message ***
+
+ Revision 1.297 2019/05/22 17:56:10 brouard
+ Summary: Fix bug by moving date2dmy and nhstepm which gaefin=-1
+
+ Revision 1.296 2019/05/20 13:03:18 brouard
+ Summary: Projection syntax simplified
+
+
+ We can now start projections, forward or backward, from the mean date
+ of inteviews up to or down to a number of years of projection:
+ prevforecast=1 yearsfproj=15.3 mobil_average=0
+ or
+ prevforecast=1 starting-proj-date=1/1/2007 final-proj-date=12/31/2017 mobil_average=0
+ or
+ prevbackcast=1 yearsbproj=12.3 mobil_average=1
+ or
+ prevbackcast=1 starting-back-date=1/10/1999 final-back-date=1/1/1985 mobil_average=1
+
+ Revision 1.295 2019/05/18 09:52:50 brouard
+ Summary: doxygen tex bug
+
+ Revision 1.294 2019/05/16 14:54:33 brouard
+ Summary: There was some wrong lines added
+
Revision 1.293 2019/05/09 15:17:34 brouard
*** empty log message ***
@@ -1087,12 +1150,12 @@ typedef struct {
#define ODIRSEPARATOR '\\'
#endif
-/* $Id: imach.c,v 1.293 2019/05/09 15:17:34 brouard Exp $ */
+/* $Id: imach.c,v 1.308 2021/03/31 13:11:57 brouard Exp $ */
/* $State: Exp $ */
#include "version.h"
char version[]=__IMACH_VERSION__;
-char copyright[]="April 2018,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018";
-char fullversion[]="$Revision: 1.293 $ $Date: 2019/05/09 15:17:34 $";
+char copyright[]="March 2021,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021, INED 2000-2021";
+char fullversion[]="$Revision: 1.308 $ $Date: 2021/03/31 13:11:57 $";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
@@ -1117,7 +1180,7 @@ int nqtveff=0; /**< ntqveff number of ef
int cptcov=0; /* Working variable */
int nobs=10; /* Number of observations in the data lastobs-firstobs */
int ncovcombmax=NCOVMAX; /* Maximum calculated number of covariate combination = pow(2, cptcoveff) */
-int npar=NPARMAX;
+int npar=NPARMAX; /* Number of parameters (nlstate+ndeath-1)*nlstate*ncovmodel; */
int nlstate=2; /* Number of live states */
int ndeath=1; /* Number of dead states */
int ncovmodel=0, ncovcol=0; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */
@@ -1256,7 +1319,12 @@ double **pmmij, ***probs; /* Global poin
double ***mobaverage, ***mobaverages; /* New global variable */
double *ageexmed,*agecens;
double dateintmean=0;
+ double anprojd, mprojd, jprojd; /* For eventual projections */
+ double anprojf, mprojf, jprojf;
+ double anbackd, mbackd, jbackd; /* For eventual backprojections */
+ double anbackf, mbackf, jbackf;
+ double jintmean,mintmean,aintmean;
double *weight;
int **s; /* Status */
double *agedc;
@@ -2375,13 +2443,13 @@ void powell(double p[], double **xi, int
/* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit */
/* New value of last point Pn is not computed, P(n-1) */
for(j=1;j<=n;j++) {
- if(flatdir[j] >0){
- printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
- fprintf(ficlog," p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
- }
- /* printf("\n"); */
- /* fprintf(ficlog,"\n"); */
- }
+ if(flatdir[j] >0){
+ printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
+ fprintf(ficlog," p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
+ }
+ /* printf("\n"); */
+ /* fprintf(ficlog,"\n"); */
+ }
/* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /\* Did we reach enough precision? *\/ */
if (2.0*fabs(fp-(*fret)) <= ftol) { /* Did we reach enough precision? */
/* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */
@@ -2969,24 +3037,8 @@ double **pmij(double **ps, double *cov,
ps[ii][ii]=1;
}
}
- /* Added for backcast */ /* Transposed matrix too */
- for(jj=1; jj<= nlstate+ndeath; jj++){
- s1=0.;
- for(ii=1; ii<= nlstate+ndeath; ii++){
- s1+=ps[ii][jj];
- }
- for(ii=1; ii<= nlstate; ii++){
- ps[ii][jj]=ps[ii][jj]/s1;
- }
- }
- /* Transposition */
- for(jj=1; jj<= nlstate+ndeath; jj++){
- for(ii=jj; ii<= nlstate+ndeath; ii++){
- s1=ps[ii][jj];
- ps[ii][jj]=ps[jj][ii];
- ps[jj][ii]=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]); */
@@ -3006,7 +3058,7 @@ double **pmij(double **ps, double *cov,
/* double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, double ***dnewm, double **doldm, double **dsavm, int ij ) */
double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate, double ***prevacurrent, int ij )
{
- /* Computes the backward probability at age agefin and covariate combination ij. In fact cov is already filled and x too.
+ /* Computes the backward probability at age agefin, cov[2], and covariate combination 'ij'. In fact cov is already filled and x too.
* Call to pmij(cov and x), call to cross prevalence, sums and inverses, left multiply, and returns in **ps as well as **bmij.
*/
int i, ii, j,k;
@@ -3037,7 +3089,7 @@ double **pmij(double **ps, double *cov,
sumnew=0.;
/*for (ii=1;ii<=nlstate+ndeath;ii++){*/
for (ii=1;ii<=nlstate;ii++){ /* Only on live states */
- /* printf(" agefin=%d, ii=%d, ij=%d, prev=%f\n",(int)agefin,ii, ij, prevacurrent[(int)agefin][ii][ij]); */
+ /* printf(" agefin=%d, ii=%d, ij=%d, prev=%f\n",(int)agefin,ii, ij, prevacurrent[(int)agefin][ii][ij]); */
sumnew+=prevacurrent[(int)agefin][ii][ij];
}
if(sumnew >0.01){ /* At least some value in the prevalence */
@@ -3153,7 +3205,7 @@ double **bpmij(double **ps, double *cov,
ps[ii][ii]=1;
}
}
- /* Added for backcast */ /* Transposed matrix too */
+ /* Added for prevbcast */ /* Transposed matrix too */
for(jj=1; jj<= nlstate+ndeath; jj++){
s1=0.;
for(ii=1; ii<= nlstate+ndeath; ii++){
@@ -4412,6 +4464,20 @@ void pstamp(FILE *fichier)
fprintf(fichier,"# %s.%s\n#IMaCh version %s, %s\n#%s\n# %s", optionfilefiname,optionfilext,version,copyright, fullversion, strstart);
}
+void date2dmy(double date,double *day, double *month, double *year){
+ double yp=0., yp1=0., yp2=0.;
+
+ yp1=modf(date,&yp);/* extracts integral of date in yp and
+ fractional in yp1 */
+ *year=yp;
+ yp2=modf((yp1*12),&yp);
+ *month=yp;
+ yp1=modf((yp2*30.5),&yp);
+ *day=yp;
+ if(*day==0) *day=1;
+ if(*month==0) *month=1;
+}
+
/************ Frequencies ********************/
@@ -4986,6 +5052,7 @@ Title=%s
Datafile=%s Firstpass=%d La
}
} /* end mle=-2 */
dateintmean=dateintsum/k2cpt;
+ date2dmy(dateintmean,&jintmean,&mintmean,&aintmean);
fclose(ficresp);
fclose(ficresphtm);
@@ -5179,11 +5246,11 @@ void prevalence(double ***probs, double
void concatwav(int wav[], int **dh, int **bh, int **mw, int **s, double *agedc, double **agev, int firstpass, int lastpass, int imx, int nlstate, int stepm)
{
- /* Concatenates waves: wav[i] is the number of effective (useful waves) of individual i.
+ /* Concatenates waves: wav[i] is the number of effective (useful waves in the sense that a non interview is useless) of individual i.
Death is a valid wave (if date is known).
mw[mi][i] is the mi (mi=1 to wav[i]) effective wave of individual i
dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]
- and mw[mi+1][i]. dh depends on stepm.
+ and mw[mi+1][i]. dh depends on stepm. s[m][i] exists for any wave from firstpass to lastpass
*/
int i=0, mi=0, m=0, mli=0;
@@ -5221,10 +5288,10 @@ void concatwav(int wav[], int **dh, int
#else
if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){
if(firsthree == 0){
- printf("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 as 1-p%d%d .\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, s[m][i], nlstate+ndeath);
+ printf("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 as 1-p_{%d%d} .\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, s[m][i], nlstate+ndeath);
firsthree=1;
}
- 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 as 1-p%d%d .\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, s[m][i], nlstate+ndeath);
+ 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 as 1-p_{%d%d} .\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, s[m][i], nlstate+ndeath);
mw[++mi][i]=m;
mli=m;
}
@@ -6028,7 +6095,7 @@ void concatwav(int wav[], int **dh, int
prlim[i][i]=mobaverage[(int)age][i][ij];
}
}
- /**< Computes the shifted transition matrix \f$ {}{h}_p^{ij}_x\f$ at horizon h.
+ /**< Computes the shifted transition matrix \f$ {}{h}_p^{ij}x\f$ at horizon h.
*/
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres); /* Returns p3mat[i][j][h] for h=0 to nhstepm */
/**< And for each alive state j, sums over i \f$ w^i_x {}{h}_p^{ij}x\f$, which are the probability
@@ -6821,9 +6888,9 @@ To be simple, these graphs help to under
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 mobilav, int prevfcast, int mobilavproj, int backcast, int estepm , \
- double jprev1, double mprev1,double anprev1, double dateprev1, double dateproj1, double dateback1, \
- double jprev2, double mprev2,double anprev2, double dateprev2, double dateproj2, double dateback2){
+ int popforecast, int mobilav, int prevfcast, int mobilavproj, int prevbcast, int estepm , \
+ double jprev1, double mprev1,double anprev1, double dateprev1, double dateprojd, double dateback1, \
+ double jprev2, double mprev2,double anprev2, double dateprev2, double dateprojf, double dateback2){
int jj1, k1, i1, cpt, k4, nres;
fprintf(fichtm,"
- Result files (first order: no variance)\n \
@@ -6967,7 +7034,7 @@ divided by h: hPij
fprintf(fichtm,"
\n- Convergence to period (stable) prevalence in state %d. Or probability for a person being in state (1 to %d) at different ages, to be in state %d some years after. %s_%d-%d-%d.svg
\
", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
}
- if(backcast==1){
+ if(prevbcast==1){
/* Backward prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){
fprintf(fichtm,"
\n- Convergence to mixed (stable) back prevalence in state %d. Or probability for a person to be in state %d at a younger age, knowing that she/he was in state (1 to %d) at different older ages. %s_%d-%d-%d.svg
\
@@ -6978,10 +7045,10 @@ divided by h: hPij
/* Projection of prevalence up to period (forward stable) prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){
fprintf(fichtm,"
\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), from year %.1f up to year %.1f tending to period (stable) forward prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). %s_%d-%d-%d.svg
\
-", dateprev1, dateprev2, mobilavproj, dateproj1, dateproj2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
+", dateprev1, dateprev2, mobilavproj, dateprojd, dateprojf, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
}
}
- if(backcast==1){
+ if(prevbcast==1){
/* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){
fprintf(fichtm,"
\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), \
@@ -7094,7 +7161,7 @@ true period expectancies (those weighted
}
/******************* Gnuplot file **************/
-void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double bage, double fage , int prevfcast, int backcast, char pathc[], double p[], int offyear, int offbyear){
+void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double bage, double fage , int prevfcast, int prevbcast, char pathc[], double p[], int offyear, int offbyear){
char dirfileres[132],optfileres[132];
char gplotcondition[132], gplotlabel[132];
@@ -7248,7 +7315,7 @@ void printinggnuplot(char fileresu[], ch
} /* end covariate */
} /* end if no covariate */
- if(backcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */
+ if(prevbcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */
/* fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt); */
fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1, nres in 2 to be fixed */
if(cptcoveff ==0){
@@ -7275,7 +7342,7 @@ void printinggnuplot(char fileresu[], ch
}
} /* end covariate */
} /* end if no covariate */
- if(backcast == 1){
+ if(prevbcast == 1){
fprintf(ficgp,", \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres);
/* k1-1 error should be nres-1*/
for (i=1; i<= nlstate ; i ++) {
@@ -7294,7 +7361,7 @@ void printinggnuplot(char fileresu[], ch
}
fprintf(ficgp,"\" t\"\" w l lt 4");
} /* end if backprojcast */
- } /* end if backcast */
+ } /* end if prevbcast */
/* fprintf(ficgp,"\nset out ;unset label;\n"); */
fprintf(ficgp,"\nset out ;unset title;\n");
} /* nres */
@@ -7582,7 +7649,7 @@ set ter svg size 640, 480\nunset log y\n
/* 7eme */
- if(backcast == 1){
+ if(prevbcast == 1){
/* CV backward prevalence for each covariate */
for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
@@ -7634,7 +7701,7 @@ set ter svg size 640, 480\nunset log y\n
fprintf(ficgp,"\nset out; unset label;\n");
} /* end cpt state*/
} /* end covariate */
- } /* End if backcast */
+ } /* End if prevbcast */
/* 8eme */
if(prevfcast==1){
@@ -7750,7 +7817,7 @@ set ter svg size 640, 480\nunset log y\n
} /* end covariate */
} /* End if prevfcast */
- if(backcast==1){
+ if(prevbcast==1){
/* Back projection from cross-sectional to stable (mixed) for each covariate */
for (k1=1; k1<= m ; k1 ++) /* For each covariate combination if any */
@@ -7863,7 +7930,7 @@ set ter svg size 640, 480\nunset log y\n
fprintf(ficgp,"\nset out; unset label;\n");
} /* end cpt state*/
} /* end covariate */
- } /* End if backcast */
+ } /* End if prevbcast */
/* 9eme writing MLE parameters */
@@ -8311,16 +8378,21 @@ set ter svg size 640, 480\nunset log y\n
}/* End movingaverage */
+
/************** Forecasting ******************/
- void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double ***prev, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){
- /* proj1, year, month, day of starting projection
+/* void prevforecast(char fileres[], double dateintmean, double anprojd, double mprojd, double jprojd, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double ***prev, double bage, double fage, int firstpass, int lastpass, double anprojf, double p[], int cptcoveff)*/
+void prevforecast(char fileres[], double dateintmean, double dateprojd, double dateprojf, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double ***prev, double bage, double fage, int firstpass, int lastpass, double p[], int cptcoveff){
+ /* dateintemean, mean date of interviews
+ dateprojd, year, month, day of starting projection
+ dateprojf date of end of projection;year of end of projection (same day and month as proj1).
agemin, agemax range of age
dateprev1 dateprev2 range of dates during which prevalence is computed
- anproj2 year of en of projection (same day and month as proj1).
*/
+ /* double anprojd, mprojd, jprojd; */
+ /* double anprojf, mprojf, jprojf; */
int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
double agec; /* generic age */
- double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
+ double agelim, ppij, yp,yp1,yp2;
double *popeffectif,*popcount;
double ***p3mat;
/* double ***mobaverage; */
@@ -8357,22 +8429,27 @@ set ter svg size 640, 480\nunset log y\n
if(estepm > stepm){ /* Yes every two year */
stepsize=2;
}
+ hstepm=hstepm/stepm;
- 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;
+
+ /* yp1=modf(dateintmean,&yp);/\* extracts integral of datemean in yp and */
+ /* fractional in yp1 *\/ */
+ /* aintmean=yp; */
+ /* yp2=modf((yp1*12),&yp); */
+ /* mintmean=yp; */
+ /* yp1=modf((yp2*30.5),&yp); */
+ /* jintmean=yp; */
+ /* if(jintmean==0) jintmean=1; */
+ /* if(mintmean==0) mintmean=1; */
+
+ /* date2dmy(dateintmean,&jintmean,&mintmean,&aintmean); */
+ /* date2dmy(dateprojd,&jprojd, &mprojd, &anprojd); */
+ /* date2dmy(dateprojf,&jprojf, &mprojf, &anprojf); */
i1=pow(2,cptcoveff);
if (cptcovn < 1){i1=1;}
- fprintf(ficresf,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2);
+ fprintf(ficresf,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jintmean,mintmean,aintmean,dateintmean,dateprev1,dateprev2);
fprintf(ficresf,"#****** Routine prevforecast **\n");
@@ -8398,9 +8475,9 @@ set ter svg size 640, 480\nunset log y\n
fprintf(ficresf," p%d%d",i,j);
fprintf(ficresf," wp.%d",j);
}
- for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {
+ for (yearp=0; yearp<=(anprojf-anprojd);yearp +=stepsize) {
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 ",jprojd,mprojd,anprojd+yearp);
/* for (agec=fage; agec>=(ageminpar-1); agec--){ */
for (agec=fage; agec>=(bage); agec--){
nhstepm=(int) rint((agelim-agec)*YEARM/stepm);
@@ -8418,7 +8495,7 @@ set ter svg size 640, 480\nunset log y\n
fprintf(ficresf,"\n");
for(j=1;j<=cptcoveff;j++)
fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
- fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm);
+ fprintf(ficresf,"%.f %.f ",anprojd+yearp,agec+h*hstepm/YEARM*stepm);
for(j=1; j<=nlstate+ndeath;j++) {
ppij=0.;
@@ -8446,8 +8523,9 @@ set ter svg size 640, 480\nunset log y\n
}
/************** Back Forecasting ******************/
- void prevbackforecast(char fileres[], double ***prevacurrent, double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){
- /* back1, year, month, day of starting backection
+ /* void prevbackforecast(char fileres[], double ***prevacurrent, double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){ */
+ void prevbackforecast(char fileres[], double ***prevacurrent, double dateintmean, double dateprojd, double dateprojf, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double p[], int cptcoveff){
+ /* back1, year, month, day of starting backprojection
agemin, agemax range of age
dateprev1 dateprev2 range of dates during which prevalence is computed
anback2 year of end of backprojection (same day and month as back1).
@@ -8455,7 +8533,7 @@ set ter svg size 640, 480\nunset log y\n
*/
int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
double agec; /* generic age */
- double agelim, ppij, ppi, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
+ double agelim, ppij, ppi, yp,yp1,yp2; /* ,jintmean,mintmean,aintmean;*/
double *popeffectif,*popcount;
double ***p3mat;
/* double ***mobaverage; */
@@ -8498,21 +8576,21 @@ set ter svg size 640, 480\nunset log y\n
}
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;
+ /* yp1=modf(dateintmean,&yp);/\* extracts integral of datemean in yp and */
+ /* fractional in yp1 *\/ */
+ /* aintmean=yp; */
+ /* yp2=modf((yp1*12),&yp); */
+ /* mintmean=yp; */
+ /* yp1=modf((yp2*30.5),&yp); */
+ /* jintmean=yp; */
+ /* if(jintmean==0) jintmean=1; */
+ /* if(mintmean==0) jintmean=1; */
i1=pow(2,cptcoveff);
if (cptcovn < 1){i1=1;}
- fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2);
- printf("# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2);
+ fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jintmean,mintmean,aintmean,dateintmean,dateprev1,dateprev2);
+ printf("# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jintmean,mintmean,aintmean,dateintmean,dateprev1,dateprev2);
fprintf(ficresfb,"#****** Routine prevbackforecast **\n");
@@ -8537,10 +8615,10 @@ set ter svg size 640, 480\nunset log y\n
fprintf(ficresfb," b%d%d",i,j);
fprintf(ficresfb," b.%d",j);
}
- for (yearp=0; yearp>=(anback2-anback1);yearp -=stepsize) {
+ for (yearp=0; yearp>=(anbackf-anbackd);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);
+ fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jbackd,mbackd,anbackd+yearp);
/* printf("\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp); */
/* for (agec=bage; agec<=agemax-1; agec++){ /\* testing *\/ */
for (agec=bage; agec<=fage; agec++){ /* testing */
@@ -8563,7 +8641,7 @@ set ter svg size 640, 480\nunset log y\n
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);
+ fprintf(ficresfb,"%.f %.f ",anbackd+yearp,agec-h*hstepm/YEARM*stepm);
for(i=1; i<=nlstate+ndeath;i++) {
ppij=0.;ppi=0.;
for(j=1; j<=nlstate;j++) {
@@ -9011,7 +9089,7 @@ void prwizard(int ncovmodel, int nlstate
/******************* Gompertz Likelihood ******************************/
double gompertz(double x[])
{
- double A,B,L=0.0,sump=0.,num=0.;
+ double A=0.0,B=0.,L=0.0,sump=0.,num=0.;
int i,n=0; /* n is the size of the sample */
for (i=1;i<=imx ; i++) {
@@ -9019,28 +9097,34 @@ double gompertz(double x[])
/* sump=sump+1;*/
num=num+1;
}
-
-
+ L=0.0;
+ /* agegomp=AGEGOMP; */
/* for (i=0; i<=imx; i++)
if (wav[i]>0) printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d %d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]);*/
- for (i=1;i<=imx ; i++)
- {
- if (cens[i] == 1 && wav[i]>1)
- A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)));
-
- if (cens[i] == 0 && wav[i]>1)
+ for (i=1;i<=imx ; i++) {
+ /* mu(a)=mu(agecomp)*exp(teta*(age-agegomp))
+ mu(a)=x[1]*exp(x[2]*(age-agegomp)); x[1] and x[2] are per year.
+ * L= Product mu(agedeces)exp(-\int_ageexam^agedc mu(u) du ) for a death between agedc (in month)
+ * and agedc +1 month, cens[i]=0: log(x[1]/YEARM)
+ * +
+ * exp(-\int_ageexam^agecens mu(u) du ) when censored, cens[i]=1
+ */
+ if (wav[i] > 1 || agedc[i] < AGESUP) {
+ if (cens[i] == 1){
+ A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)));
+ } else if (cens[i] == 0){
A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)))
- +log(x[1]/YEARM)+x[2]*(agedc[i]-agegomp)+log(YEARM);
-
+ +log(x[1]/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM);
+ } else
+ printf("Gompertz cens[%d] neither 1 nor 0\n",i);
/*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */
- if (wav[i] > 1 ) { /* ??? */
- L=L+A*weight[i];
+ L=L+A*weight[i];
/* printf("\ni=%d A=%f L=%lf x[1]=%lf x[2]=%lf ageex=%lf agecens=%lf cens=%d agedc=%lf weight=%lf\n",i,A,L,x[1],x[2],ageexmed[i]*12,agecens[i]*12,cens[i],agedc[i]*12,weight[i]);*/
- }
- }
+ }
+ }
- /*printf("x1=%2.9f x2=%2.9f x3=%2.9f L=%f\n",x[1],x[2],x[3],L);*/
+ /*printf("x1=%2.9f x2=%2.9f x3=%2.9f L=%f\n",x[1],x[2],x[3],L);*/
return -2*L*num/sump;
}
@@ -9049,7 +9133,7 @@ double gompertz(double x[])
/******************* Gompertz_f Likelihood ******************************/
double gompertz_f(const gsl_vector *v, void *params)
{
- double A,B,LL=0.0,sump=0.,num=0.;
+ double A=0.,B=0.,LL=0.0,sump=0.,num=0.;
double *x= (double *) v->data;
int i,n=0; /* n is the size of the sample */
@@ -9142,6 +9226,7 @@ int readdata(char datafile[], int firsto
int i=0, j=0, n=0, iv=0, v;
int lstra;
int linei, month, year,iout;
+ int noffset=0; /* This is the offset if BOM data file */
char line[MAXLINE], linetmp[MAXLINE];
char stra[MAXLINE], strb[MAXLINE];
char *stratrunc;
@@ -9175,8 +9260,53 @@ int readdata(char datafile[], int firsto
fprintf(ficlog,"Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(ficlog);return 1;
}
- i=1;
+ /* Is it a BOM UTF-8 Windows file? */
+ /* First data line */
linei=0;
+ while(fgets(line, MAXLINE, fic)) {
+ noffset=0;
+ if( line[0] == (char)0xEF && line[1] == (char)0xBB) /* EF BB BF */
+ {
+ noffset=noffset+3;
+ printf("# Data file '%s' is an UTF8 BOM file, please convert to UTF8 or ascii file and rerun.\n",datafile);fflush(stdout);
+ fprintf(ficlog,"# Data file '%s' is an UTF8 BOM file, please convert to UTF8 or ascii file and rerun.\n",datafile);
+ fflush(ficlog); return 1;
+ }
+ /* else if( line[0] == (char)0xFE && line[1] == (char)0xFF)*/
+ else if( line[0] == (char)0xFF && line[1] == (char)0xFE)
+ {
+ noffset=noffset+2;
+ printf("# Error Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile);fflush(stdout);
+ fprintf(ficlog,"# Error Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile);
+ fflush(ficlog); return 1;
+ }
+ else if( line[0] == 0 && line[1] == 0)
+ {
+ if( line[2] == (char)0xFE && line[3] == (char)0xFF){
+ noffset=noffset+4;
+ printf("# Error Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile);fflush(stdout);
+ fprintf(ficlog,"# Error Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile);
+ fflush(ficlog); return 1;
+ }
+ } else{
+ ;/*printf(" Not a BOM file\n");*/
+ }
+ /* If line starts with a # it is a comment */
+ if (line[noffset] == '#') {
+ linei=linei+1;
+ break;
+ }else{
+ break;
+ }
+ }
+ fclose(fic);
+ if((fic=fopen(datafile,"r"))==NULL) {
+ printf("Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(stdout);
+ fprintf(ficlog,"Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(ficlog);return 1;
+ }
+ /* Not a Bom file */
+
+ i=1;
while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs))) {
linei=linei+1;
for(j=strlen(line); j>=0;j--){ /* Untabifies line */
@@ -9297,7 +9427,11 @@ int readdata(char datafile[], int firsto
return 1;
}
anint[j][i]= (double) year;
- mint[j][i]= (double)month;
+ mint[j][i]= (double)month;
+ /* if( (int)anint[j][i]+ (int)(mint[j][i])/12. < (int) (moisnais[i]/12.+annais[i])){ */
+ /* printf("Warning reading data around '%s' at line number %d for individual %d, '%s'\nThe date of interview (%2d/%4d) at wave %d occurred before the date of birth (%2d/%4d).\n",strb, linei,i, line, mint[j][i],anint[j][i], moisnais[i],annais[i]); */
+ /* fprintf(ficlog,"Warning reading data around '%s' at line number %d for individual %d, '%s'\nThe date of interview (%2d/%4d) at wave %d occurred before the date of birth (%2d/%4d).\n",strb, linei,i, line, mint[j][i],anint[j][i], moisnais[i],annais[i]); */
+ /* } */
strcpy(line,stra);
} /* End loop on waves */
@@ -9336,7 +9470,14 @@ int readdata(char datafile[], int firsto
}
annais[i]=(double)(year);
- moisnais[i]=(double)(month);
+ moisnais[i]=(double)(month);
+ for (j=1;j<=maxwav;j++){
+ if( (int)anint[j][i]+ (int)(mint[j][i])/12. < (int) (moisnais[i]/12.+annais[i])){
+ printf("Warning reading data around '%s' at line number %d for individual %d, '%s'\nThe date of interview (%2d/%4d) at wave %d occurred before the date of birth (%2d/%4d).\n",strb, linei,i, line, (int)mint[j][i],(int)anint[j][i], j,(int)moisnais[i],(int)annais[i]);
+ fprintf(ficlog,"Warning reading data around '%s' at line number %d for individual %d, '%s'\nThe date of interview (%2d/%4d) at wave %d occurred before the date of birth (%2d/%4d).\n",strb, linei,i, line, (int)mint[j][i],(int)anint[j][i], j, (int)moisnais[i],(int)annais[i]);
+ }
+ }
+
strcpy(line,stra);
/* Sample weight */
@@ -9459,7 +9600,6 @@ int decoderesult ( char resultline[], in
char stra[80], strb[80], strc[80], strd[80],stre[80];
removefirstspace(&resultline);
- printf("decoderesult:%s\n",resultline);
if (strstr(resultline,"v") !=0){
printf("Error. 'v' must be in upper case 'V' result: %s ",resultline);
@@ -9474,7 +9614,6 @@ int decoderesult ( char resultline[], in
TKresult[nres]=0; /* Combination for the nresult and the model */
return (0);
}
-
if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */
printf("ERROR: the number of variable in the resultline, %d, differs from the number of variable used in the model line, %d.\n",j, cptcovs);
fprintf(ficlog,"ERROR: the number of variable in the resultline, %d, differs from the number of variable used in the model line, %d.\n",j, cptcovs);
@@ -10713,8 +10852,8 @@ int hPijx(double *p, int bage, int fage)
/* for (agedeb=fage; agedeb>=bage; agedeb--){ /\* If stepm=6 months *\/ */
for (agedeb=bage; agedeb<=fage; agedeb++){ /* If stepm=6 months and estepm=24 (2 years) */
/* nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /\* Typically 20 years = 20*12/6=40 *\/ */
- nhstepm=(int) rint((agedeb-ageminl)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
- nhstepm = nhstepm/hstepm; /* Typically 40/4=10, because estepm=24 stepm=6 => hstepm=24/6=4 */
+ nhstepm=(int) rint((agedeb-ageminl)*YEARM/stepm+0.1)-1; /* Typically 20 years = 20*12/6=40 or 55*12/24=27.5-1.1=>27 */
+ nhstepm = nhstepm/hstepm; /* Typically 40/4=10, because estepm=24 stepm=6 => hstepm=24/6=4 or 28*/
/* nhstepm=nhstepm*YEARM; aff par mois*/
@@ -10806,7 +10945,14 @@ int main(int argc, char *argv[])
int *tab;
int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */
- int backcast=0; /* defined as global for mlikeli and mle*/
+ /* double anprojd, mprojd, jprojd; /\* For eventual projections *\/ */
+ /* double anprojf, mprojf, jprojf; */
+ /* double jintmean,mintmean,aintmean; */
+ int prvforecast = 0; /* Might be 1 (date of beginning of projection is a choice or 2 is the dateintmean */
+ int prvbackcast = 0; /* Might be 1 (date of beginning of projection is a choice or 2 is the dateintmean */
+ double yrfproj= 10.0; /* Number of years of forward projections */
+ double yrbproj= 10.0; /* Number of years of backward projections */
+ int prevbcast=0; /* defined as global for mlikeli and mle, replacing backcast */
int mobilav=0,popforecast=0;
int hstepm=0, nhstepm=0;
int agemortsup;
@@ -10831,8 +10977,9 @@ int main(int argc, char *argv[])
double *epj, vepp;
double dateprev1, dateprev2;
- double jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000, dateproj1=0, dateproj2=0;
- double jback1=1,mback1=1,anback1=2000,jback2=1,mback2=1,anback2=2000, dateback1=0, dateback2=0;
+ double jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000, dateproj1=0, dateproj2=0, dateprojd=0, dateprojf=0;
+ double jback1=1,mback1=1,anback1=2000,jback2=1,mback2=1,anback2=2000, dateback1=0, dateback2=0, datebackd=0, datebackf=0;
+
double **ximort;
char *alph[]={"a","a","b","c","d","e"}, str[4]="1234";
@@ -11027,7 +11174,8 @@ int main(int argc, char *argv[])
noffset=noffset+3;
printf("# File is an UTF8 Bom.\n"); // 0xBF
}
- else if( line[0] == (char)0xFE && line[1] == (char)0xFF)
+/* else if( line[0] == (char)0xFE && line[1] == (char)0xFF)*/
+ else if( line[0] == (char)0xFF && line[1] == (char)0xFE)
{
noffset=noffset+2;
printf("# File is an UTF16BE BOM file\n");
@@ -11115,8 +11263,8 @@ int main(int argc, char *argv[])
}
if((num_filled=sscanf(line,"model=1+age%[^.\n]", model)) !=EOF){
if (num_filled != 1){
- printf("ERROR %d: Model should be at minimum 'model=1+age' %s\n",num_filled, line);
- fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age' %s\n",num_filled, line);
+ printf("ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line);
+ fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line);
model[0]='\0';
goto end;
}
@@ -11802,8 +11950,8 @@ Interval (in months) between two waves:
ximort[i][j]=(i == j ? 1.0 : 0.0);
}
- /*p[1]=0.0268; p[NDIM]=0.083;*/
- /*printf("%lf %lf", p[1], p[2]);*/
+ p[1]=0.0268; p[NDIM]=0.083;
+ /* printf("%lf %lf", p[1], p[2]); */
#ifdef GSL
@@ -11929,9 +12077,9 @@ Interval (in months) between two waves:
printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i]));
fprintf(ficlog,"%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i]));
}
- lsurv=vector(1,AGESUP);
- lpop=vector(1,AGESUP);
- tpop=vector(1,AGESUP);
+ lsurv=vector(agegomp,AGESUP);
+ lpop=vector(agegomp,AGESUP);
+ tpop=vector(agegomp,AGESUP);
lsurv[agegomp]=100000;
for (k=agegomp;k<=AGESUP;k++) {
@@ -11978,9 +12126,9 @@ Please run with mle=-1 to get a correct
stepm, weightopt,\
model,imx,p,matcov,agemortsup);
- free_vector(lsurv,1,AGESUP);
- free_vector(lpop,1,AGESUP);
- free_vector(tpop,1,AGESUP);
+ free_vector(lsurv,agegomp,AGESUP);
+ free_vector(lpop,agegomp,AGESUP);
+ free_vector(tpop,agegomp,AGESUP);
free_matrix(ximort,1,NDIM,1,NDIM);
free_ivector(dcwave,firstobs,lastobs);
free_vector(agecens,firstobs,lastobs);
@@ -12181,6 +12329,7 @@ Please run with mle=-1 to get a correct
fputs(line,stdout);
fputs(line,ficparo);
fputs(line,ficlog);
+ fputs(line,ficres);
continue;
}else
break;
@@ -12226,6 +12375,7 @@ Please run with mle=-1 to get a correct
fputs(line,stdout);
fputs(line,ficparo);
fputs(line,ficlog);
+ fputs(line,ficres);
continue;
}else
break;
@@ -12251,6 +12401,7 @@ Please run with mle=-1 to get a correct
fputs(line,stdout);
fputs(line,ficparo);
fputs(line,ficlog);
+ fputs(line,ficres);
continue;
}else
break;
@@ -12273,95 +12424,108 @@ Please run with mle=-1 to get a correct
}
/* Results */
+ endishere=0;
nresult=0;
+ parameterline=0;
do{
if(!fgets(line, MAXLINE, ficpar)){
endishere=1;
- parameterline=14;
+ parameterline=15;
}else if (line[0] == '#') {
/* If line starts with a # it is a comment */
numlinepar++;
fputs(line,stdout);
fputs(line,ficparo);
fputs(line,ficlog);
+ fputs(line,ficres);
continue;
}else if(sscanf(line,"prevforecast=%[^\n]\n",modeltemp))
parameterline=11;
- else if(sscanf(line,"backcast=%[^\n]\n",modeltemp))
+ else if(sscanf(line,"prevbackcast=%[^\n]\n",modeltemp))
parameterline=12;
- else if(sscanf(line,"result:%[^\n]\n",modeltemp))
+ else if(sscanf(line,"result:%[^\n]\n",modeltemp)){
parameterline=13;
+ }
else{
parameterline=14;
}
- switch (parameterline){
+ switch (parameterline){ /* =0 only if only comments */
case 11:
- if((num_filled=sscanf(line,"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)) !=EOF){
- if (num_filled != 8) {
- printf("Error: Not 8 (data)parameters in line but %d, for example:prevforecast=1 starting-proj-date=1/1/1990 final-proj-date=1/1/2000 mobil_average=0\n, your line=%s . Probably you are running an older format.\n",num_filled,line);
- fprintf(ficlog,"Error: Not 8 (data)parameters in line but %d, for example:prevforecast=1 starting-proj-date=1/1/1990 final-proj-date=1/1/2000 mov_average=0\n, your line=%s . Probably you are running an older format.\n",num_filled,line);
- goto end;
- }
- fprintf(ficparo,"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);
+ if((num_filled=sscanf(line,"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)) !=EOF && (num_filled == 8)){
+ fprintf(ficparo,"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);
printf("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);
fprintf(ficlog,"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);
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.*/
dateproj1=anproj1+(mproj1-1)/12.+(jproj1-1)/365.;
dateproj2=anproj2+(mproj2-1)/12.+(jproj2-1)/365.;
-
+ prvforecast = 1;
+ }
+ else if((num_filled=sscanf(line,"prevforecast=%d yearsfproj=%lf mobil_average=%d\n",&prevfcast,&yrfproj,&mobilavproj)) !=EOF){/* && (num_filled == 3))*/
+ printf("prevforecast=%d yearsfproj=%lf.2 mobil_average=%d\n",prevfcast,yrfproj,mobilavproj);
+ fprintf(ficlog,"prevforecast=%d yearsfproj=%lf.2 mobil_average=%d\n",prevfcast,yrfproj,mobilavproj);
+ fprintf(ficres,"prevforecast=%d yearsfproj=%lf.2 mobil_average=%d\n",prevfcast,yrfproj,mobilavproj);
+ prvforecast = 2;
+ }
+ else {
+ printf("Error: Not 8 (data)parameters in line but %d, for example:prevforecast=1 starting-proj-date=1/1/1990 final-proj-date=1/1/2000 mobil_average=0\nnor 3 (data)parameters, for example:prevforecast=1 yearsfproj=10 mobil_average=0. Your line=%s . You are running probably an older format.\n, ",num_filled,line);
+ fprintf(ficlog,"Error: Not 8 (data)parameters in line but %d, for example:prevforecast=1 starting-proj-date=1/1/1990 final-proj-date=1/1/2000 mobil_average=0\nnor 3 (data)parameters, for example:prevforecast=1 yearproj=10 mobil_average=0. Your line=%s . You are running probably an older format.\n, ",num_filled,line);
+ goto end;
}
break;
case 12:
- /*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);*/
- if((num_filled=sscanf(line,"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)) !=EOF){
- if (num_filled != 8) {
- printf("Error: Not 8 (data)parameters in line but %d, for example:backcast=1 starting-back-date=1/1/1990 final-back-date=1/1/1970 mobil_average=1\n, your line=%s . Probably you are running an older format.\n",num_filled,line);
- fprintf(ficlog,"Error: Not 8 (data)parameters in line but %d, for example:backcast=1 starting-back-date=1/1/1990 final-back-date=1/1/1970 mobil_average=1\n, your line=%s . Probably you are running an older format.\n",num_filled,line);
- goto end;
- }
- printf("backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
- fprintf(ficparo,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
- fprintf(ficlog,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
- fprintf(ficres,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
- /* day and month of proj2 are not used but only year anproj2.*/
+ if((num_filled=sscanf(line,"prevbackcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&prevbcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj)) !=EOF && (num_filled == 8)){
+ fprintf(ficparo,"prevbackcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevbcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
+ printf("prevbackcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevbcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
+ fprintf(ficlog,"prevbackcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevbcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
+ fprintf(ficres,"prevbackcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevbcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
+ /* day and month of back2 are not used but only year anback2.*/
dateback1=anback1+(mback1-1)/12.+(jback1-1)/365.;
dateback2=anback2+(mback2-1)/12.+(jback2-1)/365.;
+ prvbackcast = 1;
+ }
+ else if((num_filled=sscanf(line,"prevbackcast=%d yearsbproj=%lf mobil_average=%d\n",&prevbcast,&yrbproj,&mobilavproj)) ==3){/* && (num_filled == 3))*/
+ printf("prevbackcast=%d yearsbproj=%lf.2 mobil_average=%d\n",prevbcast,yrbproj,mobilavproj);
+ fprintf(ficlog,"prevbackcast=%d yearsbproj=%lf.2 mobil_average=%d\n",prevbcast,yrbproj,mobilavproj);
+ fprintf(ficres,"prevbackcast=%d yearsbproj=%lf.2 mobil_average=%d\n",prevbcast,yrbproj,mobilavproj);
+ prvbackcast = 2;
+ }
+ else {
+ printf("Error: Not 8 (data)parameters in line but %d, for example:prevbackcast=1 starting-back-date=1/1/1990 final-back-date=1/1/2000 mobil_average=0\nnor 3 (data)parameters, for example:prevbackcast=1 yearsbproj=10 mobil_average=0. Your line=%s . You are running probably an older format.\n, ",num_filled,line);
+ fprintf(ficlog,"Error: Not 8 (data)parameters in line but %d, for example:prevbackcast=1 starting-back-date=1/1/1990 final-back-date=1/1/2000 mobil_average=0\nnor 3 (data)parameters, for example:prevbackcast=1 yearbproj=10 mobil_average=0. Your line=%s . You are running probably an older format.\n, ",num_filled,line);
+ goto end;
}
break;
case 13:
- if((num_filled=sscanf(line,"result:%[^\n]\n",resultline)) !=EOF){
- if (num_filled == 0){
- resultline[0]='\0';
- printf("Warning %d: no result line! It should be at minimum 'result: V2=0 V1=1 or result:.\n%s\n", num_filled, line);
- fprintf(ficlog,"Warning %d: no result line! It should be at minimum 'result: V2=0 V1=1 or result:.\n%s\n", num_filled, line);
- break;
- } else if (num_filled != 1){
- printf("ERROR %d: result line! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",num_filled, line);
- fprintf(ficlog,"ERROR %d: result line! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",num_filled, line);
- }
- nresult++; /* Sum of resultlines */
- printf("Result %d: result=%s\n",nresult, resultline);
- if(nresult > MAXRESULTLINES){
- printf("ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\n",MAXRESULTLINES,nresult);
- fprintf(ficlog,"ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\n",MAXRESULTLINES,nresult);
- goto end;
- }
- decoderesult(resultline, nresult); /* Fills TKresult[nresult] combination and Tresult[nresult][k4+1] combination values */
- fprintf(ficparo,"result: %s\n",resultline);
- fprintf(ficres,"result: %s\n",resultline);
- fprintf(ficlog,"result: %s\n",resultline);
- break;
- case 14:
- if(ncovmodel >2 && nresult==0 ){
- printf("ERROR: no result lines! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",line);
- goto end;
- }
- break;
- default:
- nresult=1;
- decoderesult(".",nresult ); /* No covariate */
+ num_filled=sscanf(line,"result:%[^\n]\n",resultline);
+ nresult++; /* Sum of resultlines */
+ printf("Result %d: result:%s\n",nresult, resultline);
+ if(nresult > MAXRESULTLINES){
+ printf("ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\nYou can use the 'r' parameter file '%s' which uses option mle=0 to get other results. ",MAXRESULTLINES,nresult,rfileres);
+ fprintf(ficlog,"ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\nYou can use the 'r' parameter file '%s' which uses option mle=0 to get other results. ",MAXRESULTLINES,nresult,rfileres);
+ goto end;
+ }
+ decoderesult(resultline, nresult); /* Fills TKresult[nresult] combination and Tresult[nresult][k4+1] combination values */
+ fprintf(ficparo,"result: %s\n",resultline);
+ fprintf(ficres,"result: %s\n",resultline);
+ fprintf(ficlog,"result: %s\n",resultline);
+ break;
+ case 14:
+ printf("Error: Unknown command '%s'\n",line);
+ fprintf(ficlog,"Error: Unknown command '%s'\n",line);
+ if(ncovmodel >=2 && nresult==0 ){
+ printf("ERROR: no result lines! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",line);
+ fprintf(ficlog,"ERROR: no result lines! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",line);
}
+ /* goto end; */
+ break;
+ case 15:
+ printf("End of resultlines.\n");
+ fprintf(ficlog,"End of resultlines.\n");
+ break;
+ default: /* parameterline =0 */
+ nresult=1;
+ decoderesult(".",nresult ); /* No covariate */
} /* End switch parameterline */
}while(endishere==0); /* End do */
@@ -12378,11 +12542,44 @@ This is probably because your parameter
Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
}else{
/* printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, backcast, pathc,p, (int)anproj1-(int)agemin, (int)anback1-(int)agemax+1); */
- printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,bage, fage, prevfcast, backcast, pathc,p, (int)anproj1-bage, (int)anback1-fage);
+ /* It seems that anprojd which is computed from the mean year at interview which is known yet because of freqsummary */
+ /* date2dmy(dateintmean,&jintmean,&mintmean,&aintmean); */ /* Done in freqsummary */
+ if(prvforecast==1){
+ dateprojd=(jproj1+12*mproj1+365*anproj1)/365;
+ jprojd=jproj1;
+ mprojd=mproj1;
+ anprojd=anproj1;
+ dateprojf=(jproj2+12*mproj2+365*anproj2)/365;
+ jprojf=jproj2;
+ mprojf=mproj2;
+ anprojf=anproj2;
+ } else if(prvforecast == 2){
+ dateprojd=dateintmean;
+ date2dmy(dateprojd,&jprojd, &mprojd, &anprojd);
+ dateprojf=dateintmean+yrfproj;
+ date2dmy(dateprojf,&jprojf, &mprojf, &anprojf);
+ }
+ if(prvbackcast==1){
+ datebackd=(jback1+12*mback1+365*anback1)/365;
+ jbackd=jback1;
+ mbackd=mback1;
+ anbackd=anback1;
+ datebackf=(jback2+12*mback2+365*anback2)/365;
+ jbackf=jback2;
+ mbackf=mback2;
+ anbackf=anback2;
+ } else if(prvbackcast == 2){
+ datebackd=dateintmean;
+ date2dmy(datebackd,&jbackd, &mbackd, &anbackd);
+ datebackf=dateintmean-yrbproj;
+ date2dmy(datebackf,&jbackf, &mbackf, &anbackf);
+ }
+
+ printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,bage, fage, prevfcast, prevbcast, pathc,p, (int)anprojd-bage, (int)anbackd-fage);
}
printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt, \
- model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,backcast, estepm, \
- jprev1,mprev1,anprev1,dateprev1, dateproj1, dateback1,jprev2,mprev2,anprev2,dateprev2,dateproj2, dateback2);
+ model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,prevbcast, estepm, \
+ jprev1,mprev1,anprev1,dateprev1, dateprojd, datebackd,jprev2,mprev2,anprev2,dateprev2,dateprojf, datebackf);
/*------------ free_vector -------------*/
/* chdir(path); */
@@ -12455,13 +12652,23 @@ Please run with mle=-1 to get a correct
}/* end if moving average */
/*---------- Forecasting ------------------*/
- if(prevfcast==1){
- /* if(stepm ==1){*/
- prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, mobaverage, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
+ if(prevfcast==1){
+ /* /\* if(stepm ==1){*\/ */
+ /* /\* anproj1, mproj1, jproj1 either read explicitly or yrfproj *\/ */
+ /*This done previously after freqsummary.*/
+ /* dateprojd=(jproj1+12*mproj1+365*anproj1)/365; */
+ /* dateprojf=(jproj2+12*mproj2+365*anproj2)/365; */
+
+ /* } else if (prvforecast==2){ */
+ /* /\* if(stepm ==1){*\/ */
+ /* /\* anproj1, mproj1, jproj1 either read explicitly or yrfproj *\/ */
+ /* } */
+ /*prevforecast(fileresu, dateintmean, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, mobaverage, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);*/
+ prevforecast(fileresu,dateintmean, dateprojd, dateprojf, agemin, agemax, dateprev1, dateprev2, mobilavproj, mobaverage, bage, fage, firstpass, lastpass, p, cptcoveff);
}
- /* Backcasting */
- if(backcast==1){
+ /* Prevbcasting */
+ if(prevbcast==1){
ddnewms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);
ddoldms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);
ddsavms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);
@@ -12476,8 +12683,14 @@ Please run with mle=-1 to get a correct
hBijx(p, bage, fage, mobaverage);
fclose(ficrespijb);
- prevbackforecast(fileresu, mobaverage, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2,
- mobilavproj, bage, fage, firstpass, lastpass, anback2, p, cptcoveff);
+ /* /\* prevbackforecast(fileresu, mobaverage, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, *\/ */
+ /* /\* mobilavproj, bage, fage, firstpass, lastpass, anback2, p, cptcoveff); *\/ */
+ /* prevbackforecast(fileresu, mobaverage, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, */
+ /* mobilavproj, bage, fage, firstpass, lastpass, anback2, p, cptcoveff); */
+ prevbackforecast(fileresu, mobaverage, dateintmean, dateprojd, dateprojf, agemin, agemax, dateprev1, dateprev2,
+ mobilavproj, bage, fage, firstpass, lastpass, p, cptcoveff);
+
+
varbprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, bprlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff);
@@ -12485,7 +12698,7 @@ Please run with mle=-1 to get a correct
free_matrix(ddnewms, 1, nlstate+ndeath, 1, nlstate+ndeath);
free_matrix(ddsavms, 1, nlstate+ndeath, 1, nlstate+ndeath);
free_matrix(ddoldms, 1, nlstate+ndeath, 1, nlstate+ndeath);
- } /* end Backcasting */
+ } /* end Prevbcasting */
/* ------ Other prevalence ratios------------ */