--- imach/src/imach.c	2002/02/20 17:05:44	1.14
+++ imach/src/imach.c	2002/02/20 17:17:09	1.18
@@ -73,7 +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, fprev,lprev;
+int popbased=0;
 
 int *wav; /* Number of waves for this individuual 0 is possible */
 int maxwav; /* Maxim number of waves */
@@ -85,7 +85,7 @@ double jmean; /* Mean space between 2 wa
 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,*ficresf;
-FILE *ficgp, *fichtm,*ficresprob;
+FILE *ficgp, *fichtm,*ficresprob,*ficpop;
 FILE *ficreseij;
   char filerese[FILENAMELENGTH];
  FILE  *ficresvij;
@@ -1150,18 +1150,18 @@ 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,double **mint,double **anint, int boolprev, double dateprev1,double dateprev2)
 {  /* Some frequencies */
  
-  int i, m, jk, k1, i1, j1, bool, z1,z2,j;
+  int i, m, jk, k1,i1, j1, bool, z1,z2,j;
   double ***freq; /* Frequencies */
   double *pp;
-  double pos;
+  double pos, k2;
   FILE *ficresp;
   char fileresp[FILENAMELENGTH];
 
   pp=vector(1,nlstate);
- probs= ma3x(1,130 ,1,8, 1,8);
+  probs= ma3x(1,130 ,1,8, 1,8);
   strcpy(fileresp,"p");
   strcat(fileresp,fileres);
   if((ficresp=fopen(fileresp,"w"))==NULL) {
@@ -1191,14 +1191,27 @@ void  freqsummary(char fileres[], int ag
 	     if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) 
 	       bool=0;
 	 }
-	  if (bool==1) {
-	   for(m=fprev; m<=lprev; m++){
+	 if (bool==1) {
+	   if (boolprev==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];
+	     }
+	   }
+	   else {
+	    for(m=firstpass; m<=lastpass; m++){
+	     k2=anint[m][i]+(mint[m][i]/12.);
+	     if ((k2>=dateprev1) && (k2<=dateprev2)) {
 	     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];
+	     }
+	    }
 	   }
-	 }
+	  }
        }
         if  (cptcovn>0) {
 	 fprintf(ficresp, "\n#********** Variable "); 
@@ -1265,6 +1278,96 @@ 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, double **mint,double **anint,int boolprev, double dateprev1, double dateprev2)
+{  /* Some frequencies */
+ 
+  int i, m, jk, k1, i1, j1, bool, z1,z2,j;
+  double ***freq; /* Frequencies */
+  double *pp;
+  double pos, k2;
+
+  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) {
+	  if (boolprev==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];
+	    }
+	  }
+	  else {
+	    for(m=firstpass; m<=lastpass; m++){
+	      k2=anint[m][i]+(mint[m][i]/12.);
+	      if ((k2>=dateprev1) && (k2<=dateprev2)) {
+		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)
@@ -1783,6 +1886,7 @@ int main()
   char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], fileresf[FILENAMELENGTH];
   char filerest[FILENAMELENGTH];
   char fileregp[FILENAMELENGTH];
+  char popfile[FILENAMELENGTH];
   char path[80],pathc[80],pathcd[80],pathtot[80],model[20];
   int firstobs=1, lastobs=10;
   int sdeb, sfin; /* Status at beginning and end */
@@ -1790,8 +1894,9 @@ 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 mobilav=0, fprevfore=1, lprevfore=1;
+  int mobilav=0, fprev, lprev ,fprevfore=1, lprevfore=1,nforecast,popforecast=0;
   int hstepm, nhstepm;
+  int *popage,boolprev=0;/*boolprev=0 if date and zero if wave*/
 
   double bage, fage, age, agelim, agebase;
   double ftolpl=FTOL;
@@ -1805,7 +1910,9 @@ int main()
   double ***eij, ***vareij;
   double **varpl; /* Variances of prevalence limits by age */
   double *epj, vepp;
-  double kk1;
+  double kk1, kk2;
+  double *popeffectif,*popcount;
+  double dateprev1, dateprev2;
 
   char version[80]="Imach version 64b, May 2001, INED-EUROREVES ";
   char *alph[]={"a","a","b","c","d","e"}, str[4];
@@ -1815,13 +1922,15 @@ int main()
 #include <sys/time.h>
 #include <time.h>
   char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80];
+  char strfprev[10], strlprev[10];
+  char strfprevfore[10], strlprevfore[10];
   /* long total_usecs;
   struct timeval start_time, end_time;
   
   gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
 
 
-  printf("\nIMACH, Version 0.64b");
+  printf("\nIMACH, Version 0.7");
   printf("\nEnter the parameter file name: ");
 
 #ifdef windows
@@ -1877,7 +1986,11 @@ while((c=getc(ficpar))=='#' && c!= EOF){
   }
   ungetc(c,ficpar);
   
-  fscanf(ficpar,"fprevalence=%d lprevalence=%d pop_based=%d\n",&fprev,&lprev,&popbased);
+  fscanf(ficpar,"fprevalence=%s lprevalence=%s pop_based=%d\n",strfprev,strlprev,&popbased);
+  fprintf(ficparo,"fprevalence=%s lprevalence=%s pop_based=%d\n",strfprev,strlprev,popbased);
+ 
+  /* printf("%s %s",strfprev,strlprev);
+     exit(0);*/
  while((c=getc(ficpar))=='#' && c!= EOF){
     ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);
@@ -1886,7 +1999,19 @@ while((c=getc(ficpar))=='#' && c!= EOF){
   }
   ungetc(c,ficpar);
   
-  fscanf(ficpar,"fprevalence=%d lprevalence=%d mob_average=%d\n",&fprevfore,&lprevfore,&mobilav);
+  fscanf(ficpar,"fprevalence=%s lprevalence=%s nforecast=%d mob_average=%d\n",strfprevfore,strlprevfore,&nforecast,&mobilav);
+  fprintf(ficparo,"fprevalence=%s lprevalence=%s nforecast=%d mob_average=%d\n",strfprevfore,strlprevfore,nforecast,mobilav);
+     
+  
+while((c=getc(ficpar))=='#' && c!= EOF){
+    ungetc(c,ficpar);
+    fgets(line, MAXLINE, ficpar);
+    puts(line);
+    fputs(line,ficparo);
+  }
+  ungetc(c,ficpar);
+  
+  fscanf(ficpar,"popforecast=%d popfile=%s\n",&popforecast,popfile);
  
   covar=matrix(0,NCOVMAX,1,n); 
   cptcovn=0; 
@@ -2196,8 +2321,8 @@ printf("Total number of individuals= %d,
     free_imatrix(outcome,1,maxwav+1,1,n);
     free_vector(moisnais,1,n);
     free_vector(annais,1,n);
-    free_matrix(mint,1,maxwav,1,n);
-    free_matrix(anint,1,maxwav,1,n);
+    /* free_matrix(mint,1,maxwav,1,n);
+       free_matrix(anint,1,maxwav,1,n);*/
     free_vector(moisdc,1,n);
     free_vector(andc,1,n);
 
@@ -2229,21 +2354,30 @@ printf("Total number of individuals= %d,
        }
      }
    } 
-
-
-   /*for(i=1; i <=m ;i++){ 
-     for(k=1; k <=cptcovn; k++){
-       printf("i=%d k=%d %d %d",i,k,codtab[i][k], cptcoveff);
-     }
-     printf("\n");
-   }
-   scanf("%d",i);*/
     
    /* 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);
 
-    pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
+    if ((nbocc(strfprev,'/')==1) && (nbocc(strlprev,'/')==1)){
+     boolprev=0;
+     cutv(stra,strb,strfprev,'/');
+     dateprev1=(double)(atoi(strb)+atoi(stra)/12.);
+     cutv(stra,strb,strlprev,'/');
+     dateprev2=(double)(atoi(strb)+atoi(stra)/12.);
+   }
+   
+   else if ((nbocc(strfprev,'/')==0) &&(nbocc(strlprev,'/')==0)){
+     boolprev=1;
+     fprev=atoi(strfprev); lprev=atoi(strlprev);
+   }
+    else {
+      printf("Error in statement lprevalence or fprevalence\n");
+      goto end;
+    }
+   
+  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax, fprev, lprev,mint,anint,boolprev,dateprev1,dateprev2); 
+  
+  pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
     oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
     newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
     savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
@@ -2258,8 +2392,10 @@ printf("Total number of individuals= %d,
     }
     
     /*--------- results files --------------*/
-    fprintf(ficres,"\ntitle=%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(ficres,"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(ficres,"fprevalence=%s lprevalence=%s pop_based=%d\n",strfprev,strlprev,popbased); 
+   fprintf(ficres,"fprevalence=%s lprevalence=%s nforecast=%d mob_average=%d\n",strfprevfore,strlprevfore,nforecast,mobilav);
+
    jk=1;
    fprintf(ficres,"# Parameters\n");
    printf("# Parameters\n");
@@ -2300,7 +2436,7 @@ printf("Total number of individuals= %d,
 	  fprintf(ficres,"\n");
 	}
       }
-      }
+     }
     
     k=1;
     fprintf(ficres,"# Covariance\n");
@@ -2505,17 +2641,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);
@@ -2546,7 +2677,7 @@ chdir(path);
     printf("Problem with %s \n",optionfilehtm);goto end;
   }
 
- fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.64b </font> <hr size=\"2\" color=\"#EC5E5E\"> 
+ fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.7 </font> <hr size=\"2\" color=\"#EC5E5E\"> 
 Titre=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>
 Total number of observations=%d <br>
 Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>
@@ -2561,7 +2692,7 @@ Interval (in months) between two waves:
         - 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>
-        - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>
+        - Prevalences and population 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>");
@@ -2709,12 +2840,29 @@ fclose(fichtm);
   }
   printf("Computing forecasting: result on file '%s' \n", fileresf);
 
-  /* Mobile average */
+ if ((nbocc(strfprevfore,'/')==1) && (nbocc(strlprevfore,'/')==1)){
+     boolprev=0;
+     cutv(stra,strb,strfprevfore,'/');
+     dateprev1=(double)(atoi(strb)+atoi(stra)/12.);
+     cutv(stra,strb,strlprevfore,'/');
+     dateprev2=(double)(atoi(strb)+atoi(stra)/12.);
+   }
+   
+   else if ((nbocc(strfprevfore,'/')==0) &&(nbocc(strlprevfore,'/')==0)){
+     boolprev=1;
+     fprev=atoi(strfprevfore); lprev=atoi(strlprevfore);
+   }
+    else {
+      printf("Error in statement lprevalence or fprevalence\n");
+      goto end;
+    }
 
-  /* 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);*/
+  prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax, fprevfore, lprevfore,mint,anint,boolprev,dateprev1,dateprev2);
+  
+  free_matrix(mint,1,maxwav,1,n);
+  free_matrix(anint,1,maxwav,1,n);
+  free_matrix(agev,1,maxwav,1,imx);
+  /* Mobile average */
 
   if (cptcoveff==0) ncodemax[cptcoveff]=1;
 
@@ -2738,13 +2886,30 @@ fclose(fichtm);
   }
 
   stepsize=(int) (stepm+YEARM-1)/YEARM;
-  if (stepm<=24) stepsize=2;
+  if (stepm<=12) stepsize=1;
 
   agelim=AGESUP;
   hstepm=stepsize*YEARM; /* Every year of age */
-  hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ 
-  hstepm=12;
-   k=0;
+  hstepm=hstepm/stepm; /* Typically 2 years, = 2 years/6 months = 4 */ 
+  
+  if (popforecast==1) {
+    if((ficpop=fopen(popfile,"r"))==NULL)    {
+      printf("Problem with population file : %s\n",popfile);goto end;
+    } 
+    popage=ivector(0,AGESUP);
+    popeffectif=vector(0,AGESUP);
+    popcount=vector(0,AGESUP);
+
+    i=1;   
+    while ((c=fscanf(ficpop,"%d %lf\n",&popage[i],&popcount[i])) != EOF)
+      {
+	i=i+1;
+      }
+    imx=i;
+  
+  for (i=1; i<imx;i++) popeffectif[popage[i]]=popcount[i];
+  }
+
   for(cptcov=1;cptcov<=i1;cptcov++){
     for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
       k=k+1;
@@ -2752,31 +2917,32 @@ fclose(fichtm);
       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);
+      if (popforecast==1)  fprintf(ficresf," [Population]");
 
       for (agedeb=fage; agedeb>=bage; agedeb--){ 
-	fprintf(ficresf,"\n%d %.f %.f 0 ",k,agedeb, agedeb);
+	fprintf(ficresf,"\n%.f %.f 0",agedeb, agedeb);
        if (mobilav==1) {
 	for(j=1; j<=nlstate;j++) 
-	  fprintf(ficresf,"%.5f ",mobaverage[(int)agedeb][j][cptcod]);
+	  fprintf(ficresf," %.3f",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.");
+	  fprintf(ficresf," %.3f",probs[(int)agedeb][j][cptcod]);
+	}  
+
+       for(j=1; j<=ndeath;j++) fprintf(ficresf," 0.00000");
+       if (popforecast==1) fprintf(ficresf," [%.f] ",popeffectif[(int)agedeb]);
       }
-      for (cpt=1; cpt<=NCOVMAX;cpt++)   
+      
+      for (cpt=1; cpt<=nforecast;cpt++) { 
+	fprintf(ficresf,"\n");
       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);*/
+	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;
@@ -2784,24 +2950,38 @@ fclose(fichtm);
 		
 	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);
+	 if (h*hstepm/YEARM*stepm==cpt)
+	    fprintf(ficresf,"\n%.f %.f %.f",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);
+	 for(j=1; j<=nlstate+ndeath;j++) {
+	   kk1=0.;kk2=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 (popforecast==1) kk2=kk1*popeffectif[(int)agedeb];
+	    }
+	   if (h*hstepm/YEARM*stepm==cpt) {
+	     fprintf(ficresf," %.3f", kk1);
+	       if (popforecast==1) fprintf(ficresf," [%.f]", kk2);
+	   }
 	  }
 	}
 	free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+	
+      }
       }
     }
   }
   if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+  if (popforecast==1) {
+    free_ivector(popage,0,AGESUP);
+    free_vector(popeffectif,0,AGESUP);
+    free_vector(popcount,0,AGESUP);
+  }
+  free_imatrix(s,1,maxwav+1,1,n);
+  free_vector(weight,1,n);
   fclose(ficresf);
   /*---------- Health expectancies and variances ------------*/