--- imach/src/imach.c	2002/02/20 16:57:00	1.12
+++ imach/src/imach.c	2002/02/20 17:08:52	1.15
@@ -73,6 +73,7 @@ int npar=NPARMAX;
 int nlstate=2; /* Number of live states */
 int ndeath=1; /* Number of dead states */
 int ncovmodel, ncov;     /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */
+int popbased=0;
 
 int *wav; /* Number of waves for this individuual 0 is possible */
 int maxwav; /* Maxim number of waves */
@@ -83,8 +84,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 +128,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 */
@@ -760,7 +761,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
@@ -1149,7 +1150,7 @@ void lubksb(double **a, int n, int *indx
 } 
 
 /************ Frequencies ********************/
-void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax)
+void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax, int fprev1,int lprev1)
 {  /* Some frequencies */
  
   int i, m, jk, k1, i1, j1, bool, z1,z2,j;
@@ -1160,7 +1161,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) {
@@ -1191,7 +1192,7 @@ void  freqsummary(char fileres[], int ag
 	       bool=0;
 	 }
 	  if (bool==1) {
-	   for(m=firstpass; m<=lastpass-1; m++){
+	   for(m=fprev1; m<=lprev1; m++){
 	     if(agev[m][i]==0) agev[m][i]=agemax+1;
 	     if(agev[m][i]==1) agev[m][i]=agemax+2;
 	     freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
@@ -1215,7 +1216,7 @@ void  freqsummary(char fileres[], int ag
       printf("Age %d", i);
     for(jk=1; jk <=nlstate ; jk++){
       for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
-	pp[jk] += freq[jk][m][i];
+	pp[jk] += freq[jk][m][i]; 
     }
     for(jk=1; jk <=nlstate ; jk++){
       for(m=-1, pos=0; m <=0 ; m++)
@@ -1225,10 +1226,12 @@ void  freqsummary(char fileres[], int ag
       else
         printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
     }
-    for(jk=1; jk <=nlstate ; jk++){
-      for(m=1, pp[jk]=0; m <=nlstate+ndeath; m++)
+
+     for(jk=1; jk <=nlstate ; jk++){
+      for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)
 	pp[jk] += freq[jk][m][i];
-    }
+     }
+
     for(jk=1,pos=0; jk <=nlstate ; jk++)
       pos += pp[jk];
     for(jk=1; jk <=nlstate ; jk++){
@@ -1237,8 +1240,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);
       }
@@ -1259,6 +1265,83 @@ void  freqsummary(char fileres[], int ag
 
 }  /* End of Freq */
 
+/************ Prevalence ********************/
+void prevalence(int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax, int fprev1,int lprev1)
+{  /* Some frequencies */
+ 
+  int i, m, jk, k1, i1, j1, bool, z1,z2,j;
+  double ***freq; /* Frequencies */
+  double *pp;
+  double pos;
+
+  pp=vector(1,nlstate);
+  probs= ma3x(1,130 ,1,8, 1,8);
+  
+  freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);
+  j1=0;
+  
+  j=cptcoveff;
+  if (cptcovn<1) {j=1;ncodemax[1]=1;}
+  
+ for(k1=1; k1<=j;k1++){
+    for(i1=1; i1<=ncodemax[k1];i1++){
+      j1++;
+ 
+      for (i=-1; i<=nlstate+ndeath; i++)  
+	for (jk=-1; jk<=nlstate+ndeath; jk++)  
+	  for(m=agemin; m <= agemax+3; m++)
+	  freq[i][jk][m]=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 (bool==1) {
+	  for(m=fprev1; m<=lprev1; m++){
+	    if(agev[m][i]==0) agev[m][i]=agemax+1;
+	    if(agev[m][i]==1) agev[m][i]=agemax+2;
+	    freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
+	    freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];
+	  }
+	}
+      }
+       for(i=(int)agemin; i <= (int)agemax+3; i++){ 
+	for(jk=1; jk <=nlstate ; jk++){
+	  for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
+	    pp[jk] += freq[jk][m][i]; 
+	}
+	for(jk=1; jk <=nlstate ; jk++){
+	  for(m=-1, pos=0; m <=0 ; m++)
+	    pos += freq[jk][m][i];
+	}
+	
+	 for(jk=1; jk <=nlstate ; jk++){
+	   for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)
+	     pp[jk] += freq[jk][m][i];
+	 }
+	 
+	 for(jk=1,pos=0; jk <=nlstate ; jk++) pos += pp[jk];
+
+	 for(jk=1; jk <=nlstate ; jk++){	   
+	   if( i <= (int) agemax){
+	     if(pos>=1.e-5){
+	       probs[i][jk][j1]= pp[jk]/pos;
+	     }
+	   }
+	 }
+	 
+	 }
+    }
+  }
+  
+  
+  free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
+  free_vector(pp,1,nlstate);
+  
+}  /* End of Freq */
 /************* Waves Concatenation ***************/
 
 void  concatwav(int wav[], int **dh, int **mw, int **s, double *agedc, double **agev, int  firstpass, int lastpass, int imx, int nlstate, int stepm)
@@ -1489,6 +1572,12 @@ void varevsij(char fileres[], double ***
       }
       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
+
+      if (popbased==1) {
+	for(i=1; i<=nlstate;i++)
+	  prlim[i][i]=probs[(int)age][i][ij];
+      }
+      
       for(j=1; j<= nlstate; j++){
 	for(h=0; h<=nhstepm; h++){
 	  for(i=1, gp[h][j]=0.;i<=nlstate;i++)
@@ -1500,12 +1589,19 @@ void varevsij(char fileres[], double ***
 	xp[i] = x[i] - (i==theta ?delti[theta]:0);
       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
+
+      if (popbased==1) {
+	for(i=1; i<=nlstate;i++)
+	  prlim[i][i]=probs[(int)age][i][ij];
+      }
+
       for(j=1; j<= nlstate; j++){
 	for(h=0; h<=nhstepm; h++){
 	  for(i=1, gm[h][j]=0.;i<=nlstate;i++)
 	    gm[h][j] += prlim[i][i]*p3mat[i][j][h];
 	}
       }
+
       for(j=1; j<= nlstate; j++)
 	for(h=0; h<=nhstepm; h++){
 	  gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
@@ -1635,7 +1731,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 *****************/
@@ -1658,7 +1857,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];
@@ -1667,9 +1866,10 @@ int main()
   int c,  h , cpt,l;
   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 jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab; 
+  int mobilav=0, fprev, lprev ,fprevfore=1, lprevfore=1,nforecast;
   int hstepm, nhstepm;
+
   double bage, fage, age, agelim, agebase;
   double ftolpl=FTOL;
   double **prlim;
@@ -1682,9 +1882,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>
@@ -1743,7 +1946,25 @@ split(pathtot, path,optionfile);
   fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncov, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);
   printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate,ndeath, maxwav, mle, weightopt,model);
   fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncov,nlstate,ndeath,maxwav, mle, weightopt,model);
-
+while((c=getc(ficpar))=='#' && c!= EOF){
+    ungetc(c,ficpar);
+    fgets(line, MAXLINE, ficpar);
+    puts(line);
+    fputs(line,ficparo);
+  }
+  ungetc(c,ficpar);
+  
+  fscanf(ficpar,"fprevalence=%d lprevalence=%d pop_based=%d\n",&fprev,&lprev,&popbased);
+ while((c=getc(ficpar))=='#' && c!= EOF){
+    ungetc(c,ficpar);
+    fgets(line, MAXLINE, ficpar);
+    puts(line);
+    fputs(line,ficparo);
+  }
+  ungetc(c,ficpar);
+  
+  fscanf(ficpar,"fprevalence=%d lprevalence=%d nforecast=%d mob_average=%d\n",&fprevfore,&lprevfore,&nforecast,&mobilav);
+ 
   covar=matrix(0,NCOVMAX,1,n); 
   cptcovn=0; 
   if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;
@@ -1905,8 +2126,8 @@ split(pathtot, path,optionfile);
     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]));*/
+    }
+    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); 
@@ -2097,7 +2318,7 @@ printf("Total number of individuals= %d,
     
    /* Calculates basic frequencies. Computes observed prevalence at single age
        and prints on file fileres'p'. */
-  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax);
+  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax, fprev, lprev);
 
     pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
     oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
@@ -2300,7 +2521,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++){
@@ -2361,17 +2582,12 @@ ij=1;
   fclose(ficgp);
    
 chdir(path);
-    free_matrix(agev,1,maxwav,1,imx);
+   
     free_ivector(wav,1,imx);
     free_imatrix(dh,1,lastpass-firstpass+1,1,imx);
-    free_imatrix(mw,1,lastpass-firstpass+1,1,imx);
-    
-    free_imatrix(s,1,maxwav+1,1,n);
-    
-    
+    free_imatrix(mw,1,lastpass-firstpass+1,1,imx);   
     free_ivector(num,1,n);
     free_vector(agedc,1,n);
-    free_vector(weight,1,n);
     /*free_matrix(covar,1,NCOVMAX,1,n);*/
     fclose(ficparo);
     fclose(ficres);
@@ -2416,7 +2632,9 @@ Interval (in months) between two waves:
         - Life expectancies by age and initial health status: <a href=\"e%s\">e%s</a> <br>
         - Variances of life expectancies by age and initial health status: <a href=\"v%s\">v%s</a><br>
         - Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>
-        - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br><br>",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);
+        - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>
+        - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>
+<br>",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);
 
  fprintf(fichtm," <li>Graphs</li><p>");
 
@@ -2501,6 +2719,7 @@ fclose(fichtm);
       }
     }
   fclose(ficrespl);
+
   /*------------- h Pij x at various ages ------------*/
   
   strcpy(filerespij,"pij");  strcat(filerespij,fileres);
@@ -2510,7 +2729,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 */
@@ -2549,8 +2768,111 @@ fclose(fichtm);
     }
   }
 
+  /* varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k);*/
+
   fclose(ficrespij);
 
+  /*---------- 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);
+
+  prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax, fprevfore, lprevfore);
+
+ free_matrix(agev,1,maxwav,1,imx);
+  /* Mobile average */
+
+  if (cptcoveff==0) ncodemax[cptcoveff]=1;
+
+  if (mobilav==1) {
+    mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+    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;
+	}
+      }
+    }   
+  }
+
+  stepsize=(int) (stepm+YEARM-1)/YEARM;
+  if (stepm<=12) stepsize=1;
+
+  agelim=AGESUP;
+  hstepm=stepsize*YEARM; /* Every year of age */
+  hstepm=hstepm/stepm; /* Typically 2 years, = 2 years/6 months = 4 */ 
+
+   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);
+       if (mobilav==1) {
+	for(j=1; j<=nlstate;j++) 
+	  fprintf(ficresf," %.5f ",mobaverage[(int)agedeb][j][cptcod]);
+	}
+	else {
+	  for(j=1; j<=nlstate;j++) 
+	  fprintf(ficresf," %.5f ",probs[(int)agedeb][j][cptcod]);
+	}  
+      for(j=1; j<=ndeath;j++) fprintf(ficresf," 0.00000");
+      }
+      for (cpt=1; cpt<=nforecast;cpt++)   
+      for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
+	
+	nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); 
+	nhstepm = nhstepm/hstepm; 
+	/*printf("agedeb=%.lf stepm=%d hstepm=%d nhstepm=%d \n",agedeb,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++) {	      
+	      if (mobilav==1)
+	      kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb][i][cptcod];
+	      else kk1=kk1+p3mat[i][j][h]*probs[(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);
+      }
+    }
+  }
+  if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+  free_imatrix(s,1,maxwav+1,1,n);
+  free_vector(weight,1,n);
+  fclose(ficresf);
   /*---------- Health expectancies and variances ------------*/
 
   strcpy(filerest,"t");
@@ -2610,6 +2932,11 @@ fclose(fichtm);
       epj=vector(1,nlstate+1);
       for(age=bage; age <=fage ;age++){
 	prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
+	if (popbased==1) {
+	  for(i=1; i<=nlstate;i++)
+	    prlim[i][i]=probs[(int)age][i][k];
+	}
+	
 	fprintf(ficrest," %.0f",age);
 	for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
 	  for(i=1, epj[j]=0.;i <=nlstate;i++) {
@@ -2629,6 +2956,9 @@ fclose(fichtm);
     }
   }
 	
+       
+
+
  fclose(ficreseij);
  fclose(ficresvij);
   fclose(ficrest);
@@ -2674,7 +3004,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);
   
@@ -2692,9 +3022,7 @@ strcpy(fileresvpl,"vpl");
 #ifdef windows
  chdir(pathcd);
 #endif 
- /*system("wgnuplot graph.plt");*/
- /*system("../gp37mgw/wgnuplot graph.plt");*/
- /*system("cd ../gp37mgw");*/
+ 
  system("..\\gp37mgw\\wgnuplot graph.plt");
 
 #ifdef windows