--- imach/src/imach.c	2010/04/29 18:11:38	1.137
+++ imach/src/imach.c	2010/04/30 18:19:40	1.138
@@ -1,6 +1,9 @@
-/* $Id: imach.c,v 1.137 2010/04/29 18:11:38 brouard Exp $
+/* $Id: imach.c,v 1.138 2010/04/30 18:19:40 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.138  2010/04/30 18:19:40  brouard
+  *** empty log message ***
+
   Revision 1.137  2010/04/29 18:11:38  brouard
   (Module): Checking covariates for more complex models
   than V1+V2. A lot of change to be done. Unstable.
@@ -424,11 +427,11 @@ extern int errno;
 #define ODIRSEPARATOR '/'
 #endif
 
-/* $Id: imach.c,v 1.137 2010/04/29 18:11:38 brouard Exp $ */
+/* $Id: imach.c,v 1.138 2010/04/30 18:19:40 brouard Exp $ */
 /* $State: Exp $ */
 
 char version[]="Imach version 0.98m, April 2010, INED-EUROREVES-Institut de longevite ";
-char fullversion[]="$Revision: 1.137 $ $Date: 2010/04/29 18:11:38 $"; 
+char fullversion[]="$Revision: 1.138 $ $Date: 2010/04/30 18:19:40 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -1260,21 +1263,21 @@ double **prevalim(double **prlim, int nl
   for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){
     newm=savm;
     /* Covariates have to be included here again */
-     cov[2]=agefin;
-  
-      for (k=1; k<=cptcovn;k++) {
-	cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
-	/*	printf("ij=%d k=%d Tvar[k]=%d nbcode=%d cov=%lf codtab[ij][Tvar[k]]=%d \n",ij,k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+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++)
-	cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]] * nbcode[Tvard[k][2]][codtab[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]);*/
-      /*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/
+    cov[2]=agefin;
+    
+    for (k=1; k<=cptcovn;k++) {
+      cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
+      /*	printf("ij=%d k=%d Tvar[k]=%d nbcode=%d cov=%lf codtab[ij][Tvar[k]]=%d \n",ij,k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+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++)
+      cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]] * nbcode[Tvard[k][2]][codtab[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]);*/
+    /*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/
     out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);
-
+    
     savm=oldm;
     oldm=newm;
     maxmax=0.;
@@ -1301,41 +1304,56 @@ double **prevalim(double **prlim, int nl
 
 double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
 {
-  double s1, s2;
+  /* According to parameters values stored in x and the covariate's values stored in cov,
+     computes the probability to be observed in state j being in state i by appying the
+     model to the ncovmodel covariates (including constant and age).
+     lnpijopii=ln(pij/pii)= aij+bij*age+cij*v1+dij*v2+... = sum_nc=1^ncovmodel xij(nc)*cov[nc]
+     and, according on how parameters are entered, the position of the coefficient xij(nc) of the
+     ncth covariate in the global vector x is given by the formula:
+     j<i nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel
+     j>=i nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel
+     Computes ln(pij/pii) (lnpijopii), deduces pij/pii by exponentiation,
+     sums on j different of i to get 1-pii/pii, deduces pii, and then all pij.
+     Outputs ps[i][j] the probability to be observed in j being in j according to
+     the values of the covariates cov[nc] and corresponding parameter values x[nc+shiftij]
+  */
+  double s1, lnpijopii;
   /*double t34;*/
   int i,j,j1, nc, ii, jj;
 
     for(i=1; i<= nlstate; i++){
       for(j=1; j<i;j++){
-	for (nc=1, s2=0.;nc <=ncovmodel; nc++){
-	  /*s2 += param[i][j][nc]*cov[nc];*/
-	  s2 += x[(i-1)*nlstate*ncovmodel+(j-1)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
-/* 	 printf("Int j<i s1=%.17e, s2=%.17e\n",s1,s2); */
+	for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
+	  /*lnpijopii += param[i][j][nc]*cov[nc];*/
+	  lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];
+/* 	 printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
 	}
-	ps[i][j]=s2;
-/* 	printf("s1=%.17e, s2=%.17e\n",s1,s2); */
+	ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
+/* 	printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
       }
       for(j=i+1; j<=nlstate+ndeath;j++){
-	for (nc=1, s2=0.;nc <=ncovmodel; nc++){
-	  s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
-/* 	  printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2); */
+	for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
+	  /*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/
+	  lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];
+/* 	  printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */
 	}
-	ps[i][j]=s2;
+	ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
       }
     }
-    /*ps[3][2]=1;*/
     
     for(i=1; i<= nlstate; i++){
       s1=0;
       for(j=1; j<i; j++){
-	s1+=exp(ps[i][j]);
+	s1+=exp(ps[i][j]); /* In fact sums pij/pii */
 	/*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]);
+	s1+=exp(ps[i][j]); /* In fact sums pij/pii */
 	/*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
       }
+      /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
       ps[i][i]=1./(s1+1.);
+      /* Computing other pijs */
       for(j=1; j<i; j++)
 	ps[i][j]= exp(ps[i][j])*ps[i][i];
       for(j=i+1; j<=nlstate+ndeath; j++)
@@ -1466,9 +1484,14 @@ double func( double *x)
 
   if(mle==1){
     for (i=1,ipmx=0, sw=0.; i<=imx; i++){
+      /* Computes the values of the ncovmodel covariates of the model
+	 depending if the covariates are fixed or variying (age dependent) and stores them in cov[]
+	 Then computes with function pmij which return a matrix p[i][j] giving the elementary probability
+	 to be observed in j being in i according to the model.
+       */
       for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
       /* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] 
-	 is 6, Tvar[3=age*V3] should not been computed because of age Tvar[4=V3*V2] 
+	 is 6, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2] 
 	 has been calculated etc */
       for(mi=1; mi<= wav[i]-1; mi++){
 	for (ii=1;ii<=nlstate+ndeath;ii++)