--- imach/src/imach.c	2001/05/17 16:07:14	1.11
+++ imach/src/imach.c	2002/02/20 17:02:08	1.13
@@ -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 
@@ -83,8 +83,8 @@ int **dh; /* dh[mi][i] is number of step
 double jmean; /* Mean space between 2 waves */
 double **oldm, **newm, **savm; /* Working pointers to matrices */
 double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
-FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest;
-FILE *ficgp, *fichtm;
+FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest,*ficresf;
+FILE *ficgp, *fichtm,*ficresprob;
 FILE *ficreseij;
   char filerese[FILENAMELENGTH];
  FILE  *ficresvij;
@@ -127,7 +127,7 @@ int stepm;
 int m,nb;
 int *num, firstpass=0, lastpass=4,*cod, *ncodemax, *Tage;
 double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;
-double **pmmij;
+double **pmmij, ***probs, ***mobaverage;
 
 double *weight;
 int **s; /* Status */
@@ -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]);
@@ -757,7 +760,7 @@ double **pmij(double **ps, double *cov,
 
 double **matprod2(double **out, double **in,long nrl, long nrh, long ncl, long nch, long ncolol, long ncoloh, double **b)
 {
-  /* Computes the matric product of in(1,nrh-nrl+1)(1,nch-ncl+1) times
+  /* Computes the matrix product of in(1,nrh-nrl+1)(1,nch-ncl+1) times
      b(1,nch-ncl+1)(1,ncoloh-ncolol+1) into out(...) */
   /* in, b, out are matrice of pointers which should have been initialized 
      before: only the contents of out is modified. The function returns
@@ -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; 
   
 }
 
@@ -1156,7 +1160,7 @@ void  freqsummary(char fileres[], int ag
   char fileresp[FILENAMELENGTH];
 
   pp=vector(1,nlstate);
-
+ probs= ma3x(1,130 ,1,8, 1,8);
   strcpy(fileresp,"p");
   strcat(fileresp,fileres);
   if((ficresp=fopen(fileresp,"w"))==NULL) {
@@ -1233,8 +1237,11 @@ void  freqsummary(char fileres[], int ag
       else
 	printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
       if( i <= (int) agemax){
-	if(pos>=1.e-5)
+	if(pos>=1.e-5){
 	  fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos);
+	  probs[i][jk][j1]= pp[jk]/pos;
+	  /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/
+	}
       else
 	  fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos);
       }
@@ -1311,7 +1318,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 +1326,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 +1343,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 +1381,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 +1456,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;
@@ -1630,7 +1638,110 @@ 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)
+{
+  int i, j;
+  int k=0, cptcode;
+  double **dnewm,**doldm;
+  double *xp;
+  double *gp, *gm;
+  double **gradg, **trgradg;
+  double age,agelim, cov[NCOVMAX];
+  int theta;
+  char fileresprob[FILENAMELENGTH];
+
+  strcpy(fileresprob,"prob"); 
+  strcat(fileresprob,fileres);
+  if((ficresprob=fopen(fileresprob,"w"))==NULL) {
+    printf("Problem with resultfile: %s\n", fileresprob);
+  }
+  printf("Computing variance of one-step probabilities: result on file '%s' \n",fileresprob);
+  
 
+  xp=vector(1,npar);
+  dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
+  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));
+    
+    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);
+
+     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]);
+  }
+
+    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);
+ exit(0);
+}
 
 /***********************************************/
 /**************** Main Program *****************/
@@ -1653,7 +1764,7 @@ int main()
   char line[MAXLINE], linepar[MAXLINE];
   char title[MAXLINE];
   char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH];
-  char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH];
+  char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], fileresf[FILENAMELENGTH];
   char filerest[FILENAMELENGTH];
   char fileregp[FILENAMELENGTH];
   char path[80],pathc[80],pathcd[80],pathtot[80],model[20];
@@ -1663,7 +1774,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;
@@ -1677,9 +1788,12 @@ int main()
   double ***eij, ***vareij;
   double **varpl; /* Variances of prevalence limits by age */
   double *epj, vepp;
+  double kk1;
+
   char version[80]="Imach version 64b, May 2001, INED-EUROREVES ";
   char *alph[]={"a","a","b","c","d","e"}, str[4];
 
+
   char z[1]="c", occ;
 #include <sys/time.h>
 #include <time.h>
@@ -1772,7 +1886,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 +1977,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 +1999,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 +2098,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 +2213,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) */
@@ -2279,7 +2409,7 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
       fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);   
       fprintf(ficgp,"set out \"p%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
     } 
-  }
+  }  
 
   /* proba elementaires */
    for(i=1,jk=1; i <=nlstate; i++){
@@ -2480,6 +2610,7 @@ fclose(fichtm);
       }
     }
   fclose(ficrespl);
+
   /*------------- h Pij x at various ages ------------*/
   
   strcpy(filerespij,"pij");  strcat(filerespij,fileres);
@@ -2489,7 +2620,7 @@ fclose(fichtm);
   printf("Computing pij: result on file '%s' \n", filerespij);
   
   stepsize=(int) (stepm+YEARM-1)/YEARM;
-  if (stepm<=24) stepsize=2;
+  /*if (stepm<=24) stepsize=2;*/
 
   agelim=AGESUP;
   hstepm=stepsize*YEARM; /* Every year of age */
@@ -2528,8 +2659,110 @@ fclose(fichtm);
     }
   }
 
+  /* varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k);*/
+
   fclose(ficrespij);
 
+  exit(0);
+  /*---------- Forecasting ------------------*/
+
+  strcpy(fileresf,"f"); 
+  strcat(fileresf,fileres);
+  if((ficresf=fopen(fileresf,"w"))==NULL) {
+    printf("Problem with forecast resultfile: %s\n", fileresf);goto end;
+  }
+  printf("Computing forecasting: result on file '%s' \n", fileresf);
+
+  /* Mobile average */
+
+  /* for (agedeb=bage; agedeb<=fage; agedeb++)
+    for (i=1; i<=nlstate;i++)
+      for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++)
+      printf("%f %d i=%d j1=%d\n", probs[(int)agedeb][i][cptcod],(int) agedeb,i,cptcod);*/
+
+  if (cptcoveff==0) ncodemax[cptcoveff]=1;
+
+  mobaverage= ma3x(1,130 ,1,8, 1,8);
+  for (agedeb=bage+3; agedeb<=fage-2; agedeb++)
+    for (i=1; i<=nlstate;i++)
+      for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)
+	mobaverage[(int)agedeb][i][cptcod]=0.;
+  
+  for (agedeb=bage+4; agedeb<=fage; agedeb++){
+    for (i=1; i<=nlstate;i++){
+      for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
+   	for (cpt=0;cpt<=4;cpt++){
+	  mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]+probs[(int)agedeb-cpt][i][cptcod];
+	  }
+	  mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]/5;
+      }
+    }
+  }
+
+/* if (cptcod==2) printf("m=%f p=%f %d age=%d ",mobaverage[(int)agedeb-2][i][cptcod],probs[(int)agedeb-cpt][i][cptcod],cpt,(int)agedeb-2);*/
+
+
+  stepsize=(int) (stepm+YEARM-1)/YEARM;
+  if (stepm<=24) stepsize=2;
+
+  agelim=AGESUP;
+  hstepm=stepsize*YEARM; /* Every year of age */
+  hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ 
+  hstepm=12;
+   k=0;
+  for(cptcov=1;cptcov<=i1;cptcov++){
+    for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
+      k=k+1;
+      fprintf(ficresf,"\n#****** ");
+      for(j=1;j<=cptcoveff;j++) {
+	fprintf(ficresf,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+      }
+      
+      fprintf(ficresf,"******\n");
+
+      fprintf(ficresf,"# StartingAge FinalAge Horizon(in years)");
+      for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);
+
+      for (agedeb=fage; agedeb>=bage; agedeb--){ 
+	fprintf(ficresf,"\n%d %.f %.f 0 ",k,agedeb, agedeb);
+	for(j=1; j<=nlstate;j++) 
+	  fprintf(ficresf,"%.3f ",mobaverage[(int)agedeb][j][cptcod]);
+      }
+      for(j=1; j<=ndeath;j++) fprintf(ficresf,"0.");
+
+      for (cpt=1; cpt<=8;cpt++)   
+      for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
+	
+	nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ 
+	nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
+	/*printf("stepm=%d hstepm=%d nhstepm=%d \n",stepm,hstepm,nhstepm);*/
+
+	p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+	oldm=oldms;savm=savms;
+	hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  
+		
+	for (h=0; h<=nhstepm; h++){
+	
+	  if (h*hstepm/YEARM*stepm==cpt)
+ fprintf(ficresf,"\n%d %.f %.f %.f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm, h*hstepm/YEARM*stepm);
+	  
+	  for(j=1; j<=nlstate+ndeath;j++) {
+	    kk1=0.;
+	    for(i=1; i<=nlstate;i++) {	      
+	      /*   kk1=kk1+p3mat[i][j][h]*probs[(int)agedeb][i][cptcod];*/
+		kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb][i][cptcod];
+	    }
+	    
+	    if (h*hstepm/YEARM*stepm==cpt) 
+	      fprintf(ficresf," %.5f ", kk1);
+	  }
+	  }
+	free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+	}
+      }
+    }
+  fclose(ficresf);
+
   /*---------- Health expectancies and variances ------------*/
 
   strcpy(filerest,"t");
@@ -2608,6 +2841,9 @@ fclose(fichtm);
     }
   }
 	
+       
+
+
  fclose(ficreseij);
  fclose(ficresvij);
   fclose(ficrest);
@@ -2653,7 +2889,7 @@ strcpy(fileresvpl,"vpl");
   free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
   free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
   free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);
-  
+ 
   free_matrix(matcov,1,npar,1,npar);
   free_vector(delti,1,npar);
   
@@ -2666,13 +2902,12 @@ strcpy(fileresvpl,"vpl");
   /*printf("Total time was %d uSec.\n", total_usecs);*/
   /*------ End -----------*/
 
+
  end:
 #ifdef windows
  chdir(pathcd);
 #endif 
- /*system("wgnuplot graph.plt");*/
- /*system("../gp37mgw/wgnuplot graph.plt");*/
- /*system("cd ../gp37mgw");*/
+ 
  system("..\\gp37mgw\\wgnuplot graph.plt");
 
 #ifdef windows