--- imach/src/imach.c 2006/06/30 13:02:05 1.128 +++ imach/src/imach.c 2009/10/29 13:18:53 1.134 @@ -1,6 +1,26 @@ -/* $Id: imach.c,v 1.128 2006/06/30 13:02:05 brouard Exp $ +/* $Id: imach.c,v 1.134 2009/10/29 13:18:53 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.134 2009/10/29 13:18:53 brouard + (Module): Now imach stops if date of birth, at least year of birth, is not given. Some cleaning of the code. + + Revision 1.133 2009/07/06 10:21:25 brouard + just nforces + + Revision 1.132 2009/07/06 08:22:05 brouard + Many tings + + Revision 1.131 2009/06/20 16:22:47 brouard + Some dimensions resccaled + + Revision 1.130 2009/05/26 06:44:34 brouard + (Module): Max Covariate is now set to 20 instead of 8. A + lot of cleaning with variables initialized to 0. Trying to make + V2+V3*age+V1+V4 strb=V3*age+V1+V4 working better. + + Revision 1.129 2007/08/31 13:49:27 lievre + Modification of the way of exiting when the covariate is not binary in order to see on the window the error message before exiting + Revision 1.128 2006/06/30 13:02:05 brouard (Module): Clarifications on computing e.j @@ -157,7 +177,7 @@ The same imach parameter file can be used but the option for mle should be -3. - Agnès, who wrote this part of the code, tried to keep most of the + Agnès, who wrote this part of the code, tried to keep most of the former routines in order to include the new code within the former code. The output is very simple: only an estimate of the intercept and of @@ -290,8 +310,8 @@ Also this programme outputs the covariance matrix of the parameters but also of the life expectancies. It also computes the period (stable) prevalence. - Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr). - Institut national d'études démographiques, Paris. + Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr). + Institut national d'études démographiques, Paris. This software have been partly granted by Euro-REVES, a concerted action from the European Union. It is copyrighted identically to a GNU software product, ie programme and @@ -364,13 +384,13 @@ extern int errno; #define GLOCK_ERROR_NOPATH -1 /* empty path */ #define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */ -#define MAXPARM 30 /* Maximum number of parameters for the optimization */ +#define MAXPARM 128 /* Maximum number of parameters for the optimization */ #define NPARMAX 64 /* (nlstate+ndeath-1)*nlstate*ncovmodel */ #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 8 /* Maximum number of covariates */ +#define NCOVMAX 20 /* Maximum number of covariates */ #define MAXN 20000 #define YEARM 12. /* Number of months per year */ #define AGESUP 130 @@ -386,41 +406,41 @@ extern int errno; #define ODIRSEPARATOR '/' #endif -/* $Id: imach.c,v 1.128 2006/06/30 13:02:05 brouard Exp $ */ +/* $Id: imach.c,v 1.134 2009/10/29 13:18:53 brouard Exp $ */ /* $State: Exp $ */ -char version[]="Imach version 0.98i, June 2006, INED-EUROREVES-Institut de longevite "; -char fullversion[]="$Revision: 1.128 $ $Date: 2006/06/30 13:02:05 $"; +char version[]="Imach version 0.98l, October 2009, INED-EUROREVES-Institut de longevite "; +char fullversion[]="$Revision: 1.134 $ $Date: 2009/10/29 13:18:53 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; -int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ -int nvar; -int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; +int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ +int nvar=0, nforce=0; /* Number of variables, number of forces */ +int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov=0; /* Number of covariates, of covariates with '*age' */ int npar=NPARMAX; int nlstate=2; /* Number of live states */ int ndeath=1; /* Number of dead states */ -int ncovmodel, ncovcol; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */ +int ncovmodel=0, ncovcol=0; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */ int popbased=0; int *wav; /* Number of waves for this individuual 0 is possible */ -int maxwav; /* Maxim number of waves */ -int jmin, jmax; /* min, max spacing between 2 waves */ -int ijmin, ijmax; /* Individuals having jmin and jmax */ -int gipmx, gsw; /* Global variables on the number of contributions +int maxwav=0; /* Maxim number of waves */ +int jmin=0, jmax=0; /* min, max spacing between 2 waves */ +int ijmin=0, ijmax=0; /* Individuals having jmin and jmax */ +int gipmx=0, gsw=0; /* Global variables on the number of contributions to the likelihood and the sum of weights (done by funcone)*/ -int mle, weightopt; +int mle=1, weightopt=0; int **mw; /* mw[mi][i] is number of the mi wave for this individual */ int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */ int **bh; /* bh[mi][i] is the bias (+ or -) for this individual if the delay between * wave mi and wave mi+1 is not an exact multiple of stepm. */ -double jmean; /* Mean space between 2 waves */ +double jmean=1; /* Mean space between 2 waves */ double **oldm, **newm, **savm; /* Working pointers to matrices */ double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ FILE *fic,*ficpar, *ficparo,*ficres, *ficresp, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop; FILE *ficlog, *ficrespow; -int globpr; /* Global variable for printing or not */ +int globpr=0; /* Global variable for printing or not */ double fretone; /* Only one call to likelihood */ -long ipmx; /* Number of contributions */ +long ipmx=0; /* Number of contributions */ double sw; /* Sum of weights */ char filerespow[FILENAMELENGTH]; char fileresilk[FILENAMELENGTH]; /* File of individual contributions to the likelihood */ @@ -513,7 +533,7 @@ double dateintmean=0; double *weight; int **s; /* Status */ double *agedc, **covar, idx; -int **nbcode, *Tcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff; +int **nbcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff; double *lsurv, *lpop, *tpop; double ftol=FTOL; /* Tolerance for computing Max Likelihood */ @@ -586,6 +606,20 @@ void replace_back_to_slash(char *s, char } } +char *trimbb(char *out, char *in) +{ /* Trim multiple blanks in line */ + char *s; + s=out; + while (*in != '\0'){ + while( *in == ' ' && *(in+1) == ' ' && *(in+1) != '\0'){ + in++; + } + *out++ = *in++; + } + *out='\0'; + return s; +} + int nbocc(char *s, char occ) { int i,j=0; @@ -1163,7 +1197,7 @@ double **prevalim(double **prlim, int nl int i, ii,j,k; double min, max, maxmin, maxmax,sumnew=0.; double **matprod2(); - double **out, cov[NCOVMAX], **pmij(); + double **out, cov[NCOVMAX+1], **pmij(); double **newm; double agefin, delaymax=50 ; /* Max number of years to converge */ @@ -1245,10 +1279,14 @@ double **pmij(double **ps, double *cov, for(i=1; i<= nlstate; i++){ s1=0; - for(j=1; j=1.e-10){ if(first==1){ - printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); + printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); } fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); }else{ @@ -2210,7 +2249,7 @@ void prevalence(double ***probs, double We still use firstpass and lastpass as another selection. */ - int i, m, jk, k1, i1, j1, bool, z1,z2,j; + int i, m, jk, k1, i1, j1, bool, z1,j; double ***freq; /* Frequencies */ double *pp, **prop; double pos,posprop; @@ -2427,35 +2466,41 @@ void concatwav(int wav[], int **dh, int void tricode(int *Tvar, int **nbcode, int imx) { - int Ndum[20],ij=1, k, j, i, maxncov=19; + /* Tvar[i]=atoi(stre); /* find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 */ + + int Ndum[20],ij=1, k=0, j=0, i=0, maxncov=NCOVMAX; int cptcode=0; cptcoveff=0; for (k=0; k cptcode) cptcode=ij; /* getting the maximum of covariable + if (ij > cptcode) cptcode=ij; /* getting the maximum value of the modality of the covariate (should be 0 or 1 now) Tvar[j]. If V=sex and male is 0 and female is 1, then cptcode=1.*/ } - for (i=0; i<=cptcode; i++) { - if(Ndum[i]!=0) ncodemax[j]++; /* Nomber of modalities of the j th covariates. In fact ncodemax[j]=2 (dichotom. variables) but it can be more */ - } + for (i=0; i<=cptcode; i++) { /* i=-1 ?*/ + if(Ndum[i]!=0) ncodemax[j]++; /* Nomber of modalities of the j + th covariate. In fact + ncodemax[j]=2 + (dichotom. variables only) but + it can be more */ + } /* Ndum[-1] number of undefined modalities */ ij=1; - for (i=1; i<=ncodemax[j]; i++) { - for (k=0; k<= maxncov; k++) { - if (Ndum[k] != 0) { - nbcode[Tvar[j]][ij]=k; - /* store the modality in an array. k is a modality. If we have model=V1+V1*sex then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */ - + for (i=1; i<=ncodemax[j]; i++) { /* i= 1 to 2 */ + for (k=0; k<= maxncov; k++) { /* k=-1 ?*/ + if (Ndum[k] != 0) { /* If at least one individual responded to this modality k */ + nbcode[Tvar[j]][ij]=k; /* stores the modality in an array nbcode. + k is a modality. If we have model=V1+V1*sex + then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */ ij++; } if (ij > ncodemax[j]) break; @@ -2465,9 +2510,9 @@ void tricode(int *Tvar, int **nbcode, in for (k=0; k< maxncov; k++) Ndum[k]=0; - for (i=1; i<=ncovmodel-2; i++) { + for (i=1; i<=ncovmodel-2; i++) { /* -2, cste and age */ /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ - ij=Tvar[i]; + ij=Tvar[i]; /* Tvar might be -1 if status was unknown */ Ndum[ij]++; } @@ -2478,8 +2523,8 @@ void tricode(int *Tvar, int **nbcode, in ij++; } } - - cptcoveff=ij-1; /*Number of simple covariates*/ + ij--; + cptcoveff=ij; /*Number of simple covariates*/ } /*********** Health Expectancies ****************/ @@ -3077,9 +3122,9 @@ void varevsij(char optionfilefiname[], d free_vector(gmp,nlstate+1,nlstate+ndeath); free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath); free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/ - fprintf(ficgp,"\nset noparametric;set nolabel; set ter png small;set size 0.65, 0.65"); + fprintf(ficgp,"\nunset parametric;unset label; set ter png small;set size 0.65, 0.65"); /* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */ - fprintf(ficgp,"\n set log y; set nolog x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";"); + fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";"); /* fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */ /* fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */ /* fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */ @@ -3256,7 +3301,7 @@ void varprob(char optionfilefiname[], do fprintf(ficresprobcov,"\n"); fprintf(ficresprobcor,"\n"); */ - xp=vector(1,npar); + xp=vector(1,npar); dnewm=matrix(1,(nlstate)*(nlstate+ndeath),1,npar); doldm=matrix(1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath)); mu=matrix(1,(nlstate)*(nlstate+ndeath), (int) bage, (int)fage); @@ -3409,7 +3454,7 @@ To be simple, these graphs help to under /* Confidence intervalle of pij */ /* - fprintf(ficgp,"\nset noparametric;unset label"); + fprintf(ficgp,"\nunset parametric;unset label"); fprintf(ficgp,"\nset log y;unset log x; set xlabel \"Age\";set ylabel \"probability (year-1)\""); fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65"); fprintf(fichtm,"\n
Probability with confidence intervals expressed in year-1 :pijgr%s.png, ",optionfilefiname,optionfilefiname); @@ -3644,8 +3689,8 @@ true period expectancies (those weighted void printinggnuplot(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){ char dirfileres[132],optfileres[132]; - int m,cpt,k1,i,k,j,jk,k2,k3,ij,l; - int ng; + int m0,cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0; + int ng=0; /* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */ /* printf("Problem with file %s",optionfilegnuplot); */ /* fprintf(ficlog,"Problem with file %s",optionfilegnuplot); */ @@ -3671,7 +3716,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u for (i=1; i<= nlstate ; i ++) { if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); - else fprintf(ficgp," \%%*lf (\%%*lf)"); + else fprintf(ficgp," \%%*lf (\%%*lf)"); } fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); for (i=1; i<= nlstate ; i ++) { @@ -4413,7 +4458,8 @@ int main(int argc, char *argv[]) double ***mobaverage; int *indx; char line[MAXLINE], linepar[MAXLINE]; - char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE]; + char linetmp[MAXLINE]; + char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE]; char pathr[MAXLINE], pathimach[MAXLINE]; char **bp, *tok, *val; /* pathtot */ int firstobs=1, lastobs=10; @@ -4617,11 +4663,18 @@ int main(int argc, char *argv[]) covar=matrix(0,NCOVMAX,1,n); cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement*/ if (strlen(model)>1) cptcovn=nbocc(model,'+')+1; - - ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age */ + /* where is ncovprod ?*/ + ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age : v1+v2+v3+v2*v4+v5*age makes 5+2=7*/ nvar=ncovmodel-1; /* Suppressing age as a basic covariate */ - npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/ - + nforce= (nlstate+ndeath-1)*nlstate; /* Number of forces ij from state i to j */ + npar= nforce*ncovmodel; /* Number of parameters like aij*/ + if(npar >MAXPARM || nlstate >NLSTATEMAX || ndeath >NDEATHMAX || ncovmodel>NCOVMAX){ + printf("Too complex model for current IMaCh: npar=(nlstate+ndeath-1)*nlstate*ncovmodel=%d >= %d(MAXPARM) or nlstate=%d >= %d(NLSTATEMAX) or ndeath=%d >= %d(NDEATHMAX) or ncovmodel=(k+age+#of+signs)=%d(NCOVMAX) >= %d\n",npar, MAXPARM, nlstate, NLSTATEMAX, ndeath, NDEATHMAX, ncovmodel, NCOVMAX); + fprintf(ficlog,"Too complex model for current IMaCh: %d >=%d(MAXPARM) or %d >=%d(NLSTATEMAX) or %d >=%d(NDEATHMAX) or %d(NCOVMAX) >=%d\n",npar, MAXPARM, nlstate, NLSTATEMAX, ndeath, NDEATHMAX, ncovmodel, NCOVMAX); + fflush(stdout); + fclose (ficlog); + goto end; + } delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); delti=delti3[1][1]; /*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/ @@ -4747,6 +4800,9 @@ run imach with mle=-1 to get a correct t ungetc(c,ficpar); matcov=matrix(1,npar,1,npar); + for(i=1; i <=npar; i++) + for(j=1; j <=npar; j++) matcov[i][j]=0.; + for(i=1; i <=npar; i++){ fscanf(ficpar,"%s",&str); if(mle==1) @@ -4810,7 +4866,7 @@ run imach with mle=-1 to get a correct t for(i=1;i<=n;i++) weight[i]=1.0; /* Equal weights, 1 by default */ mint=matrix(1,maxwav,1,n); anint=matrix(1,maxwav,1,n); - s=imatrix(1,maxwav+1,1,n); + s=imatrix(1,maxwav+1,1,n); /* s[i][j] health state for wave i and individual j */ tab=ivector(1,NCOVMAX); ncodemax=ivector(1,8); @@ -4831,15 +4887,25 @@ run imach with mle=-1 to get a correct t printf("Comment line\n%s\n",line); continue; } + trimbb(linetmp,line); /* Trims multiple blanks in line */ + for (j=0; line[j]!='\0';j++){ + line[j]=linetmp[j]; + } + for (j=maxwav;j>=1;j--){ cutv(stra, strb,line,' '); - errno=0; - lval=strtol(strb,&endptr,10); + if(strb[0]=='.') { /* Missing status */ + lval=-1; + }else{ + errno=0; + lval=strtol(strb,&endptr,10); /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/ - if( strb[0]=='\0' || (*endptr != '\0')){ - printf("Error reading data around '%d' at line number %d %s for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav); - exit(1); + if( strb[0]=='\0' || (*endptr != '\0')){ + printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav); + fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog); + goto end; + } } s[j][i]=lval; @@ -4851,8 +4917,9 @@ run imach with mle=-1 to get a correct t month=99; year=9999; }else{ - printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j); - exit(1); + printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j); + fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j);fflush(ficlog); + goto end; } anint[j][i]= (double) year; mint[j][i]= (double)month; @@ -4866,8 +4933,9 @@ run imach with mle=-1 to get a correct t month=99; year=9999; }else{ - printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line); - exit(1); + printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line); + fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog); + goto end; } andc[i]=(double) year; moisdc[i]=(double) month; @@ -4880,8 +4948,15 @@ run imach with mle=-1 to get a correct t month=99; year=9999; }else{ - printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line,j); - exit(1); + printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line); + fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog); + goto end; + } + if (year==9999) { + printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given. Exiting.\n",strb, linei,i,line); + fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given. Exiting.\n",strb, linei,i,line);fflush(ficlog); + goto end; + } annais[i]=(double)(year); moisnais[i]=(double)(month); @@ -4892,18 +4967,25 @@ run imach with mle=-1 to get a correct t dval=strtod(strb,&endptr); if( strb[0]=='\0' || (*endptr != '\0')){ printf("Error reading data around '%f' at line number %ld, \"%s\" for individual %d\nShould be a weight. Exiting.\n",dval, i,line,linei); - exit(1); + fprintf(ficlog,"Error reading data around '%f' at line number %ld, \"%s\" for individual %d\nShould be a weight. Exiting.\n",dval, i,line,linei); + fflush(ficlog); + goto end; } weight[i]=dval; strcpy(line,stra); for (j=ncovcol;j>=1;j--){ cutv(stra, strb,line,' '); - errno=0; - lval=strtol(strb,&endptr,10); - if( strb[0]=='\0' || (*endptr != '\0')){ - printf("Error reading data around '%d' at line number %ld %s for individual %d, '%s'\nShould be a covar (meaning 0 for the reference or 1). Exiting.\n",lval, linei,i, line); - exit(1); + if(strb[0]=='.') { /* Missing status */ + lval=-1; + }else{ + errno=0; + lval=strtol(strb,&endptr,10); + if( strb[0]=='\0' || (*endptr != '\0')){ + printf("Error reading data around '%d' at line number %ld for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative). Exiting.\n",lval, linei,i, line); + fprintf(ficlog,"Error reading data around '%d' at line number %ld for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative). Exiting.\n",lval, linei,i, line);fflush(ficlog); + goto end; + } } if(lval <-1 || lval >1){ printf("Error reading data around '%d' at line number %ld for individual %d, '%s'\n \ @@ -4915,13 +4997,22 @@ run imach with mle=-1 to get a correct t and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \ output of IMaCh is often meaningless.\n \ Exiting.\n",lval,linei, i,line,j); - exit(1); + fprintf(ficlog,"Error reading data around '%d' at line number %ld for individual %d, '%s'\n \ + Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \ + for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \ + For example, for multinomial values like 1, 2 and 3,\n \ + build V1=0 V2=0 for the reference value (1),\n \ + V1=1 V2=0 for (2) \n \ + and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \ + output of IMaCh is often meaningless.\n \ + Exiting.\n",lval,linei, i,line,j);fflush(ficlog); + goto end; } covar[j][i]=(double)(lval); strcpy(line,stra); - } + } lstra=strlen(stra); - + if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */ stratrunc = &(stra[lstra-9]); num[i]=atol(stratrunc); @@ -4953,7 +5044,7 @@ run imach with mle=-1 to get a correct t else weight[i]=1;*/ /* Calculation of the number of parameters from char model */ - Tvar=ivector(1,15); /* stores the number n of the covariates in Vm+Vn at 1 and m at 2 */ + Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. Stores the number n of the covariates in Vm+Vn at 1 and m at 2 */ Tprod=ivector(1,15); Tvaraff=ivector(1,15); Tvard=imatrix(1,15,1,2); @@ -4963,32 +5054,35 @@ run imach with mle=-1 to get a correct t j=0, j1=0, k1=1, k2=1; j=nbocc(model,'+'); /* j=Number of '+' */ j1=nbocc(model,'*'); /* j1=Number of '*' */ - cptcovn=j+1; - cptcovprod=j1; /*Number of products */ + cptcovn=j+1; /* Number of covariates V1+V2+V3 =>2+1=3 */ + cptcovprod=j1; /*Number of products V1*V2 =1 */ strcpy(modelsav,model); if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){ printf("Error. Non available option model=%s ",model); - fprintf(ficlog,"Error. Non available option model=%s ",model); + fprintf(ficlog,"Error. Non available option model=%s ",model);fflush(ficlog); goto end; } /* This loop fills the array Tvar from the string 'model'.*/ - + /* j is the number of + signs in the model V1+V2+V3 j=2 i=3 to 1 */ for(i=(j+1); i>=1;i--){ - cutv(stra,strb,modelsav,'+'); /* keeps in strb after the last + */ + cutv(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' + modelsav=V2+V3*age+V1+V4 strb=V3*age+V1+V4 + stra=V2 + */ if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */ /* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/ /*scanf("%d",i);*/ - if (strchr(strb,'*')) { /* Model includes a product */ - cutv(strd,strc,strb,'*'); /* strd*strc Vm*Vn (if not *age)*/ + if (strchr(strb,'*')) { /* Model includes a product V1+V3*age+V2 strb=V3*age*/ + cutv(strd,strc,strb,'*'); /* strd*strc Vm*Vn: V3*age strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */ if (strcmp(strc,"age")==0) { /* Vn*age */ cptcovprod--; cutv(strb,stre,strd,'V'); - Tvar[i]=atoi(stre); /* computes n in Vn and stores in Tvar*/ - cptcovage++; - Tage[cptcovage]=i; - /*printf("stre=%s ", stre);*/ + Tvar[i]=atoi(stre); /* V1+V3*age+V2 Tvar[2]=3 */ + cptcovage++; /* Sums the number of covariates including age as a product */ + Tage[cptcovage]=i; /* Tage[1] =2 */ + /*printf("stre=%s ", stre);*/ } else if (strcmp(strd,"age")==0) { /* or age*Vn */ cptcovprod--; @@ -4997,11 +5091,12 @@ run imach with mle=-1 to get a correct t cptcovage++; Tage[cptcovage]=i; } - else { /* Age is not in the model */ - cutv(strb,stre,strc,'V'); /* strc= Vn, stre is n*/ - Tvar[i]=ncovcol+k1; + else { /* Age is not in the model V1+V3*V2+V2 strb=V3*V2*/ + cutv(strb,stre,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/ + Tvar[i]=ncovcol+k1; /* find 'n' in Vn and stores in Tvar. + If already ncovcol=2 and model=V2*V1 Tvar[1]=2+1 and Tvar[2]=2+2 etc */ cutv(strb,strc,strd,'V'); /* strd was Vm, strc is m */ - Tprod[k1]=i; + Tprod[k1]=i; /* Tprod[1] */ Tvard[k1][1]=atoi(strc); /* m*/ Tvard[k1][2]=atoi(stre); /* n */ Tvar[cptcovn+k2]=Tvard[k1][1]; @@ -5018,7 +5113,7 @@ run imach with mle=-1 to get a correct t cutv(strd,strc,strb,'V'); Tvar[i]=atoi(strc); } - strcpy(modelsav,stra); + strcpy(modelsav,stra); /* modelsav=V2+V3*age+V1+V4 strb=V3*age+V1+V4 */ /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav); scanf("%d",i);*/ } /* end of loop + */ @@ -5086,7 +5181,7 @@ run imach with mle=-1 to get a correct t agev[m][i]=1; else if(agev[m][i] agemax){ agemax=agev[m][i]; @@ -5148,7 +5243,6 @@ run imach with mle=-1 to get a correct t /* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */ - Tcode=ivector(1,100); nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); ncodemax[1]=1; if (cptcovn > 0) tricode(Tvar,nbcode,imx); @@ -5158,13 +5252,14 @@ run imach with mle=-1 to get a correct t h=0; m=pow(2,cptcoveff); - for(k=1;k<=cptcoveff; k++){ - for(i=1; i <=(m/pow(2,k));i++){ - for(j=1; j <= ncodemax[k]; j++){ - for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){ + for(k=1;k<=cptcoveff; k++){ /* scans any effective covariate */ + for(i=1; i <=(m/pow(2,k));i++){ /* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 */ + for(j=1; j <= ncodemax[k]; j++){ /* For each modality of this covariate */ + for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){ /* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 */ h++; - if (h>m) h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j; - /* printf("h=%d k=%d j=%d codtab[h][k]=%d tvar[k]=%d \n",h, k,j,codtab[h][k],Tvar[k]);*/ + if (h>m) + h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j; + printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][Tvar[k]]); } } } @@ -5173,7 +5268,7 @@ run imach with mle=-1 to get a correct t codtab[1][2]=1;codtab[2][2]=2; */ /* for(i=1; i <=m ;i++){ for(k=1; k <=cptcovn; k++){ - printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff); + printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff); } printf("\n"); } @@ -5201,7 +5296,8 @@ run imach with mle=-1 to get a correct t strcat(optionfilehtm,"-mort"); strcat(optionfilehtm,".htm"); if((fichtm=fopen(optionfilehtm,"w"))==NULL) { - printf("Problem with %s \n",optionfilehtm), exit(0); + printf("Problem with %s \n",optionfilehtm); + exit(0); } strcpy(optionfilehtmcov,optionfilefiname); /* Only for matrix of covariance */ @@ -5380,7 +5476,7 @@ Interval (in months) between two waves: } /* Endof if mle==-3 */ else{ /* For mle >=1 */ - + globpr=0;/* debug */ likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); for (k=1; k<=npar;k++) @@ -5653,7 +5749,8 @@ Interval (in months) between two waves: for(cptcov=1,k=0;cptcov<=i1;cptcov++){ for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ k=k+1; - /*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/ + /* to clean */ + printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,codtab[cptcod][cptcov],nbcode); fprintf(ficrespl,"\n#******"); printf("\n#******"); fprintf(ficlog,"\n#******"); @@ -5965,6 +6062,7 @@ Interval (in months) between two waves: free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX); } /* mle==-3 arrives here for freeing */ + endfree: free_matrix(prlim,1,nlstate,1,nlstate); free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); @@ -5982,7 +6080,6 @@ Interval (in months) between two waves: free_ivector(Tprod,1,15); free_ivector(Tvaraff,1,15); free_ivector(Tage,1,15); - free_ivector(Tcode,1,100); free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX); free_imatrix(codtab,1,100,1,10);