--- imach/src/imach.c	2014/01/26 09:45:38	1.143
+++ imach/src/imach.c	2014/02/10 22:17:31	1.144
@@ -1,6 +1,9 @@
-/* $Id: imach.c,v 1.143 2014/01/26 09:45:38 brouard Exp $
+/* $Id: imach.c,v 1.144 2014/02/10 22:17:31 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.144  2014/02/10 22:17:31  brouard
+  *** empty log message ***
+
   Revision 1.143  2014/01/26 09:45:38  brouard
   Summary: Version 0.98nR (to be improved, but gives same optimization results as 0.98k. Nice, promising
 
@@ -429,18 +432,18 @@ extern int errno;
 #define	GLOCK_ERROR_NOPATH		-1	/* empty path */
 #define	GLOCK_ERROR_GETCWD		-2	/* cannot get cwd */
 
-#define MAXPARM 128 /* Maximum number of parameters for the optimization */
-#define NPARMAX 64 /* (nlstate+ndeath-1)*nlstate*ncovmodel */
+#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 20 /* Maximum number of covariates */
+#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 MAXN 20000
-#define YEARM 12. /* Number of months per year */
+#define YEARM 12. /**< Number of months per year */
 #define AGESUP 130
 #define AGEBASE 40
-#define AGEGOMP 10. /* Minimal age for Gompertz adjustment */
+#define AGEGOMP 10. /**< Minimal age for Gompertz adjustment */
 #ifdef UNIX
 #define DIRSEPARATOR '/'
 #define CHARSEPARATOR "/"
@@ -451,11 +454,11 @@ extern int errno;
 #define ODIRSEPARATOR '/'
 #endif
 
-/* $Id: imach.c,v 1.143 2014/01/26 09:45:38 brouard Exp $ */
+/* $Id: imach.c,v 1.144 2014/02/10 22:17:31 brouard Exp $ */
 /* $State: Exp $ */
 
 char version[]="Imach version 0.98nR2, January 2014,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121)";
-char fullversion[]="$Revision: 1.143 $ $Date: 2014/01/26 09:45:38 $"; 
+char fullversion[]="$Revision: 1.144 $ $Date: 2014/02/10 22:17:31 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -2209,7 +2212,7 @@ void  freqsummary(char fileres[], int ia
 
   first=1;
 
-  for(k1=1; k1<=j;k1++){   /* Loop on covariates */
+  for(k1=1; k1<=j ; k1++){   /* Loop on covariates */
     for(i1=1; i1<=ncodemax[k1];i1++){ /* Now it is 2 */
       j1++;
       /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
@@ -2227,11 +2230,17 @@ void  freqsummary(char fileres[], int ia
       k2cpt=0;
       for (i=1; i<=imx; i++) {
 	bool=1;
-	if  (cptcovn>0) {
-	  for (z1=1; z1<=cptcoveff; z1++) 
-	    if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) 
-	      bool=0;
+	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]]){
+              bool=0;
+              printf("bool=%d i=%d, z1=%d, i1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtab[%d][%d]=%d, nbcode[Tvaraff][codtab[%d][%d]=%d, j1=%d\n", 
+                bool,i,z1, i1, 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*/
+            } 
 	}
+ 
 	if (bool==1){
 	  for(m=firstpass; m<=lastpass; m++){
 	    k2=anint[m][i]+(mint[m][i]/12.);
@@ -2573,15 +2582,18 @@ void  concatwav(int wav[], int **dh, int
 /*********** Tricode ****************************/
 void tricode(int *Tvar, int **nbcode, int imx)
 {
-  /* Uses cptcovn+2*cptcovprod as the number of covariates */
-  /*	  Tvar[i]=atoi(stre); /* find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 */
+  /**< Uses cptcovn+2*cptcovprod as the number of covariates */
+  /*	  Tvar[i]=atoi(stre);  find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 
+  /* Boring subroutine which should only output nbcode[Tvar[j]][k]
+  /* nbcode[Tvar[j][1]= 
+  */
 
   int Ndum[20],ij=1, k=0, j=0, i=0, maxncov=NCOVMAX;
   int modmaxcovj=0; /* Modality max of covariates j */
   cptcoveff=0; 
  
-  for (k=0; k<maxncov; k++) Ndum[k]=0;
-  for (k=1; k<=7; k++) ncodemax[k]=0; /* Horrible constant again */
+  for (k=0; k < maxncov; k++) Ndum[k]=0;
+  for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */
 
   for (j=1; j<=(cptcovn+2*cptcovprod); j++) { /* For each covariate j */
     for (i=1; i<=imx; i++) { /*reads the data file to get the maximum value of the 
@@ -2635,7 +2647,7 @@ void tricode(int *Tvar, int **nbcode, in
    }
  }
  ij--;
- cptcoveff=ij; /*Number of simple covariates*/
+ cptcoveff=ij; /*Number of total covariates*/
 }
 
 /*********** Health Expectancies ****************/
@@ -5420,7 +5432,7 @@ run imach with mle=-1 to get a correct t
   anint=matrix(1,maxwav,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); /* hard coded ? */
+  ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */
 
   /* Reads data from file datafile */
   if (readdata(datafile, firstobs, lastobs, &imx)==1)
@@ -5499,8 +5511,8 @@ run imach with mle=-1 to get a correct t
   ncodemax[1]=1;
   if (cptcovn > 0) tricode(Tvar,nbcode,imx);
       
-  codtab=imatrix(1,100,1,10); /* Cross tabulation to get the order of 
-				 the estimations*/
+  codtab=imatrix(1,100,1,10); /**< codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) 
+			       */
   h=0;
   m=pow(2,cptcoveff);
  
@@ -5511,7 +5523,7 @@ run imach with mle=-1 to get a correct t
 	  h++;
 	  if (h>m) 
 	    h=1;
-	  /**< codtab(h,k)  k
+	  /**< codtab(h,k)  k   = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1
 	   *     h     1     2     3     4
 	   *______________________________  
 	   *     1 i=1 1 i=1 1 i=1 1 i=1 1