--- imach/src/imach.c	2015/07/16 16:49:02	1.192
+++ imach/src/imach.c	2015/09/03 07:14:39	1.198
@@ -1,6 +1,29 @@
-/* $Id: imach.c,v 1.192 2015/07/16 16:49:02 brouard Exp $
+/* $Id: imach.c,v 1.198 2015/09/03 07:14:39 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.198  2015/09/03 07:14:39  brouard
+  Summary: 0.98q5 Flavia
+
+  Revision 1.197  2015/09/01 18:24:39  brouard
+  *** empty log message ***
+
+  Revision 1.196  2015/08/18 23:17:52  brouard
+  Summary: 0.98q5
+
+  Revision 1.195  2015/08/18 16:28:39  brouard
+  Summary: Adding a hack for testing purpose
+
+  After reading the title, ftol and model lines, if the comment line has
+  a q, starting with #q, the answer at the end of the run is quit. It
+  permits to run test files in batch with ctest. The former workaround was
+  $ echo q | imach foo.imach
+
+  Revision 1.194  2015/08/18 13:32:00  brouard
+  Summary:  Adding error when the covariance matrix doesn't contain the exact number of lines required by the model line.
+
+  Revision 1.193  2015/08/04 07:17:42  brouard
+  Summary: 0.98q4
+
   Revision 1.192  2015/07/16 16:49:02  brouard
   Summary: Fixing some outputs
 
@@ -677,11 +700,12 @@ typedef struct {
 #define NLSTATEMAX 8 /**< Maximum number of live states (for func) */
 #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */
 #define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */
-#define codtabm(h,k)  1 & (h-1) >> (k-1) ;
+#define codtabm(h,k)  (1 & (h-1) >> (k-1))+1
 #define MAXN 20000
 #define YEARM 12. /**< Number of months per year */
 #define AGESUP 130
 #define AGEBASE 40
+#define AGEOVERFLOW 1.e20
 #define AGEGOMP 10 /**< Minimal age for Gompertz adjustment */
 #ifdef _WIN32
 #define DIRSEPARATOR '\\'
@@ -693,11 +717,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.192 2015/07/16 16:49:02 brouard Exp $ */
+/* $Id: imach.c,v 1.198 2015/09/03 07:14:39 brouard Exp $ */
 /* $State: Exp $ */
-
-char version[]="Imach version 0.98q3, July 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
-char fullversion[]="$Revision: 1.192 $ $Date: 2015/07/16 16:49:02 $"; 
+#include "version.h"
+char version[]=__IMACH_VERSION__;
+char copyright[]="September 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
+char fullversion[]="$Revision: 1.198 $ $Date: 2015/09/03 07:14:39 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -831,7 +856,7 @@ int estepm;
 
 int m,nb;
 long *num;
-int firstpass=0, lastpass=4,*cod, *Tage,*cens;
+int firstpass=0, lastpass=4,*cod, *cens;
 int *ncodemax;  /* ncodemax[j]= Number of modalities of the j th
 		   covariate for which somebody answered excluding 
 		   undefined. Usually 2: 0 and 1. */
@@ -851,6 +876,7 @@ double  **covar; /**< covar[j,i], value
 		  * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*age; */
 double  idx; 
 int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
+int *Tage;
 int *Ndum; /** Freq of modality (tricode */
 int **codtab; /**< codtab=imatrix(1,100,1,10); */
 int **Tvard, *Tprod, cptcovprod, *Tvaraff;
@@ -1888,13 +1914,13 @@ double **prevalim(double **prlim, int nl
     if(nagesqr==1)
       cov[3]= agefin*agefin;;
     for (k=1; k<=cptcovn;k++) {
-      cov[2+nagesqr+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
-      /*printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtab[%d][Tvar[%d]]=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], ij, k, codtab[ij][Tvar[k]]);*/
+      cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])];
+      /* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); */
     }
     /*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
-    for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]*cov[2];
+    for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]*cov[2];
     for (k=1; k<=cptcovprod;k++) /* Useless */
-      cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]] * nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
+      cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])];
     
     /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
     /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
@@ -2067,12 +2093,12 @@ double ***hpxij(double ***po, int nhstep
       if(nagesqr==1)
 	cov[3]= agexact*agexact;
       for (k=1; k<=cptcovn;k++) 
-	cov[2+nagesqr+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
+	cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])];
       for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */
 	/* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
-	cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2];
+	cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2];
       for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */
-	cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
+	cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])];
 
 
       /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
@@ -2891,7 +2917,7 @@ void lubksb(double **a, int n, int *indx
 
 void pstamp(FILE *fichier)
 {
-  fprintf(fichier,"# %s.%s\n#%s\n#%s\n# %s", optionfilefiname,optionfilext,version,fullversion,strstart);
+  fprintf(fichier,"# %s.%s\n#IMaCh version %s, %s\n#%s\n# %s", optionfilefiname,optionfilext,version,copyright, fullversion, strstart);
 }
 
 /************ Frequencies ********************/
@@ -2943,13 +2969,13 @@ void  freqsummary(char fileres[], int ia
 	bool=1;
 	if  (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
 	  for (z1=1; z1<=cptcoveff; z1++)       
-            if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]){
+            if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){
                 /* Tests if the value of each of the covariates of i is equal to filter j1 */
               bool=0;
-              /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtab[%d][%d]=%d, nbcode[Tvaraff][codtab[%d][%d]=%d, j1=%d\n", 
-                bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtab[j1][z1],
-                j1,z1,nbcode[Tvaraff[z1]][codtab[j1][z1]],j1);*/
-              /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtab[7][3]=1 and nbcde[3][?]=1*/
+              /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", 
+                bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),
+                j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/
+              /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/
             } 
 	}
  
@@ -2978,10 +3004,10 @@ void  freqsummary(char fileres[], int ia
       pstamp(ficresp);
       if  (cptcovn>0) {
 	fprintf(ficresp, "\n#********** Variable "); 
-	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
 	fprintf(ficresp, "**********\n#");
 	fprintf(ficlog, "\n#********** Variable "); 
-	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
 	fprintf(ficlog, "**********\n#");
       }
       for(i=1; i<=nlstate;i++) 
@@ -3109,7 +3135,7 @@ void prevalence(double ***probs, double
 	bool=1;
 	if  (cptcovn>0) {
 	  for (z1=1; z1<=cptcoveff; z1++) 
-	    if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) 
+	    if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) 
 	      bool=0;
 	} 
 	if (bool==1) { 
@@ -3375,14 +3401,15 @@ void tricode(int *Tvar, int **nbcode, in
        nbcode[Tvar[j]][1]=0;
        nbcode[Tvar[j]][2]=1;
        nbcode[Tvar[j]][3]=2;
+       To be continued (not working yet).
     */
-    ij=0; /* ij is similar to i but can jumps over null modalities */
-    for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 to 1*/
-    	if (Ndum[i] == 0) { /* If at least one individual responded to this modality k */
+    ij=0; /* ij is similar to i but can jump over null modalities */
+    for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/
+    	if (Ndum[i] == 0) { /* If nobody responded to this modality k */
 	  break;
 	}
 	ij++;
-	nbcode[Tvar[j]][ij]=i;  /* stores the original modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/
+	nbcode[Tvar[j]][ij]=i;  /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/
 	cptcode = ij; /* New max modality for covar j */
     } /* end of loop on modality i=-1 to 1 or more */
       
@@ -4211,10 +4238,9 @@ void varprob(char optionfilefiname[], do
   fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n");
   fprintf(fichtm,"\n");
 
-  fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of pairs of step probabilities (drawings)</a></h4></li>\n",optionfilehtmcov);
-  fprintf(fichtmcov,"\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n\
-  file %s<br>\n",optionfilehtmcov);
-  fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated\
+  fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4></br>this page is important in order to visualize confidence intervals and especially correlation between disability and recovery</li>\n",optionfilehtmcov);
+  fprintf(fichtmcov,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n",optionfilehtmcov, optionfilehtmcov);
+  fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated \
 and drawn. It helps understanding how is the covariance between two incidences.\
  They are expressed in year<sup>-1</sup> in order to be less dependent of stepm.<br>\n");
   fprintf(fichtmcov,"\n<br> Contour plot corresponding to x'cov<sup>-1</sup>x = 4 (where x is the column vector (pij,pkl)) are drawn. \
@@ -4235,23 +4261,23 @@ To be simple, these graphs help to under
     /*j1++;*/
       if  (cptcovn>0) {
 	fprintf(ficresprob, "\n#********** Variable "); 
-	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
 	fprintf(ficresprob, "**********\n#\n");
 	fprintf(ficresprobcov, "\n#********** Variable "); 
-	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
 	fprintf(ficresprobcov, "**********\n#\n");
 	
 	fprintf(ficgp, "\n#********** Variable "); 
-	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
 	fprintf(ficgp, "**********\n#\n");
 	
 	
 	fprintf(fichtmcov, "\n<hr  size=\"2\" color=\"#EC5E5E\">********** Variable "); 
-	for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+	for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
 	fprintf(fichtmcov, "**********\n<hr size=\"2\" color=\"#EC5E5E\">");
 	
 	fprintf(ficresprobcor, "\n#********** Variable ");    
-	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
 	fprintf(ficresprobcor, "**********\n#");    
       }
       
@@ -4264,7 +4290,7 @@ To be simple, these graphs help to under
 	if(nagesqr==1)
 	  cov[3]= age*age;
 	for (k=1; k<=cptcovn;k++) {
-	  cov[2+nagesqr+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];/* j1 1 2 3 4
+	  cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];/* j1 1 2 3 4
 							 * 1  1 1 1 1
 							 * 2  2 1 1 1
 							 * 3  1 2 1 1
@@ -4272,9 +4298,9 @@ To be simple, these graphs help to under
 	  /* nbcode[1][1]=0 nbcode[1][2]=1;*/
 	}
 	/* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
-	for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2];
+	for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2];
 	for (k=1; k<=cptcovprod;k++)
-	  cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
+	  cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])];
 	
     
 	for(theta=1; theta <=npar; theta++){
@@ -4511,8 +4537,8 @@ fprintf(fichtm," \n<ul><li><b>Graphs</b>
      if (cptcovn > 0) {
        fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
        for (cpt=1; cpt<=cptcoveff;cpt++){ 
-	 fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
-	 printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);fflush(stdout);
+	 fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);
+	 printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);fflush(stdout);
        }
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
      }
@@ -4536,12 +4562,19 @@ fprintf(fichtm," \n<ul><li><b>Graphs</b>
  }/* End k1 */
  fprintf(fichtm,"</ul>");
 
-
  fprintf(fichtm,"\
 \n<br><li><h4> <a name='secondorder'>Result files (second order: variances)</a></h4>\n\
- - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n", rfileres,rfileres);
+ - Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br> \
+ - 95%% confidence intervals and Wald tests of the estimated parameters are in the log file.<br> \
+But because parameters are usually highly correlated (a higher incidence of disability \
+and a higher incidence of recovery can give very close observed transition) it might \
+be very useful to look not only at linear confidence intervals estimated from the \
+variances but at the covariance matrix. And instead of looking at the estimated coefficients \
+(parameters) of the logistic regression, it might be more meaningful to visualize the \
+covariance matrix of the one-step probabilities. \
+See page 'Matrix of variance-covariance of one-step probabilities' below. \n", rfileres,rfileres);
 
- fprintf(fichtm," - Variance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
+ fprintf(fichtm," - Standard deviation of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
 	 subdirf2(fileres,"prob"),subdirf2(fileres,"prob"));
  fprintf(fichtm,"\
  - Variance-covariance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
@@ -4587,7 +4620,7 @@ fprintf(fichtm," \n<ul><li><b>Graphs</b>
      if (cptcovn > 0) {
        fprintf(fichtm,"<hr  size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
        for (cpt=1; cpt<=cptcoveff;cpt++) 
-	 fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
+	 fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
      }
      for(cpt=1; cpt<=nlstate;cpt++) {
@@ -4795,12 +4828,15 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
 		 fprintf(ficgp," exp(p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);
 	     ij=1;/* To be checked else nbcode[0][0] wrong */
 	     for(j=3; j <=ncovmodel-nagesqr; j++) {
-	       if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { /* Bug valgrind */
-	       	 fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
-	       	 ij++;
+	       /* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */
+	       if(ij <=cptcovage) { /* Bug valgrind */
+		 if((j-2)==Tage[ij]) { /* Bug valgrind */
+		   fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]);
+		   ij++;
+		 }
 	       }
 	       else
-		 fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
+		 fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
 	     }
 	     fprintf(ficgp,")/(1");
 	     
@@ -4812,12 +4848,14 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
   
 	       ij=1;
 	       for(j=3; j <=ncovmodel-nagesqr; j++){
-		 if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
-		   fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
-		   ij++;
+		 if(ij <=cptcovage) { /* Bug valgrind */
+		   if((j-2)==Tage[ij]) { /* Bug valgrind */
+		     fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]);
+		     ij++;
+		   }
 		 }
 		 else
-		   fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
+		   fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
 	       }
 	       fprintf(ficgp,")");
 	     }
@@ -4944,7 +4982,7 @@ void prevforecast(char fileres[], double
       k=k+1;
       fprintf(ficresf,"\n#******");
       for(j=1;j<=cptcoveff;j++) {
-	fprintf(ficresf," V%d=%d, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	fprintf(ficresf," V%d=%d, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
       }
       fprintf(ficresf,"******\n");
       fprintf(ficresf,"# Covariate valuofcovar yearproj age");
@@ -4968,7 +5006,7 @@ void prevforecast(char fileres[], double
 	    if (h*hstepm/YEARM*stepm ==yearp) {
               fprintf(ficresf,"\n");
               for(j=1;j<=cptcoveff;j++) 
-                fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+                fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	      fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm);
 	    } 
 	    for(j=1; j<=nlstate+ndeath;j++) {
@@ -5066,7 +5104,7 @@ void populforecast(char fileres[], doubl
       k=k+1;
       fprintf(ficrespop,"\n#******");
       for(j=1;j<=cptcoveff;j++) {
-	fprintf(ficrespop," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	fprintf(ficrespop," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
       }
       fprintf(ficrespop,"******\n");
       fprintf(ficrespop,"# Age");
@@ -5433,8 +5471,8 @@ int readdata(char datafile[], int firsto
 
 
   if((fic=fopen(datafile,"r"))==NULL)    {
-    printf("Problem while opening datafile: %s\n", datafile);return 1;
-    fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);return 1;
+    printf("Problem while opening datafile: %s\n", datafile);fflush(stdout);
+    fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);fflush(ficlog);return 1;
   }
 
   i=1;
@@ -5723,9 +5761,9 @@ int decodemodel ( char model[], int last
       /* 	k=1 Tvar[1]=2 (from V2) */
       /* 	k=5 Tvar[5] */
       /* for (k=1; k<=cptcovn;k++) { */
-      /* 	cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; */
+      /* 	cov[2+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
       /* 	} */
-      /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2]; */
+      /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k])]]*cov[2]; */
       /*
        * Treating invertedly V2+V1+V3*age+V2*V4 is as if written V2*V4 +V3*age + V1 + V2 */
       for(k=cptcovt; k>=1;k--) /**< Number of covariates */
@@ -6136,14 +6174,14 @@ int prevalence_limit(double *p, double *
       //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
 	k=k+1;
 	/* to clean */
-	//printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtab[cptcod][cptcov]);
+	//printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
 	fprintf(ficrespl,"\n#******");
 	printf("\n#******");
 	fprintf(ficlog,"\n#******");
 	for(j=1;j<=cptcoveff;j++) {
-	  fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
-	  printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
-	  fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	  fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+	  printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
+	  fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	}
 	fprintf(ficrespl,"******\n");
 	printf("******\n");
@@ -6151,7 +6189,7 @@ int prevalence_limit(double *p, double *
 
 	fprintf(ficrespl,"#Age ");
 	for(j=1;j<=cptcoveff;j++) {
-	  fprintf(ficrespl,"V%d %d",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	  fprintf(ficrespl,"V%d %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	}
 	for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
 	fprintf(ficrespl,"\n");
@@ -6161,7 +6199,7 @@ int prevalence_limit(double *p, double *
 	  prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
 	  fprintf(ficrespl,"%.0f ",age );
 	  for(j=1;j<=cptcoveff;j++)
-	    fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	    fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 	  for(i=1; i<=nlstate;i++)
 	    fprintf(ficrespl," %.5f", prlim[i][i]);
 	  fprintf(ficrespl,"\n");
@@ -6208,7 +6246,7 @@ int hPijx(double *p, int bage, int fage)
     for (k=1; k <= (int) pow(2,cptcoveff); k++){
       fprintf(ficrespij,"\n#****** ");
       for(j=1;j<=cptcoveff;j++) 
-	fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+	fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
       fprintf(ficrespij,"******\n");
       
       for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
@@ -6260,6 +6298,7 @@ int main(int argc, char *argv[])
 
   int jj, ll, li, lj, lk;
   int numlinepar=0; /* Current linenumber of parameter file */
+  int num_filled;
   int itimes;
   int NDIM=2;
   int vpopbased=0;
@@ -6269,7 +6308,8 @@ int main(int argc, char *argv[])
   /* FILE *ficgp;*/ /*Gnuplot File */
   struct stat info;
   double agedeb=0.;
-  double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;
+
+  double ageminpar=AGEOVERFLOW,agemin=AGEOVERFLOW, agemaxpar=-AGEOVERFLOW, agemax=-AGEOVERFLOW;
 
   double fret;
   double dum=0.; /* Dummy variable */
@@ -6277,13 +6317,17 @@ int main(int argc, char *argv[])
   double ***mobaverage;
 
   char line[MAXLINE];
-  char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE];
+  char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE];
+
+  char model[MAXLINE], modeltemp[MAXLINE];
   char pathr[MAXLINE], pathimach[MAXLINE]; 
   char *tok, *val; /* pathtot */
   int firstobs=1, lastobs=10;
-  int c,  h , cpt;
+  int c,  h , cpt, c2;
   int jl=0;
   int i1, j1, jk, stepsize=0;
+  int count=0;
+
   int *tab; 
   int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */
   int mobilav=0,popforecast=0;
@@ -6355,7 +6399,7 @@ int main(int argc, char *argv[])
   getcwd(pathcd, size);
 #endif
   syscompilerinfo(0);
-  printf("\n%s\n%s",version,fullversion);
+  printf("\nIMaCh version %s, %s\n%s",version, copyright, fullversion);
   if(argc <=1){
     printf("\nEnter the parameter file name: ");
     fgets(pathr,FILENAMELENGTH,stdin);
@@ -6419,7 +6463,7 @@ int main(int argc, char *argv[])
     goto end;
   }
   fprintf(ficlog,"Log filename:%s\n",filelog);
-  fprintf(ficlog,"\n%s\n%s",version,fullversion);
+  fprintf(ficlog,"Version %s %s",version,fullversion);
   fprintf(ficlog,"\nEnter the parameter file name: \n");
   fprintf(ficlog,"pathimach=%s\npathtot=%s\n\
  path=%s \n\
@@ -6427,7 +6471,7 @@ int main(int argc, char *argv[])
  optionfilext=%s\n\
  optionfilefiname='%s'\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname);
 
-  syscompilerinfo(0);
+  syscompilerinfo(1);
 
   printf("Local time (at start):%s",strstart);
   fprintf(ficlog,"Local time (at start): %s",strstart);
@@ -6463,23 +6507,82 @@ int main(int argc, char *argv[])
 
   /* Reads comments: lines beginning with '#' */
   numlinepar=0;
-  while((c=getc(ficpar))=='#' && c!= EOF){
-    ungetc(c,ficpar);
-    fgets(line, MAXLINE, ficpar);
+
+    /* First parameter line */
+  while(fgets(line, MAXLINE, ficpar)) {
+    /* If line starts with a # it is a comment */
+    if (line[0] == '#') {
+      numlinepar++;
+      fputs(line,stdout);
+      fputs(line,ficparo);
+      fputs(line,ficlog);
+      continue;
+    }else
+      break;
+  }
+  if((num_filled=sscanf(line,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", \
+			title, datafile, &lastobs, &firstpass,&lastpass)) !=EOF){
+    if (num_filled != 5) {
+      printf("Should be 5 parameters\n");
+    }
     numlinepar++;
-    fputs(line,stdout);
-    fputs(line,ficparo);
-    fputs(line,ficlog);
+    printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", title, datafile, lastobs, firstpass,lastpass);
+  }
+  /* Second parameter line */
+  while(fgets(line, MAXLINE, ficpar)) {
+    /* If line starts with a # it is a comment */
+    if (line[0] == '#') {
+      numlinepar++;
+      fputs(line,stdout);
+      fputs(line,ficparo);
+      fputs(line,ficlog);
+      continue;
+    }else
+      break;
+  }
+  if((num_filled=sscanf(line,"ftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n", \
+			&ftol, &stepm, &ncovcol, &nlstate, &ndeath, &maxwav, &mle, &weightopt)) !=EOF){
+    if (num_filled != 8) {
+      printf("Not 8\n");
+    }
+    printf("ftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt);
   }
-  ungetc(c,ficpar);
 
-  fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);
-  numlinepar++;
-  printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model);
+  /* Third parameter line */
+  while(fgets(line, MAXLINE, ficpar)) {
+    /* If line starts with a # it is a comment */
+    if (line[0] == '#') {
+      numlinepar++;
+      fputs(line,stdout);
+      fputs(line,ficparo);
+      fputs(line,ficlog);
+      continue;
+    }else
+      break;
+  }
+  if((num_filled=sscanf(line,"model=1+age%[^.\n]\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);
+      model[0]='\0';
+      goto end;
+    }
+    else{
+      if (model[0]=='+'){
+	for(i=1; i<=strlen(model);i++)
+	  modeltemp[i-1]=model[i];
+      }
+      strcpy(model,modeltemp); 
+    }
+    printf(" model=1+age%s modeltemp= %s, model=%s\n",model, modeltemp, model);fflush(stdout);
+  }
+  /* fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); */
+  /* numlinepar=numlinepar+3; /\* In general *\/ */
+  /* printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); */
   if(model[strlen(model)-1]=='.') /* Suppressing leading dot in the model */
     model[strlen(model)-1]='\0';
-  fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
-  fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
+  fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
+  fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
   fflush(ficlog);
   /* if(model[0]=='#'|| model[0]== '\0'){ */
   if(model[0]=='#'){
@@ -6495,6 +6598,10 @@ int main(int argc, char *argv[])
     ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);
     numlinepar++;
+    if(line[1]=='q'){ /* This #q will quit imach (the answer is q) */
+      z[0]=line[1];
+    }
+    /* printf("****line [1] = %c \n",line[1]); */
     fputs(line, stdout);
     //puts(line);
     fputs(line,ficparo);
@@ -6561,7 +6668,7 @@ int main(int argc, char *argv[])
 	if(jj==i) continue;
 	j++;
 	fscanf(ficpar,"%1d%1d",&i1,&j1);
-	if ((i1 != i) && (j1 != j)){
+	if ((i1 != i) || (j1 != jj)){
 	  printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n \
 It might be a problem of design; if ncovcol and the model are correct\n \
 run imach with mle=-1 to get a correct template of the parameter file.\n",numlinepar, i,j, i1, j1);
@@ -6569,8 +6676,8 @@ run imach with mle=-1 to get a correct t
 	}
 	fprintf(ficparo,"%1d%1d",i1,j1);
 	if(mle==1)
-	  printf("%1d%1d",i,j);
-	fprintf(ficlog,"%1d%1d",i,j);
+	  printf("%1d%1d",i,jj);
+	fprintf(ficlog,"%1d%1d",i,jj);
 	for(k=1; k<=ncovmodel;k++){
 	  fscanf(ficpar," %lf",&param[i][j][k]);
 	  if(mle==1){
@@ -6651,12 +6758,22 @@ run imach with mle=-1 to get a correct t
     for(i=1; i <=npar; i++)
       for(j=1; j <=npar; j++) matcov[i][j]=0.;
       
+    /* Scans npar lines */
     for(i=1; i <=npar; i++){
-      fscanf(ficpar,"%s",str);
+      count=fscanf(ficpar,"%1d%1d%1d",&i1,&j1,&jk);
+      if(count != 3){
+	printf("Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\
+This is probably because your covariance matrix doesn't \n  contain exactly %d lines corresponding to your model line '1+age+%s'.\n\
+Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model);
+	fprintf(ficlog,"Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\
+This is probably because your covariance matrix doesn't \n  contain exactly %d lines corresponding to your model line '1+age+%s'.\n\
+Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model);
+	exit(1);
+      }else
       if(mle==1)
-	printf("%s",str);
-      fprintf(ficlog,"%s",str);
-      fprintf(ficparo,"%s",str);
+	printf("%1d%1d%1d",i1,j1,jk);
+      fprintf(ficlog,"%1d%1d%1d",i1,j1,jk);
+      fprintf(ficparo,"%1d%1d%1d",i1,j1,jk);
       for(j=1; j <=i; j++){
 	fscanf(ficpar," %le",&matcov[i][j]);
 	if(mle==1){
@@ -6672,6 +6789,7 @@ run imach with mle=-1 to get a correct t
       fprintf(ficlog,"\n");
       fprintf(ficparo,"\n");
     }
+    /* End of read covariance matrix npar lines */
     for(i=1; i <=npar; i++)
       for(j=i+1;j<=npar;j++)
 	matcov[i][j]=matcov[j][i];
@@ -6800,7 +6918,7 @@ run imach with mle=-1 to get a correct t
   /* 1 to ncodemax[j] is the maximum value of this jth covariate */
 
   codtab=imatrix(1,100,1,10); /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */
-  /*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtab[100][10]);*/
+  /*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtabm(100,10));*/
   /* codtab gives the value 1 or 2 of the hth combination of k covariates (1 or 2).*/
   h=0;
 
@@ -6810,13 +6928,6 @@ run imach with mle=-1 to get a correct t
  
   m=pow(2,cptcoveff);
  
-  for(k=1;k<=cptcoveff; k++){ /* scans any effective covariate */
-    for(i=1; i <=pow(2,cptcoveff-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 ncodemax=2*/
-	for(cpt=1; cpt <=pow(2,k-1); 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)  k   = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1
 	   * For k=4 covariates, h goes from 1 to 2**k
 	   * codtabm(h,k)=  1 & (h-1) >> (k-1) ;
@@ -6830,32 +6941,47 @@ run imach with mle=-1 to get a correct t
 	   *     6     2     1     2     1
 	   *     7 i=4 1     2     2     1
 	   *     8     2     2     2     1
-	   *     9 i=5 1 i=3 1 i=2 1     1
-	   *    10     2     1     1     1
-	   *    11 i=6 1     2     1     1
-	   *    12     2     2     1     1
-	   *    13 i=7 1 i=4 1     2     1    
-	   *    14     2     1     2     1
-	   *    15 i=8 1     2     2     1
-	   *    16     2     2     2     1
+	   *     9 i=5 1 i=3 1 i=2 1     2
+	   *    10     2     1     1     2
+	   *    11 i=6 1     2     1     2
+	   *    12     2     2     1     2
+	   *    13 i=7 1 i=4 1     2     2    
+	   *    14     2     1     2     2
+	   *    15 i=8 1     2     2     2
+	   *    16     2     2     2     2
 	   */
-	  codtab[h][k]=j;
-	  /* codtab[12][3]=1; */
-	  /*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]]);
-	} 
-      }
-    }
-  } 
+  for(h=1; h <=100 ;h++){ 
+    /* printf("h=%2d ", h); */
+     for(k=1; k <=10; k++){
+       /* printf("k=%d %d ",k,codtabm(h,k)); */
+       codtab[h][k]=codtabm(h,k);
+     }
+     /* printf("\n"); */
+  }
+  /* for(k=1;k<=cptcoveff; k++){ /\* scans any effective covariate *\/ */
+  /*   for(i=1; i <=pow(2,cptcoveff-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 ncodemax=2*\/ */
+  /* 	for(cpt=1; cpt <=pow(2,k-1); 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[12][3]=1; *\/ */
+  /* 	  /\*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]]); *\/ */
+  /* 	}  */
+  /*     } */
+  /*   } */
+  /* }  */
   /* printf("codtab[1][2]=%d codtab[2][2]=%d",codtab[1][2],codtab[2][2]); 
      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("\n");
-     }
-     scanf("%d",i);*/
+  /* 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("\n"); */
+  /* } */
+  /*   scanf("%d",i);*/
 
  free_ivector(Ndum,-1,NCOVMAX);
 
@@ -7111,9 +7237,10 @@ Interval (in months) between two waves:
     }
     
     printf("iter=%d MLE=%f Eq=%lf*exp(%lf*(age-%d))\n",iter,-gompertz(p),p[1],p[2],agegomp);
-    for (i=1;i<=NDIM;i++) 
+    for (i=1;i<=NDIM;i++) {
       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);
@@ -7145,8 +7272,15 @@ Interval (in months) between two waves:
     
     
     replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */
-    printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
-    
+    if(ageminpar == AGEOVERFLOW ||agemaxpar == AGEOVERFLOW){
+      	printf("Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\
+This is probably because your parameter file doesn't \n  contain the exact number of lines (or columns) corresponding to your model line.\n\
+Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
+      	fprintf(ficlog,"Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\
+This is probably because your parameter file doesn't \n  contain the exact number of lines (or columns) corresponding to your model line.\n\
+Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
+    }else
+      printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
     printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \
 		     stepm, weightopt,\
 		     model,imx,p,matcov,agemortsup);
@@ -7210,6 +7344,24 @@ Interval (in months) between two waves:
       ftolhess=ftol; /* Usually correct */
       hesscov(matcov, p, npar, delti, ftolhess, func);
     }
+    printf("Parameters and 95%% confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W .\n But be careful that parameters are highly correlated because incidence of disability is highly correlated to incidence of recovery.\n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n");
+    fprintf(ficlog, "Parameters, Wald tests and Wald-based confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W \n  It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n");
+    for(i=1,jk=1; i <=nlstate; i++){
+      for(k=1; k <=(nlstate+ndeath); k++){
+	if (k != i) {
+	  printf("%d%d ",i,k);
+	  fprintf(ficlog,"%d%d ",i,k);
+	  for(j=1; j <=ncovmodel; j++){
+	    printf("%12.7f W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));
+	    fprintf(ficlog,"%12.7f W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk]));
+	    jk++; 
+	  }
+	  printf("\n");
+	  fprintf(ficlog,"\n");
+	}
+      }
+    }
+
     fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");
     printf("# Scales (for hessian or gradient estimation)\n");
     fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n");
@@ -7366,6 +7518,7 @@ Interval (in months) between two waves:
     dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/365.;
     
     fscanf(ficpar,"pop_based=%d\n",&popbased);
+    fprintf(ficlog,"pop_based=%d\n",popbased);
     fprintf(ficparo,"pop_based=%d\n",popbased);   
     fprintf(ficres,"pop_based=%d\n",popbased);   
     
@@ -7390,7 +7543,15 @@ Interval (in months) between two waves:
     /* ,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); */
     
     replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */
-    printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
+    if(ageminpar == AGEOVERFLOW ||agemaxpar == -AGEOVERFLOW){
+      	printf("Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\
+This is probably because your parameter file doesn't \n  contain the exact number of lines (or columns) corresponding to your model line.\n\
+Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
+      	fprintf(ficlog,"Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\
+This is probably because your parameter file doesn't \n  contain the exact number of lines (or columns) corresponding to your model line.\n\
+Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
+    }else
+      printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
     
     printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\
 		 model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\
@@ -7793,7 +7954,7 @@ Interval (in months) between two waves:
   }
   end:
   while (z[0] != 'q') {
-    printf("\nType  q for exiting: ");
+    printf("\nType  q for exiting: "); fflush(stdout);
     scanf("%s",z);
   }
 }