--- imach/src/imach.c 2019/05/09 13:44:18 1.291 +++ imach/src/imach.c 2021/03/08 18:11:32 1.307 @@ -1,6 +1,69 @@ -/* $Id: imach.c,v 1.291 2019/05/09 13:44:18 brouard Exp $ +/* $Id: imach.c,v 1.307 2021/03/08 18:11:32 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + 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 *** + + 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 @@ -1081,12 +1144,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.291 2019/05/09 13:44:18 brouard Exp $ */ +/* $Id: imach.c,v 1.307 2021/03/08 18:11:32 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.291 $ $Date: 2019/05/09 13:44:18 $"; +char copyright[]="Feb 2021,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, INED 2000-2021"; +char fullversion[]="$Revision: 1.307 $ $Date: 2021/03/08 18:11:32 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -1111,7 +1174,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 */ @@ -1250,7 +1313,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; @@ -2369,13 +2437,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 */ @@ -2963,8 +3031,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]); */ @@ -2984,7 +3052,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; @@ -2992,7 +3060,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; @@ -3011,11 +3079,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 */ @@ -3038,10 +3106,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.; @@ -3059,7 +3127,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 */ @@ -3131,7 +3199,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++){ @@ -3891,7 +3959,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 @@ -3915,7 +3983,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) @@ -4390,6 +4458,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 ********************/ @@ -4964,6 +5046,7 @@ Title=%s
Datafile=%s Firstpass=%d La } } /* end mle=-2 */ dateintmean=dateintsum/k2cpt; + date2dmy(dateintmean,&jintmean,&mintmean,&aintmean); fclose(ficresp); fclose(ficresphtm); @@ -5157,11 +5240,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; @@ -5199,10 +5282,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; } @@ -6006,10 +6089,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++){ @@ -6799,9 +6882,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,"