\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1])); + /* mu = \frac{w x}{\sum w} + var = \frac{\sum w (x-mu)^2}{\sum w} = \frac{w x^2}{\sum w} - mu^2 + */ + printf(" fixed quantitative variable V%d on %.3g (weighted) representatives of the population : %8.5g (%8.5g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt(stdq[z1]/idq[z1]-meanq[z1]*meanq[z1]/idq[z1]/idq[z1])); + fprintf(ficlog," fixed quantitative variable V%d on %.3g (weighted) representatives of the population : %8.5g (%8.5g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt(stdq[z1]/idq[z1]-meanq[z1]*meanq[z1]/idq[z1]/idq[z1])); + fprintf(ficresphtmfr," fixed quantitative variable V%d on %.3g (weighted) representatives of the population : %8.5g (%8.5g)
\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt(stdq[z1]/idq[z1]-meanq[z1]*meanq[z1]/idq[z1]/idq[z1])); } /* for (z1=1; z1<= nqtveff; z1++) { */ /* for(m=1;m<=lastpass;m++){ */ @@ -5189,11 +5273,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; @@ -5214,33 +5298,33 @@ void concatwav(int wav[], int **dh, int for(i=1; i<=imx; i++){ /* For simple cases and if state is death */ mi=0; /* First valid wave */ mli=0; /* Last valid wave */ - m=firstpass; - while(s[m][i] <= nlstate){ /* a live state */ + m=firstpass; /* Loop on waves */ + while(s[m][i] <= nlstate){ /* a live state or unknown state */ if(m >firstpass && s[m][i]==s[m-1][i] && mint[m][i]==mint[m-1][i] && anint[m][i]==anint[m-1][i]){/* Two succesive identical information on wave m */ mli=m-1;/* mw[++mi][i]=m-1; */ }else if(s[m][i]>=1 || s[m][i]==-4 || s[m][i]==-5){ /* Since 0.98r4 if status=-2 vital status is really unknown, wave should be skipped */ - mw[++mi][i]=m; + mw[++mi][i]=m; /* Valid wave: incrementing mi and updating mi; mw[mi] is the wave number of mi_th valid transition */ mli=m; } /* else might be a useless wave -1 and mi is not incremented and mw[mi] not updated */ if(m < lastpass){ /* m < lastpass, standard case */ m++; /* mi gives the "effective" current wave, m the current wave, go to next wave by incrementing m */ } - else{ /* m >= lastpass, eventual special issue with warning */ + else{ /* m = lastpass, eventual special issue with warning */ #ifdef UNKNOWNSTATUSNOTCONTRIBUTING break; #else - if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){ + if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){ /* case -2 (vital status unknown is warned later */ 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); - mw[++mi][i]=m; + fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood 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; /* Valid transition with unknown status */ mli=m; } if(s[m][i]==-2){ /* Vital status is really unknown */ nbwarn++; - if((int)anint[m][i] == 9999){ /* Has the vital status really been verified? */ + if((int)anint[m][i] == 9999){ /* Has the vital status really been verified?not a transition */ printf("Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); fprintf(ficlog,"Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); } @@ -5265,34 +5349,35 @@ void concatwav(int wav[], int **dh, int #ifndef DISPATCHINGKNOWNDEATHAFTERLASTWAVE else if ((int) andc[i] != 9999) { /* Date of death is known */ if ((int)anint[m][i]!= 9999) { /* date of last interview is known */ - if((andc[i]+moisdc[i]/12.) <=(anint[m][i]+mint[m][i]/12.)){ /* death occured before last wave and status should have been death instead of -1 */ + if((andc[i]+moisdc[i]/12.) <=(anint[m][i]+mint[m][i]/12.)){ /* month of death occured before last wave month and status should have been death instead of -1 */ nbwarn++; if(firstfiv==0){ - printf("Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d interviewed at %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); + printf("Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d, interviewed on %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); firstfiv=1; }else{ - fprintf(ficlog,"Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d interviewed at %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); + fprintf(ficlog,"Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d, interviewed on %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); } - }else{ /* Death occured afer last wave potential bias */ + s[m][i]=nlstate+1; /* Fixing the status as death. Be careful if multiple death states */ + }else{ /* Month of Death occured afer last wave month, potential bias */ nberr++; if(firstwo==0){ - printf("Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood. Please add a new fictive wave at the date of last vital status scan, with a dead status or alive but unknown state status (-1). See documentation\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); + printf("Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d with status %d. Potential bias if other individuals are still alive on this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood. Please add a new fictitious wave at the date of last vital status scan, with a dead status. See documentation\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); firstwo=1; } - fprintf(ficlog,"Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood. Please add a new fictive wave at the date of last vital status scan, with a dead status or alive but unknown state status (-1). See documentation\n\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); + fprintf(ficlog,"Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d with status %d. Potential bias if other individuals are still alive on this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood. Please add a new fictitious wave at the date of last vital status scan, with a dead status. See documentation\n\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); } }else{ /* if date of interview is unknown */ /* death is known but not confirmed by death status at any wave */ if(firstfour==0){ - printf("Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); + printf("Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d with status %d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); firstfour=1; } - fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); + fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d with status %d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); } } /* end if date of death is known */ #endif - wav[i]=mi; /* mi should be the last effective wave (or mli) */ - /* wav[i]=mw[mi][i]; */ + wav[i]=mi; /* mi should be the last effective wave (or mli), */ + /* wav[i]=mw[mi][i]; */ if(mi==0){ nbwarn++; if(first==0){ @@ -5427,6 +5512,8 @@ void concatwav(int wav[], int **dh, int if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ switch(Fixed[k]) { case 0: /* Testing on fixed dummy covariate, simple or product of fixed */ + modmaxcovj=0; + modmincovj=0; for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the modality of this covariate Vj*/ ij=(int)(covar[Tvar[k]][i]); /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i @@ -5440,8 +5527,10 @@ void concatwav(int wav[], int **dh, int else if (ij < modmincovj) modmincovj=ij; if (ij <0 || ij >1 ){ - printf("Information, IMaCh doesn't treat covariate with missing values (-1), individual %d will be skipped.\n",i); - fprintf(ficlog,"Information, currently IMaCh doesn't treat covariate with missing values (-1), individual %d will be skipped.\n",i); + printf("ERROR, IMaCh doesn't treat covariate with missing values V%d=-1, individual %d will be skipped.\n",Tvar[k],i); + fprintf(ficlog,"ERROR, currently IMaCh doesn't treat covariate with missing values V%d=-1, individual %d will be skipped.\n",Tvar[k],i); + fflush(ficlog); + exit(1); } if ((ij < -1) || (ij > NCOVMAX)){ printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX ); @@ -5516,6 +5605,16 @@ void concatwav(int wav[], int **dh, int break; } /* end switch */ } /* end dummy test */ + if(Dummy[k]==1 && Typevar[k] !=1){ /* Dummy covariate and not age product */ + for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the modality of this covariate Vj*/ + if(isnan(covar[Tvar[k]][i])){ + printf("ERROR, IMaCh doesn't treat fixed quantitative covariate with missing values V%d=., individual %d will be skipped.\n",Tvar[k],i); + fprintf(ficlog,"ERROR, currently IMaCh doesn't treat covariate with missing values V%d=., individual %d will be skipped.\n",Tvar[k],i); + fflush(ficlog); + exit(1); + } + } + } } /* end of loop on model-covariate k. nbcode[Tvark][1]=-1, nbcode[Tvark][1]=0 and nbcode[Tvark][2]=1 sets the value of covariate k*/ for (k=-1; k< maxncov; k++) Ndum[k]=0; @@ -8321,20 +8420,6 @@ set ter svg size 640, 480\nunset log y\n }/* End movingaverage */ - int 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; - return; - } /************** Forecasting ******************/ /* 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)*/ @@ -8490,7 +8575,7 @@ void prevforecast(char fileres[], double */ 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,jintmean,mintmean,aintmean; + double agelim, ppij, ppi, yp,yp1,yp2; /* ,jintmean,mintmean,aintmean;*/ double *popeffectif,*popcount; double ***p3mat; /* double ***mobaverage; */ @@ -9046,7 +9131,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++) { @@ -9054,28 +9139,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; } @@ -9084,7 +9175,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 */ @@ -9177,6 +9268,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; @@ -9210,8 +9302,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 */ @@ -9332,7 +9469,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 */ @@ -9371,7 +9512,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 */ @@ -9391,6 +9539,8 @@ int readdata(char datafile[], int firsto cutv(stra, strb, line, ' '); if(strb[0]=='.') { /* Missing value */ lval=-1; + coqvar[iv][i]=NAN; + covar[ncovcol+iv][i]=NAN; /* including qvar in standard covar for performance reasons */ }else{ errno=0; /* what_kind_of_number(strb); */ @@ -9494,7 +9644,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); @@ -9509,15 +9658,14 @@ 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); + printf("ERROR: the number of variables in the resultline, %d, differs from the number of variables used in the model line, %d.\n",j, cptcovs); + fprintf(ficlog,"ERROR: the number of variables in the resultline, %d, differs from the number of variables used in the model line, %d.\n",j, cptcovs); } for(k=1; k<=j;k++){ /* Loop on any covariate of the result line */ if(nbocc(resultsav,'=') >1){ cutl(stra,strb,resultsav,' '); /* keeps in strb after the first ' ' - resultsav= V4=1 V5=25.1 V3=0 strb=V3=0 stra= V4=1 V5=25.1 */ + resultsav= V4=1 V5=25.1 V3=0 stra= V5=25.1 V3=0 strb= V4=1 */ cutl(strc,strd,strb,'='); /* strb:V4=1 strc=1 strd=V4 */ }else cutl(strc,strd,resultsav,'='); @@ -9542,7 +9690,9 @@ int decoderesult ( char resultline[], in } } if(match == 0){ - printf("Error in result line: %d value missing; result: %s, model=%s\n",k1, resultline, model); + printf("Error in result line: V%d is missing in result: %s according to model=%s\n",k1, resultline, model); + fprintf(ficlog,"Error in result line: V%d is missing in result: %s according to model=%s\n",k1, resultline, model); + return 1; } } } @@ -9559,8 +9709,12 @@ int decoderesult ( char resultline[], in } if(match == 0){ printf("Error in result line: %d value missing; result: %s, model=%s\n",k1, resultline, model); + fprintf(ficlog,"Error in result line: %d value missing; result: %s, model=%s\n",k1, resultline, model); + return 1; }else if(match > 1){ printf("Error in result line: %d doubled; result: %s, model=%s\n",k2, resultline, model); + fprintf(ficlog,"Error in result line: %d doubled; result: %s, model=%s\n",k2, resultline, model); + return 1; } } @@ -10748,8 +10902,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*/ @@ -11070,7 +11224,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"); @@ -11158,8 +11313,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; } @@ -11845,8 +12000,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 @@ -11972,9 +12127,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++) { @@ -12021,9 +12176,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); @@ -12224,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; @@ -12269,6 +12425,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; @@ -12294,6 +12451,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; @@ -12316,28 +12474,32 @@ 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,"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 && (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); @@ -12350,7 +12512,9 @@ Please run with mle=-1 to get a correct prvforecast = 1; } else if((num_filled=sscanf(line,"prevforecast=%d yearsfproj=%lf mobil_average=%d\n",&prevfcast,&yrfproj,&mobilavproj)) !=EOF){/* && (num_filled == 3))*/ - printf(" Num_filled=%d, yearsfproj=%lf, mobil_average=%d\n",prevfcast,yrfproj,mobilavproj); + 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 { @@ -12371,7 +12535,9 @@ Please run with mle=-1 to get a correct prvbackcast = 1; } else if((num_filled=sscanf(line,"prevbackcast=%d yearsbproj=%lf mobil_average=%d\n",&prevbcast,&yrbproj,&mobilavproj)) ==3){/* && (num_filled == 3))*/ - printf(" Num_filled=%d, yearsbproj=%lf, mobil_average=%d\n",prevbcast,yrbproj,mobilavproj); + 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 { @@ -12380,55 +12546,38 @@ Please run with mle=-1 to get a correct goto end; } break; - /* /\*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.*\/ */ - /* dateback1=anback1+(mback1-1)/12.+(jback1-1)/365.; */ - /* dateback2=anback2+(mback2-1)/12.+(jback2-1)/365.; */ - /* } */ - /* 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; + } + if(!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); + } else + goto end; + 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 */