--- imach/src/imach.c	2001/05/17 16:07:14	1.11
+++ imach/src/imach.c	2002/02/20 16:57:00	1.12
@@ -22,7 +22,7 @@
   delay between waves is not identical for each individual, or if some 
   individual missed an interview, the information is not rounded or lost, but
   taken into account using an interpolation or extrapolation.
-  hPijx is the probability to be
+  hPijx is the probability to be 
   observed in state i at age x+h conditional to the observed state i at age 
   x. The delay 'h' can be split into an exact number (nh*stepm) of 
   unobserved intermediate  states. This elementary transition (by month or 
@@ -693,7 +693,7 @@ double **prevalim(double **prlim, int nl
   }
 }
 
-/*************** transition probabilities **********/ 
+/*************** transition probabilities ***************/ 
 
 double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
 {
@@ -716,9 +716,11 @@ double **pmij(double **ps, double *cov,
 	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);*/
       }
-      ps[i][j]=s2;
+      ps[i][j]=(s2);
     }
   }
+    /*ps[3][2]=1;*/
+
   for(i=1; i<= nlstate; i++){
      s1=0;
     for(j=1; j<i; j++)
@@ -740,6 +742,7 @@ double **pmij(double **ps, double *cov,
     }
   }
 
+
   /*   for(ii=1; ii<= nlstate+ndeath; ii++){
     for(jj=1; jj<= nlstate+ndeath; jj++){
      printf("%lf ",ps[ii][jj]);
@@ -804,9 +807,9 @@ double ***hpxij(double ***po, int nhstep
       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<=cptcovage;k++)
+      for (k=1; k<=cptcovage;k++)
 	cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
-   for (k=1; k<=cptcovprod;k++)
+      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]]];
 
 
@@ -913,7 +916,6 @@ void hesscov(double **matcov, double p[]
   void lubksb(double **a, int npar, int *indx, double b[]) ;
   void ludcmp(double **a, int npar, int *indx, double *d) ;
 
-
   hess=matrix(1,npar,1,npar);
 
   printf("\nCalculation of the hessian matrix. Wait...\n");
@@ -921,14 +923,16 @@ void hesscov(double **matcov, double p[]
     printf("%d",i);fflush(stdout);
     hess[i][i]=hessii(p,ftolhess,i,delti);
     /*printf(" %f ",p[i]);*/
+    /*printf(" %lf ",hess[i][i]);*/
   }
-
+  
   for (i=1;i<=npar;i++) {
     for (j=1;j<=npar;j++)  {
       if (j>i) { 
 	printf(".%d%d",i,j);fflush(stdout);
 	hess[i][j]=hessij(p,delti,i,j);
-	hess[j][i]=hess[i][j];
+	hess[j][i]=hess[i][j];    
+	/*printf(" %lf ",hess[i][j]);*/
       }
     }
   }
@@ -1032,7 +1036,7 @@ double hessii( double x[], double delta,
     }
   }
   delti[theta]=delts;
-  return res;
+  return res; 
   
 }
 
@@ -1311,7 +1315,7 @@ void  concatwav(int wav[], int **dh, int
 	  if (j >= jmax) jmax=j;
 	  if (j <= jmin) jmin=j;
 	  sum=sum+j;
-	  if (j<0) printf("j=%d num=%d ",j,i); 
+	  /* if (j<10) printf("j=%d num=%d ",j,i); */
 	  }
 	}
 	else{
@@ -1319,6 +1323,7 @@ void  concatwav(int wav[], int **dh, int
 	  k=k+1;
 	  if (j >= jmax) jmax=j;
 	  else if (j <= jmin)jmin=j;
+	  /*   if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
 	  sum=sum+j;
 	}
 	jk= j/stepm;
@@ -1335,7 +1340,7 @@ void  concatwav(int wav[], int **dh, int
   }
   jmean=sum/k;
   printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);
-}
+ }
 /*********** Tricode ****************************/
 void tricode(int *Tvar, int **nbcode, int imx)
 {
@@ -1373,7 +1378,7 @@ void tricode(int *Tvar, int **nbcode, in
 
  for (k=0; k<19; k++) Ndum[k]=0;
 
- for (i=1; i<=ncovmodel; i++) {
+ for (i=1; i<=ncovmodel-2; i++) {
       ij=Tvar[i];
       Ndum[ij]++; 
     }
@@ -1448,7 +1453,7 @@ void varevsij(char fileres[], double ***
   double **dnewm,**doldm;
   int i, j, nhstepm, hstepm, h;
   int k, cptcode;
-   double *xp;
+  double *xp;
   double **gp, **gm;
   double ***gradg, ***trgradg;
   double ***p3mat;
@@ -1663,7 +1668,7 @@ int main()
   int ju,jl, mi;
   int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;
   int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
-  
+ 
   int hstepm, nhstepm;
   double bage, fage, age, agelim, agebase;
   double ftolpl=FTOL;
@@ -1772,7 +1777,8 @@ split(pathtot, path,optionfile);
       fprintf(ficparo,"\n");
     }
   
-  npar= (nlstate+ndeath-1)*nlstate*ncovmodel;
+    npar= (nlstate+ndeath-1)*nlstate*ncovmodel;
+
   p=param[1][1];
   
   /* Reads comments: lines beginning with '#' */
@@ -1862,7 +1868,7 @@ split(pathtot, path,optionfile);
     tab=ivector(1,NCOVMAX);
     ncodemax=ivector(1,8);
 
-    i=1; 
+    i=1;
     while (fgets(line, MAXLINE, fic) != NULL)    {
       if ((i >= firstobs) && (i <=lastobs)) {
 	
@@ -1884,16 +1890,24 @@ split(pathtot, path,optionfile);
 	  cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);
 	} 
 	num[i]=atol(stra);
-
-	/*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]));*/
+	
+	/*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){
+	  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])); ij=ij+1;}*/
 
 	i=i+1;
       }
     } 
-
-    /*scanf("%d",i);*/
+    /* printf("ii=%d", ij);
+       scanf("%d",i);*/
   imx=i-1; /* Number of individuals */
 
+  /* for (i=1; i<=imx; i++){
+    if ((s[1][i]==3) && (s[2][i]==2)) s[2][i]=3;
+    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++) 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); 
   Tprod=ivector(1,15); 
@@ -1975,6 +1989,13 @@ split(pathtot, path,optionfile);
     }
     /*-calculation of age at interview from date of interview and age at death -*/
     agev=matrix(1,maxwav,1,imx);
+
+   for (i=1; i<=imx; i++) 
+     for(m=2; (m<= maxwav); m++)
+       if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){
+	 anint[m][i]=9999;
+	 s[m][i]=-1;
+       }
     
     for (i=1; i<=imx; i++)  {
       agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
@@ -2083,7 +2104,7 @@ printf("Total number of individuals= %d,
     newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
     savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
     oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
-    
+     
     /* For Powell, parameters are in a vector p[] starting at p[1]
        so we point p on param[1][1] so that p[1] maps on param[1][1][1] */
     p=param[1][1]; /* *(*(*(param +1)+1)+0) */
@@ -2666,6 +2687,7 @@ strcpy(fileresvpl,"vpl");
   /*printf("Total time was %d uSec.\n", total_usecs);*/
   /*------ End -----------*/
 
+
  end:
 #ifdef windows
  chdir(pathcd);