--- imach/src/imach.c 2018/05/02 20:58:27 1.288 +++ imach/src/imach.c 2020/02/22 21:00:05 1.302 @@ -1,6 +1,62 @@ -/* $Id: imach.c,v 1.288 2018/05/02 20:58:27 brouard Exp $ +/* $Id: imach.c,v 1.302 2020/02/22 21:00:05 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + 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 *** + + Revision 1.292 2019/05/09 14:17:20 brouard + Summary: Some updates + + Revision 1.291 2019/05/09 13:44:18 brouard + Summary: Before ncovmax + + Revision 1.290 2019/05/09 13:39:37 brouard + Summary: 0.99r18 unlimited number of individuals + + The number n which was limited to 20,000 cases is now unlimited, from firstobs to lastobs. If the number is too for the virtual memory, probably an error will occur. + + Revision 1.289 2018/12/13 09:16:26 brouard + Summary: Bug for young ages (<-30) will be in r17 + Revision 1.288 2018/05/02 20:58:27 brouard Summary: Some bugs fixed @@ -1046,11 +1102,11 @@ typedef struct { #define NINTERVMAX 8 #define NLSTATEMAX 8 /**< Maximum number of live states (for func) */ #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */ -#define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */ +#define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */ #define codtabm(h,k) (1 & (h-1) >> (k-1))+1 /*#define decodtabm(h,k,cptcoveff)= (h <= (1<> (k-1)) & 1) +1 : -1)*/ #define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1 -#define MAXN 20000 +/*#define MAXN 20000 */ /* Should by replaced by nobs, real number of observations and unlimited */ #define YEARM 12. /**< Number of months per year */ /* #define AGESUP 130 */ /* #define AGESUP 150 */ @@ -1070,12 +1126,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.288 2018/05/02 20:58:27 brouard Exp $ */ +/* $Id: imach.c,v 1.302 2020/02/22 21:00:05 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.288 $ $Date: 2018/05/02 20:58:27 $"; +char copyright[]="May 2019,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020"; +char fullversion[]="$Revision: 1.302 $ $Date: 2020/02/22 21:00:05 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -1098,8 +1154,9 @@ int nqfveff=0; /**< nqfveff Number of Qu int ntveff=0; /**< ntveff number of effective time varying variables */ int nqtveff=0; /**< ntqveff number of effective time varying quantitative variables */ 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 */ @@ -1238,7 +1295,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; @@ -2357,13 +2419,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 */ @@ -2951,8 +3013,8 @@ double **pmij(double **ps, double *cov, ps[ii][ii]=1; } } - - + + /* 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]); */ @@ -2972,7 +3034,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; @@ -2980,7 +3042,7 @@ double **pmij(double **ps, double *cov, double **out, **pmij(); double sumnew=0.; double agefin; - double k3=0.; /* constant of the w_x diagonal matrixe (in order for B to sum to 1 even for death state) */ + double k3=0.; /* constant of the w_x diagonal matrix (in order for B to sum to 1 even for death state) */ double **dnewm, **dsavm, **doldm; double **bbmij; @@ -2999,11 +3061,11 @@ double **pmij(double **ps, double *cov, /* outputs pmmij which is a stochastic matrix in row */ /* Diag(w_x) */ - /* Problem with prevacurrent which can be zero */ + /* Rescaling the cross-sectional prevalence: Problem with prevacurrent which can be zero */ 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 */ @@ -3026,10 +3088,10 @@ double **pmij(double **ps, double *cov, } /* End doldm, At the end doldm is diag[(w_i)] */ - /* left Product of this diag matrix by pmmij=Px (dnewm=dsavm*doldm) */ - bbmij=matprod2(dnewm, doldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, pmmij); /* Bug Valgrind */ + /* Left product of this diag matrix by pmmij=Px (dnewm=dsavm*doldm): diag[(w_i)*Px */ + bbmij=matprod2(dnewm, doldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, pmmij); /* was a Bug Valgrind */ - /* Diag(Sum_i w^i_x p^ij_x */ + /* Diag(Sum_i w^i_x p^ij_x, should be the prevalence at age x+stepm */ /* w1 p11 + w2 p21 only on live states N1./N..*N11/N1. + N2./N..*N21/N2.=(N11+N21)/N..=N.1/N.. */ for (j=1;j<=nlstate+ndeath;j++){ sumnew=0.; @@ -3047,7 +3109,7 @@ double **pmij(double **ps, double *cov, } /*End ii */ } /* End j, At the end dsavm is diag[1/(w_1p1i+w_2 p2i)] for ALL states even if the sum is only for live states */ - ps=matprod2(ps, dnewm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dsavm); /* Bug Valgrind */ + ps=matprod2(ps, dnewm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, dsavm); /* was a Bug Valgrind */ /* ps is now diag[w_i] * Px * diag [1/(w_1p1i+w_2 p2i)] */ /* end bmij */ return ps; /*pointer is unchanged */ @@ -3119,7 +3181,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++){ @@ -3879,7 +3941,7 @@ return -l; /*************** function likelione ***********/ -void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpri, long *ipmx, double *sw, double *fretone, double (*funcone)(double [])) +void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpri, long *ipmx, double *sw, double *fretone, double (*func)(double [])) { /* This routine should help understanding what is done with the selection of individuals/waves and @@ -3903,7 +3965,7 @@ void likelione(FILE *ficres,double p[], fprintf(ficresilk," -2*gipw/gsw*weight*ll(total)\n"); } - *fretone=(*funcone)(p); + *fretone=(*func)(p); if(*globpri !=0){ fclose(ficresilk); if (mle ==0) @@ -4378,6 +4440,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 ********************/ @@ -4952,6 +5028,7 @@ Title=%s
Datafile=%s Firstpass=%d La } } /* end mle=-2 */ dateintmean=dateintsum/k2cpt; + date2dmy(dateintmean,&jintmean,&mintmean,&aintmean); fclose(ficresp); fclose(ficresphtm); @@ -5145,11 +5222,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; @@ -5187,10 +5264,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; } @@ -5994,10 +6071,10 @@ 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 + /**< And for each alive state j, sums over i \f$ w^i_x {}{h}_p^{ij}x\f$, which are the probability * at horizon h in state j including mortality. */ for(j=1; j<= nlstate; j++){ @@ -6787,9 +6864,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,"\n", \ @@ -11710,10 +11891,10 @@ Interval (in months) between two waves: for(j=1;j<=NDIM;j++) ximort[i][j]=0.; /* ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */ - cens=ivector(1,n); - ageexmed=vector(1,n); - agecens=vector(1,n); - dcwave=ivector(1,n); + cens=ivector(firstobs,lastobs); + ageexmed=vector(firstobs,lastobs); + agecens=vector(firstobs,lastobs); + dcwave=ivector(firstobs,lastobs); for (i=1; i<=imx; i++){ dcwave[i]=-1; @@ -11747,8 +11928,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 @@ -11874,9 +12055,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++) { @@ -11923,13 +12104,14 @@ 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(cens,1,n); - free_vector(agecens,1,n); - free_ivector(dcwave,1,n); + free_ivector(dcwave,firstobs,lastobs); + free_vector(agecens,firstobs,lastobs); + free_vector(ageexmed,firstobs,lastobs); + free_ivector(cens,firstobs,lastobs); #ifdef GSL #endif } /* Endof if mle==-3 mortality only */ @@ -12125,6 +12307,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; @@ -12170,6 +12353,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; @@ -12195,6 +12379,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; @@ -12228,10 +12413,11 @@ Please run with mle=-1 to get a correct 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)) parameterline=13; @@ -12240,37 +12426,49 @@ Please run with mle=-1 to get a correct } switch (parameterline){ 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: @@ -12322,11 +12520,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); */ @@ -12335,8 +12566,8 @@ Please run with mle=-1 to get a correct /* free_imatrix(dh,1,lastpass-firstpass+2,1,imx); */ /* free_imatrix(bh,1,lastpass-firstpass+2,1,imx); */ /* free_imatrix(mw,1,lastpass-firstpass+2,1,imx); */ - free_lvector(num,1,n); - free_vector(agedc,1,n); + free_lvector(num,firstobs,lastobs); + free_vector(agedc,firstobs,lastobs); /*free_matrix(covar,0,NCOVMAX,1,n);*/ /*free_matrix(covar,1,NCOVMAX,1,n);*/ fclose(ficparo); @@ -12399,13 +12630,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); @@ -12420,8 +12661,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); @@ -12429,7 +12676,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------------ */ @@ -12650,12 +12897,12 @@ Please run with mle=-1 to get a correct varprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, prlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff); - free_vector(weight,1,n); + free_vector(weight,firstobs,lastobs); free_imatrix(Tvard,1,NCOVMAX,1,2); - free_imatrix(s,1,maxwav+1,1,n); - free_matrix(anint,1,maxwav,1,n); - free_matrix(mint,1,maxwav,1,n); - free_ivector(cod,1,n); + free_imatrix(s,1,maxwav+1,firstobs,lastobs); + free_matrix(anint,1,maxwav,firstobs,lastobs); + free_matrix(mint,1,maxwav,firstobs,lastobs); + free_ivector(cod,firstobs,lastobs); free_ivector(tab,1,NCOVMAX); fclose(ficresstdeij); fclose(ficrescveij); @@ -12675,10 +12922,10 @@ Please run with mle=-1 to get a correct free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath); free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath); - if(ntv+nqtv>=1)free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,1,n); - if(nqtv>=1)free_ma3x(cotqvar,1,maxwav,1,nqtv,1,n); - if(nqv>=1)free_matrix(coqvar,1,nqv,1,n); - free_matrix(covar,0,NCOVMAX,1,n); + if(ntv+nqtv>=1)free_ma3x(cotvar,1,maxwav,1,ntv+nqtv,firstobs,lastobs); + if(nqtv>=1)free_ma3x(cotqvar,1,maxwav,1,nqtv,firstobs,lastobs); + if(nqv>=1)free_matrix(coqvar,1,nqv,firstobs,lastobs); + free_matrix(covar,0,NCOVMAX,firstobs,lastobs); free_matrix(matcov,1,npar,1,npar); free_matrix(hess,1,npar,1,npar); /*free_vector(delti,1,npar);*/ @@ -12798,13 +13045,16 @@ Please run with mle=-1 to get a correct sprintf(plotcmd,"%s %s",pplotcmd, optionfilegnuplot); printf("Starting graphs with: '%s'\n",plotcmd);fflush(stdout); + strcpy(pplotcmd,plotcmd); if((outcmd=system(plotcmd)) != 0){ - printf("gnuplot command might not be in your path: '%s', err=%d\n", plotcmd, outcmd); + printf("Error in gnuplot, command might not be in your path: '%s', err=%d\n", plotcmd, outcmd); printf("\n Trying if gnuplot resides on the same directory that IMaCh\n"); sprintf(plotcmd,"%sgnuplot %s", pathimach, optionfilegnuplot); - if((outcmd=system(plotcmd)) != 0) + if((outcmd=system(plotcmd)) != 0){ printf("\n Still a problem with gnuplot command %s, err=%d\n", plotcmd, outcmd); + strcpy(plotcmd,pplotcmd); + } } printf(" Successful, please wait..."); while (z[0] != 'q') {