--- imach/src/imach.c	2002/04/03 12:19:36	1.38
+++ imach/src/imach.c	2002/04/05 15:45:00	1.39
@@ -1,4 +1,4 @@
-/* $Id: imach.c,v 1.38 2002/04/03 12:19:36 lievre Exp $
+/* $Id: imach.c,v 1.39 2002/04/05 15:45:00 lievre Exp $
    Interpolated Markov Chain
 
   Short summary of the programme:
@@ -14,7 +14,7 @@
   model. More health states you consider, more time is necessary to reach the
   Maximum Likelihood of the parameters involved in the model.  The
   simplest model is the multinomial logistic model where pij is the
-  probabibility to be observed in state j at the second wave
+  probability to be observed in state j at the second wave
   conditional to be observed in state i at the first wave. Therefore
   the model is: log(pij/pii)= aij + bij*age+ cij*sex + etc , where
   'age' is age and 'sex' is a covariate. If you want to have a more
@@ -1503,7 +1503,7 @@ void tricode(int *Tvar, int **nbcode, in
       for (k=0; k<=19; k++) {
 	if (Ndum[k] != 0) {
 	  nbcode[Tvar[j]][ij]=k; 
-	  /*	 printf("nbcodeaaaaaaaaaaa=%d Tvar[j]=%d ij=%d j=%d",nbcode[Tvar[j]][ij],Tvar[j],ij,j);*/
+	  
 	  ij++;
 	}
 	if (ij > ncodemax[j]) break; 
@@ -1585,6 +1585,9 @@ void evsij(char fileres[], double ***eij
     /* Computed by stepm unit matrices, product of hstepm matrices, stored
        in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */
     hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, ij);  
+  
+    /*for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++) printf("%f %.5f\n", age*12+h, p3mat[1][1][h]);*/
+
     hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */
     for(i=1; i<=nlstate;i++)
       for(j=1; j<=nlstate;j++)
@@ -1822,9 +1825,9 @@ void varprevlim(char fileres[], double *
 }
 
 /************ Variance of one-step probabilities  ******************/
-void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij)
+void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax)
 {
-  int i, j;
+  int i, j, i1, k1, j1, z1;
   int k=0, cptcode;
   double **dnewm,**doldm;
   double *xp;
@@ -1847,83 +1850,101 @@ void varprob(char fileres[], double **ma
   doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath));
   
   cov[1]=1;
-  for (age=bage; age<=fage; age ++){ 
-    cov[2]=age;
-    gradg=matrix(1,npar,1,9);
-    trgradg=matrix(1,9,1,npar);
-    gp=vector(1,(nlstate+ndeath)*(nlstate+ndeath));
-    gm=vector(1,(nlstate+ndeath)*(nlstate+ndeath));
+  j=cptcoveff;
+  if (cptcovn<1) {j=1;ncodemax[1]=1;}
+  j1=0;
+  for(k1=1; k1<=1;k1++){
+    for(i1=1; i1<=ncodemax[k1];i1++){ 
+    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]]);
+      fprintf(ficresprob, "**********\n#");
+    }
     
-    for(theta=1; theta <=npar; theta++){
-      for(i=1; i<=npar; i++)
-	xp[i] = x[i] + (i==theta ?delti[theta]:0);
-     
-      pmij(pmmij,cov,ncovmodel,xp,nlstate);
-   
-      k=0;
-      for(i=1; i<= (nlstate+ndeath); i++){
-	for(j=1; j<=(nlstate+ndeath);j++){
-	   k=k+1;
-	  gp[k]=pmmij[i][j];
+      for (age=bage; age<=fage; age ++){ 
+	cov[2]=age;
+	for (k=1; k<=cptcovn;k++) {
+	  cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];
+	  
 	}
-      }
-
-      for(i=1; i<=npar; i++)
-	xp[i] = x[i] - (i==theta ?delti[theta]:0);
+	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]]];
+	
+	gradg=matrix(1,npar,1,9);
+	trgradg=matrix(1,9,1,npar);
+	gp=vector(1,(nlstate+ndeath)*(nlstate+ndeath));
+	gm=vector(1,(nlstate+ndeath)*(nlstate+ndeath));
     
-
-      pmij(pmmij,cov,ncovmodel,xp,nlstate);
-      k=0;
-      for(i=1; i<=(nlstate+ndeath); i++){
-	for(j=1; j<=(nlstate+ndeath);j++){
-	  k=k+1;
-	  gm[k]=pmmij[i][j];
-	}
-      }
+	for(theta=1; theta <=npar; theta++){
+	  for(i=1; i<=npar; i++)
+	    xp[i] = x[i] + (i==theta ?delti[theta]:0);
+	  
+	  pmij(pmmij,cov,ncovmodel,xp,nlstate);
+	  
+	  k=0;
+	  for(i=1; i<= (nlstate+ndeath); i++){
+	    for(j=1; j<=(nlstate+ndeath);j++){
+	      k=k+1;
+	      gp[k]=pmmij[i][j];
+	    }
+	  }
+	  
+	  for(i=1; i<=npar; i++)
+	    xp[i] = x[i] - (i==theta ?delti[theta]:0);
+    
+	  pmij(pmmij,cov,ncovmodel,xp,nlstate);
+	  k=0;
+	  for(i=1; i<=(nlstate+ndeath); i++){
+	    for(j=1; j<=(nlstate+ndeath);j++){
+	      k=k+1;
+	      gm[k]=pmmij[i][j];
+	    }
+	  }
      
-       for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++) 
-	   gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta];  
-    }
-
-     for(j=1; j<=(nlstate+ndeath)*(nlstate+ndeath);j++)
-      for(theta=1; theta <=npar; theta++)
-      trgradg[j][theta]=gradg[theta][j];
- 
-     matprod2(dnewm,trgradg,1,9,1,npar,1,npar,matcov);
-     matprod2(doldm,dnewm,1,9,1,npar,1,9,gradg);
-
-     pmij(pmmij,cov,ncovmodel,x,nlstate);
+	  for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++) 
+	    gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta];  
+	}
 
-     k=0;
-     for(i=1; i<=(nlstate+ndeath); i++){
-       for(j=1; j<=(nlstate+ndeath);j++){
-	 k=k+1;
-	 gm[k]=pmmij[i][j];
+	for(j=1; j<=(nlstate+ndeath)*(nlstate+ndeath);j++)
+	  for(theta=1; theta <=npar; theta++)
+	    trgradg[j][theta]=gradg[theta][j];
+	
+	matprod2(dnewm,trgradg,1,9,1,npar,1,npar,matcov);
+	matprod2(doldm,dnewm,1,9,1,npar,1,9,gradg);
+	
+	pmij(pmmij,cov,ncovmodel,x,nlstate);
+	
+	k=0;
+	for(i=1; i<=(nlstate+ndeath); i++){
+	  for(j=1; j<=(nlstate+ndeath);j++){
+	    k=k+1;
+	    gm[k]=pmmij[i][j];
+	  }
 	}
-     }
      
      /*printf("\n%d ",(int)age);
      for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){
-       
-
        printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));
      }*/
 
-  fprintf(ficresprob,"\n%d ",(int)age);
-
-  for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){
-    if (i== 2) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]);
-if (i== 4) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]);
-  }
+	fprintf(ficresprob,"\n%d ",(int)age);
 
+	for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++)
+	  fprintf(ficresprob,"%.3e (%.3e) ",gm[i],doldm[i][i]);
+  
+      }
+    }
     free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
     free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));
     free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
     free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
-}
- free_vector(xp,1,npar);
-fclose(ficresprob);
-
+  }
+  free_vector(xp,1,npar);
+  fclose(ficresprob);
+  
 }
 
 /******************* Printing html file ***********/
@@ -2751,11 +2772,10 @@ while((c=getc(ficpar))=='#' && c!= EOF){
     if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3;
     if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;
     }*/
- 
-  /* for (i=1; i<=imx; i++){
+   /*  for (i=1; i<=imx; i++){
      if (s[4][i]==9)  s[4][i]=-1; 
-     printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i]));}
-  */
+     printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i]));}*/
+  
  
   /* Calculation of the number of parameter from char model*/
   Tvar=ivector(1,15); 
@@ -3220,7 +3240,7 @@ while((c=getc(ficpar))=='#' && c!= EOF){
     }
   }
 
-  /* varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k);*/
+  varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k,Tvar,nbcode, ncodemax);
 
   fclose(ficrespij);