--- imach/src/imach.c 2009/05/26 06:44:34 1.130 +++ imach/src/imach.c 2009/10/29 13:18:53 1.134 @@ -1,6 +1,18 @@ -/* $Id: imach.c,v 1.130 2009/05/26 06:44:34 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 @@ -165,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 @@ -298,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 @@ -372,7 +384,7 @@ 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 @@ -394,15 +406,15 @@ extern int errno; #define ODIRSEPARATOR '/' #endif -/* $Id: imach.c,v 1.130 2009/05/26 06:44:34 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.130 $ $Date: 2009/05/26 06:44:34 $"; +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=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ -int nvar=0; +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 */ @@ -521,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 */ @@ -594,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; @@ -1171,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 */ @@ -1253,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{ @@ -2437,35 +2468,39 @@ void tricode(int *Tvar, int **nbcode, in /* 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=19; + 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; @@ -2475,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]++; } @@ -2488,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 ****************/ @@ -3087,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); */ @@ -3266,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); @@ -3419,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); @@ -3681,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 ++) { @@ -4423,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; @@ -4627,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)*/ @@ -4757,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) @@ -4820,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); @@ -4841,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; @@ -4861,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; @@ -4876,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; @@ -4890,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); @@ -4902,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 \ @@ -4925,6 +4997,15 @@ 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); + 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); @@ -4963,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); @@ -4979,7 +5060,7 @@ run imach with mle=-1 to get a correct t 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; } @@ -4999,7 +5080,7 @@ run imach with mle=-1 to get a correct t cptcovprod--; cutv(strb,stre,strd,'V'); Tvar[i]=atoi(stre); /* V1+V3*age+V2 Tvar[2]=3 */ - cptcovage++; /* Sum the number of covariates including ages as a product */ + cptcovage++; /* Sums the number of covariates including age as a product */ Tage[cptcovage]=i; /* Tage[1] =2 */ /*printf("stre=%s ", stre);*/ } @@ -5100,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]; @@ -5162,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); @@ -5172,12 +5252,13 @@ 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; + 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]]); } } @@ -5187,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"); } @@ -5215,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 */ @@ -5394,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++) @@ -5667,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#******"); @@ -5979,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); @@ -5996,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);