]> henry.ined.fr Git - .git/commitdiff
Some dimensions resccaled
authorN. Brouard <brouard@ined.fr>
Sat, 20 Jun 2009 16:22:47 +0000 (16:22 +0000)
committerN. Brouard <brouard@ined.fr>
Sat, 20 Jun 2009 16:22:47 +0000 (16:22 +0000)
src/imach.c

index 737677b9837168422972d9c73cf681fc03242a23..07d193bc6011f7d9ed8ff8dee4e9db4783cb69d5 100644 (file)
@@ -1,6 +1,11 @@
 /* $Id$
   $State$
   $Log$
+  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
 
@@ -367,7 +372,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
@@ -392,7 +397,7 @@ extern int errno;
 /* $Id$ */
 /* $State$ */
 
-char version[]="Imach version 0.98i, June 2006, INED-EUROREVES-Institut de longevite ";
+char version[]="Imach version 0.98k, June 2006, INED-EUROREVES-Institut de longevite ";
 char fullversion[]="$Revision$ $Date$"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
@@ -516,7 +521,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 */
@@ -1166,7 +1171,7 @@ double **prevalim(double **prlim, int nlstate, double x[], double age, double **
   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 */
 
@@ -1248,10 +1253,14 @@ double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
     
     for(i=1; i<= nlstate; i++){
       s1=0;
-      for(j=1; j<i; j++)
+      for(j=1; j<i; j++){
        s1+=exp(ps[i][j]);
-      for(j=i+1; j<=nlstate+ndeath; j++)
+       /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+      }
+      for(j=i+1; j<=nlstate+ndeath; j++){
        s1+=exp(ps[i][j]);
+       /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+      }
       ps[i][i]=1./(s1+1.);
       for(j=1; j<i; j++)
        ps[i][j]= exp(ps[i][j])*ps[i][i];
@@ -1317,7 +1326,7 @@ double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, in
      */
 
   int i, j, d, h, k;
-  double **out, cov[NCOVMAX];
+  double **out, cov[NCOVMAX+1];
   double **newm;
 
   /* Hstepm could be zero and should return the unit matrix */
@@ -1333,7 +1342,8 @@ double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, in
       /* Covariates have to be included here again */
       cov[1]=1.;
       cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM;
-      for (k=1; k<=cptcovn;k++) cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
+      for (k=1; k<=cptcovn;k++) 
+       cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
       for (k=1; k<=cptcovage;k++)
        cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
       for (k=1; k<=cptcovprod;k++)
@@ -1363,7 +1373,7 @@ double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, in
 double func( double *x)
 {
   int i, ii, j, k, mi, d, kk;
-  double l, ll[NLSTATEMAX], cov[NCOVMAX];
+  double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
   double **out;
   double sw; /* Sum of weights */
   double lli; /* Individual log likelihood */
@@ -1617,7 +1627,7 @@ double funcone( double *x)
 {
   /* Same as likeli but slower because of a lot of printf and if */
   int i, ii, j, k, mi, d, kk;
-  double l, ll[NLSTATEMAX], cov[NCOVMAX];
+  double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
   double **out;
   double lli; /* Individual log likelihood */
   double llt;
@@ -1679,7 +1689,7 @@ double funcone( double *x)
       ipmx +=1;
       sw += weight[i];
       ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
-/*       printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
+      printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); 
       if(globpr){
        fprintf(ficresilk,"%9d %6d %2d %2d %1d %1d %3d %11.6f %8.4f\
  %11.6f %11.6f %11.6f ", \
@@ -2432,35 +2442,35 @@ void tricode(int *Tvar, int **nbcode, int imx)
   
   /*     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<maxncov; k++) Ndum[k]=0;
-  for (k=1; k<=7; k++) ncodemax[k]=0;
+  for (k=1; k<=7; k++) ncodemax[k]=0; /* Horrible constant again */
 
-  for (j=1; j<=(cptcovn+2*cptcovprod); j++) {
+  for (j=1; j<=(cptcovn+2*cptcovprod); j++) { /* For each covariate */
     for (i=1; i<=imx; i++) { /*reads the data file to get the maximum 
                               modality*/ 
-      ij=(int)(covar[Tvar[j]][i]); /* ij is the modality of this individual*/
-      Ndum[ij]++; /*store the modality */
+      ij=(int)(covar[Tvar[j]][i]); /* ij is the modality of this individual, might be -1*/
+      Ndum[ij]++; /*counts the occurence of this modality */
       /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/
-      if (ij > 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; 
@@ -2470,9 +2480,9 @@ void tricode(int *Tvar, int **nbcode, int imx)
 
  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]++;
  }
 
@@ -2483,8 +2493,8 @@ void tricode(int *Tvar, int **nbcode, int imx)
      ij++;
    }
  }
- cptcoveff=ij-1; /*Number of simple covariates*/
+ ij--;
+ cptcoveff=ij; /*Number of simple covariates*/
 }
 
 /*********** Health Expectancies ****************/
@@ -3082,9 +3092,9 @@ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double
   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); */
@@ -3261,7 +3271,7 @@ void varprob(char optionfilefiname[], double **matcov, double x[], double delti[
   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);
@@ -3414,7 +3424,7 @@ To be simple, these graphs help to understand the significativity of each parame
 
       /* 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<br>Probability with  confidence intervals expressed in year<sup>-1</sup> :<a href=\"pijgr%s.png\">pijgr%s.png</A>, ",optionfilefiname,optionfilefiname);
@@ -3676,7 +3686,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,subdirf2(fil
 
      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 ++) {
@@ -4626,7 +4636,13 @@ int main(int argc, char *argv[])
   ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age */
   nvar=ncovmodel-1; /* Suppressing age as a basic covariate */
   npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/
-
+  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)*/
@@ -4752,6 +4768,9 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin
     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)
@@ -4815,7 +4834,7 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin
   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);
 
@@ -4839,12 +4858,17 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin
 
     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;
       
@@ -4856,8 +4880,9 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin
        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; 
@@ -4871,8 +4896,9 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin
       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; 
@@ -4886,7 +4912,8 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin
       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);
+      fprintf(ficlog,"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);fflush(ficlog);
+       goto end;
     }
     annais[i]=(double)(year);
     moisnais[i]=(double)(month); 
@@ -4897,18 +4924,25 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin
     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 \
@@ -4920,6 +4954,15 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin
  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);
@@ -4958,7 +5001,7 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin
      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);
@@ -4974,7 +5017,7 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin
     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;
     }
     
@@ -5157,7 +5200,6 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin
 
   /* 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);
@@ -5167,12 +5209,13 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin
   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]]);
        } 
       }
@@ -5182,7 +5225,7 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin
      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");
      }
@@ -5210,7 +5253,8 @@ run imach with mle=-1 to get a correct template of the parameter file.\n",numlin
     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 */
@@ -5389,8 +5433,8 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n",\
   } /* Endof if mle==-3 */
   
   else{ /* For mle >=1 */
-  
-    likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */
+    globpr=1;/* 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++)
       printf(" %d %8.5f",k,p[k]);
@@ -5662,7 +5706,8 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n",\
     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#******");
@@ -5974,6 +6019,7 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n",\
     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);
@@ -5991,7 +6037,6 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n",\
     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);