Some changes in population forecasting

diff --git a/src/imach.c b/src/imach.c
index c5f073450e45686730b7587f6d156e83a5025611..9d7b9d70f9262e7cc31c4af3f0e940c207b5a156 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -94,7 +94,7 @@ int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */
 double jmean; /* Mean space between 2 waves */\r
 double **oldm, **newm, **savm; /* Working pointers to matrices */\r
 double **oldms, **newms, **savms; /* Fixed working pointers to matrices */\r
-FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest,*ficresf;\r
+FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;\r
 FILE *ficgp,*ficresprob,*ficpop;\r
 FILE *ficreseij;\r
   char filerese[FILENAMELENGTH];\r
@@ -2146,6 +2146,261 @@ void movingaverage(double agedeb, double fage,double agemin, double ***mobaverag
     \r
 }\r
 \r
+\r
+/************** Forecasting ******************/\r
+prevforecast(char fileres[], double anproj1,double mproj1,double jproj1,double agemin, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anproj2,double p[], int i2){\r
+  \r
+  int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;\r
+  int *popage;\r
+  double calagedate, agelim, kk1, kk2, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;\r
+  double *popeffectif,*popcount;\r
+  double ***p3mat;\r
+  char fileresf[FILENAMELENGTH];\r
+\r
+ agelim=AGESUP;\r
+calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;\r
+  \r
+  prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);\r
+\r
+ \r
+  strcpy(fileresf,"f"); \r
+  strcat(fileresf,fileres);\r
+  if((ficresf=fopen(fileresf,"w"))==NULL) {\r
+    printf("Problem with forecast resultfile: %s\n", fileresf);\r
+  }\r
+  printf("Computing forecasting: result on file '%s' \n", fileresf);\r
+\r
+  if (cptcoveff==0) ncodemax[cptcoveff]=1;\r
+\r
+  if (mobilav==1) {\r
+    mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);\r
+    movingaverage(agedeb, fage, agemin, mobaverage);\r
+  }\r
+\r
+  stepsize=(int) (stepm+YEARM-1)/YEARM;\r
+  if (stepm<=12) stepsize=1;\r
+  \r
+  agelim=AGESUP;\r
+  \r
+  hstepm=1;\r
+  hstepm=hstepm/stepm; \r
+  yp1=modf(dateintmean,&yp);\r
+  anprojmean=yp;\r
+  yp2=modf((yp1*12),&yp);\r
+  mprojmean=yp;\r
+  yp1=modf((yp2*30.5),&yp);\r
+  jprojmean=yp;\r
+  if(jprojmean==0) jprojmean=1;\r
+  if(mprojmean==0) jprojmean=1;\r
+  \r
+  fprintf(ficresf,"# Estimated date of observed prevalence: %.lf/%.lf/%.lf ",jprojmean,mprojmean,anprojmean); \r
+  \r
+  for(cptcov=1;cptcov<=i2;cptcov++){\r
+    for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){\r
+      k=k+1;\r
+      fprintf(ficresf,"\n#******");\r
+      for(j=1;j<=cptcoveff;j++) {\r
+       fprintf(ficresf," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);\r
+      }\r
+      fprintf(ficresf,"******\n");\r
+      fprintf(ficresf,"# StartingAge FinalAge");\r
+      for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);\r
+      \r
+      \r
+      for (cpt=0; cpt<=(anproj2-anproj1);cpt++) { \r
+       fprintf(ficresf,"\n");\r
+       fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);   \r
+       \r
+       for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){ \r
+         nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); \r
+         nhstepm = nhstepm/hstepm; \r
+         \r
+         p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);\r
+         oldm=oldms;savm=savms;\r
+         hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  \r
+       \r
+         for (h=0; h<=nhstepm; h++){\r
+           if (h==(int) (calagedate+YEARM*cpt)) {\r
+             fprintf(ficresf,"\n %.f ",agedeb+h*hstepm/YEARM*stepm);\r
+           } \r
+           for(j=1; j<=nlstate+ndeath;j++) {\r
+             kk1=0.;kk2=0;\r
+             for(i=1; i<=nlstate;i++) {              \r
+               if (mobilav==1) \r
+                 kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod];\r
+               else {\r
+                 kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod];\r
+               }\r
+               \r
+             }\r
+             if (h==(int)(calagedate+12*cpt)){\r
+               fprintf(ficresf," %.3f", kk1);\r
+                       \r
+             }\r
+           }\r
+         }\r
+         free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);\r
+       }\r
+      }\r
+    }\r
+  }\r
+       \r
+  if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);\r
+\r
+  fclose(ficresf);\r
+}\r
+/************** Forecasting ******************/\r
+populforecast(char fileres[], double anproj1,double mproj1,double jproj1,double agemin, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anproj2,double p[], int i2){\r
+  \r
+  int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;\r
+  int *popage;\r
+  double calagedate, agelim, kk1, kk2, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;\r
+  double *popeffectif,*popcount;\r
+  double ***p3mat,***tabpop,***tabpopprev;\r
+  char filerespop[FILENAMELENGTH];\r
+\r
+tabpop= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);\r
+tabpopprev= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);\r
+ agelim=AGESUP;\r
+calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;\r
+  \r
+  prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);\r
+\r
+ \r
+  strcpy(filerespop,"pop"); \r
+  strcat(filerespop,fileres);\r
+  if((ficrespop=fopen(filerespop,"w"))==NULL) {\r
+    printf("Problem with forecast resultfile: %s\n", filerespop);\r
+  }\r
+  printf("Computing forecasting: result on file '%s' \n", filerespop);\r
+\r
+  if (cptcoveff==0) ncodemax[cptcoveff]=1;\r
+\r
+  if (mobilav==1) {\r
+    mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);\r
+    movingaverage(agedeb, fage, agemin, mobaverage);\r
+  }\r
+\r
+  stepsize=(int) (stepm+YEARM-1)/YEARM;\r
+  if (stepm<=12) stepsize=1;\r
+  \r
+  agelim=AGESUP;\r
+  \r
+  hstepm=1;\r
+  hstepm=hstepm/stepm; \r
+  \r
+  if (popforecast==1) {\r
+    if((ficpop=fopen(popfile,"r"))==NULL) {\r
+      printf("Problem with population file : %s\n",popfile);exit(0);\r
+    } \r
+    popage=ivector(0,AGESUP);\r
+    popeffectif=vector(0,AGESUP);\r
+    popcount=vector(0,AGESUP);\r
+    \r
+    i=1;   \r
+    while ((c=fscanf(ficpop,"%d %lf\n",&popage[i],&popcount[i])) != EOF) i=i+1;\r
+   \r
+    imx=i;\r
+    for (i=1; i<imx;i++) popeffectif[popage[i]]=popcount[i];\r
+  }\r
+\r
+  for(cptcov=1;cptcov<=i2;cptcov++){\r
+   for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){\r
+      k=k+1;\r
+      fprintf(ficrespop,"\n#******");\r
+      for(j=1;j<=cptcoveff;j++) {\r
+       fprintf(ficrespop," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);\r
+      }\r
+      fprintf(ficrespop,"******\n");\r
+      fprintf(ficrespop,"# StartingAge FinalAge");\r
+      for(j=1; j<=nlstate+ndeath;j++) fprintf(ficrespop," P.%d",j);\r
+      if (popforecast==1)  fprintf(ficrespop," [Population]");\r
+      \r
+      for (cpt=0; cpt<=0;cpt++) { \r
+       fprintf(ficrespop,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);   \r
+       \r
+       for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){ \r
+         nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); \r
+         nhstepm = nhstepm/hstepm; \r
+         \r
+         p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);\r
+         oldm=oldms;savm=savms;\r
+         hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  \r
+       \r
+         for (h=0; h<=nhstepm; h++){\r
+           if (h==(int) (calagedate+YEARM*cpt)) {\r
+             fprintf(ficrespop,"\n %.f ",agedeb+h*hstepm/YEARM*stepm);\r
+           } \r
+           for(j=1; j<=nlstate+ndeath;j++) {\r
+             kk1=0.;kk2=0;\r
+             for(i=1; i<=nlstate;i++) {              \r
+               if (mobilav==1) \r
+                 kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod];\r
+               else {\r
+                 kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod];\r
+               }\r
+             }\r
+             if (h==(int)(calagedate+12*cpt)){\r
+               tabpop[(int)(agedeb)][j][cptcod]=kk1;\r
+                 /*fprintf(ficrespop," %.3f", kk1);\r
+                   if (popforecast==1) fprintf(ficrespop," [%.f]", kk1*popeffectif[(int)agedeb+1]);*/\r
+             }\r
+           }\r
+           for(i=1; i<=nlstate;i++){\r
+             kk1=0.;\r
+               for(j=1; j<=nlstate;j++){\r
+                 kk1= kk1+Tabpop[(int)(agedeb)][j][cptcod]; \r
+               }\r
+                 tabpopprev[(int)(agedeb)][i][cptcod]=tabpop[(int)(agedeb)][i][cptcod]/kk1*popeffectif[(int)(agedeb+(calagedate+12*cpt)*hstepm/YEARM*stepm-1)];\r
+           }\r
+\r
+if (h==(int)(calagedate+12*cpt)) for(j=1; j<=nlstate;j++) fprintf(ficrespop," %.3f",tabpopprev[(int)(agedeb+1)][j][cptcod]);\r
+         }\r
+         free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);\r
+       }\r
+      }\r
+ \r
+  /******/\r
+\r
+      for (cpt=1; cpt<=4;cpt++) { \r
+       fprintf(ficrespop,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);   \r
+       for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){ \r
+         nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); \r
+         nhstepm = nhstepm/hstepm; \r
+         \r
+         p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);\r
+         oldm=oldms;savm=savms;\r
+         hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  \r
+         for (h=0; h<=nhstepm; h++){\r
+           if (h==(int) (calagedate+YEARM*cpt)) {\r
+             fprintf(ficresf,"\n %.f ",agedeb+h*hstepm/YEARM*stepm);\r
+           } \r
+           for(j=1; j<=nlstate+ndeath;j++) {\r
+             kk1=0.;kk2=0;\r
+             for(i=1; i<=nlstate;i++) {              \r
+               kk1=kk1+p3mat[i][j][h]*tabpopprev[(int)agedeb+1][i][cptcod];    \r
+             }\r
+             if (h==(int)(calagedate+12*cpt)) fprintf(ficresf," %.3f", kk1);   \r
+           }\r
+         }\r
+         free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);\r
+       }\r
+      }\r
+   } \r
+  }\r
+ \r
+  if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);\r
+\r
+  if (popforecast==1) {\r
+    free_ivector(popage,0,AGESUP);\r
+    free_vector(popeffectif,0,AGESUP);\r
+    free_vector(popcount,0,AGESUP);\r
+  }\r
+  free_ma3x(tabpop,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);\r
+  free_ma3x(tabpopprev,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);\r
+  fclose(ficrespop);\r
+}\r
+\r
 /***********************************************/\r
 /**************** Main Program *****************/\r
 /***********************************************/\r
@@ -2168,7 +2423,7 @@ int main(int argc, char *argv[])
   char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH];\r
   char optionfilext[10], optionfilefiname[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH], optionfilegnuplot[FILENAMELENGTH], plotcmd[FILENAMELENGTH];\r
   \r
-  char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], fileresf[FILENAMELENGTH];\r
+  char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH];\r
 \r
   char filerest[FILENAMELENGTH];\r
   char fileregp[FILENAMELENGTH];\r
@@ -2182,7 +2437,6 @@ int main(int argc, char *argv[])
   int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab; \r
   int mobilav=0,popforecast=0;\r
   int hstepm, nhstepm;\r
-  int *popage;/*boolprev=0 if date and zero if wave*/\r
   double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2;\r
 \r
   double bage, fage, age, agelim, agebase;\r
@@ -2198,9 +2452,8 @@ int main(int argc, char *argv[])
   double **varpl; /* Variances of prevalence limits by age */\r
   double *epj, vepp;\r
   double kk1, kk2;\r
-  double *popeffectif,*popcount;\r
-  double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,jprojmean,mprojmean,anprojmean, calagedate;\r
-  double yp,yp1,yp2;\r
+  double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;\r
+  \r
 \r
   char version[80]="Imach version 0.7, February 2002, INED-EUROREVES ";\r
   char *alph[]={"a","a","b","c","d","e"}, str[4];\r
@@ -2902,139 +3155,20 @@ fprintf(ficres,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf fina
 \r
   fclose(ficrespij);\r
 \r
-  if(stepm == 1) {\r
-  /*---------- Forecasting ------------------*/\r
-  calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;\r
-\r
-  \r
-  \r
-  prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);\r
-\r
- \r
-  strcpy(fileresf,"f"); \r
-  strcat(fileresf,fileres);\r
-  if((ficresf=fopen(fileresf,"w"))==NULL) {\r
-    printf("Problem with forecast resultfile: %s\n", fileresf);goto end;\r
-  }\r
-  printf("Computing forecasting: result on file '%s' \n", fileresf);\r
-  \r
-  free_matrix(mint,1,maxwav,1,n);\r
-  free_matrix(anint,1,maxwav,1,n);\r
-\r
-  if (cptcoveff==0) ncodemax[cptcoveff]=1;\r
-\r
-  if (mobilav==1) {\r
-    mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);\r
-    movingaverage(agedeb, fage, agemin, mobaverage);\r
-  }\r
-\r
-  stepsize=(int) (stepm+YEARM-1)/YEARM;\r
-  if (stepm<=12) stepsize=1;\r
-  \r
-  agelim=AGESUP;\r
-  \r
-  hstepm=1;\r
-  hstepm=hstepm/stepm; \r
-  yp1=modf(dateintmean,&yp);\r
-  anprojmean=yp;\r
-  yp2=modf((yp1*12),&yp);\r
-  mprojmean=yp;\r
-  yp1=modf((yp2*30.5),&yp);\r
-  jprojmean=yp;\r
-  if(jprojmean==0) jprojmean=1;\r
-  if(mprojmean==0) jprojmean=1;\r
-  \r
-  fprintf(ficresf,"# Estimated date of observed prevalence: %.lf/%.lf/%.lf ",jprojmean,mprojmean,anprojmean); \r
-  \r
-  if (popforecast==1) {\r
-    if((ficpop=fopen(popfile,"r"))==NULL)    {\r
-      printf("Problem with population file : %s\n",popfile);goto end;\r
-    } \r
-    popage=ivector(0,AGESUP);\r
-    popeffectif=vector(0,AGESUP);\r
-    popcount=vector(0,AGESUP);\r
-    \r
-    i=1;   \r
-    while ((c=fscanf(ficpop,"%d %lf\n",&popage[i],&popcount[i])) != EOF)\r
-      {\r
-       i=i+1;\r
-      }\r
-    imx=i;\r
-    \r
-    for (i=1; i<imx;i++) popeffectif[popage[i]]=popcount[i];\r
-  }\r
-  \r
-  for(cptcov=1;cptcov<=i1;cptcov++){\r
-    for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){\r
-      k=k+1;\r
-      fprintf(ficresf,"\n#******");\r
-      for(j=1;j<=cptcoveff;j++) {\r
-       fprintf(ficresf," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);\r
-      }\r
-      fprintf(ficresf,"******\n");\r
-      fprintf(ficresf,"# StartingAge FinalAge");\r
-      for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);\r
-      if (popforecast==1)  fprintf(ficresf," [Population]");\r
-      \r
-      for (cpt=0; cpt<=(anproj2-anproj1);cpt++) { \r
-       fprintf(ficresf,"\n");\r
-       fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);   \r
-       \r
-       for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){ \r
-         nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); \r
-         nhstepm = nhstepm/hstepm; \r
-         \r
-         p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);\r
-         oldm=oldms;savm=savms;\r
-         hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  \r
-       \r
-         for (h=0; h<=nhstepm; h++){\r
-           if (h==(int) (calagedate+YEARM*cpt)) {\r
-             fprintf(ficresf,"\n %.f ",agedeb+h*hstepm/YEARM*stepm);\r
-           } \r
-           for(j=1; j<=nlstate+ndeath;j++) {\r
-             kk1=0.;kk2=0;\r
-             for(i=1; i<=nlstate;i++) {              \r
-               if (mobilav==1) \r
-                 kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod];\r
-               else {\r
-                 kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod];\r
-               }\r
-               if (popforecast==1) kk2=kk1*popeffectif[(int)agedeb];\r
-             }\r
-             if (h==(int)(calagedate+12*cpt)){\r
-               fprintf(ficresf," %.3f", kk1);\r
-            \r
-               if (popforecast==1) fprintf(ficresf," [%.f]", kk2);\r
-             }\r
-           }\r
-         }\r
-         free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);\r
-       }\r
-      }\r
-    }\r
-  }\r
-       \r
-  if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);\r
-\r
-  if (popforecast==1) {\r
-    free_ivector(popage,0,AGESUP);\r
-    free_vector(popeffectif,0,AGESUP);\r
-    free_vector(popcount,0,AGESUP);\r
-  }\r
-  \r
-  free_imatrix(s,1,maxwav+1,1,n);\r
-  free_vector(weight,1,n);\r
-  fclose(ficresf);\r
-  }\r
 \r
-  /* End forecasting */\r
+  /*---------- Forecasting ------------------*/\r
+  if(stepm == 1) {\r
+    prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1);\r
+populforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1);\r
+    free_matrix(mint,1,maxwav,1,n);\r
+    free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n);\r
+    free_vector(weight,1,n);}\r
   else{\r
     erreur=108;\r
     printf("Error %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d\n", erreur, stepm);\r
   }\r
   \r
- \r
+\r
   /*---------- Health expectancies and variances ------------*/\r
 \r
   strcpy(filerest,"t");\r
@@ -3119,20 +3253,16 @@ fprintf(ficres,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf fina
       }\r
     }\r
   }\r
-       \r
-       \r
-\r
 \r
- fclose(ficreseij);\r
- fclose(ficresvij);\r
+  fclose(ficreseij);\r
+  fclose(ficresvij);\r
   fclose(ficrest);\r
   fclose(ficpar);\r
   free_vector(epj,1,nlstate+1);\r
-  /*  scanf("%d ",i); */\r
   \r
   /*------- Variance limit prevalence------*/   \r
 \r
-strcpy(fileresvpl,"vpl");\r
+  strcpy(fileresvpl,"vpl");\r
   strcat(fileresvpl,fileres);\r
   if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {\r
     printf("Problem with variance prev lim resultfile: %s\n", fileresvpl);\r
@@ -3140,19 +3270,19 @@ strcpy(fileresvpl,"vpl");
   }\r
   printf("Computing Variance-covariance of Prevalence limit: file '%s' \n", fileresvpl);\r
 \r
- k=0;\r
- for(cptcov=1;cptcov<=i1;cptcov++){\r
-   for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){\r
-     k=k+1;\r
-     fprintf(ficresvpl,"\n#****** ");\r
-     for(j=1;j<=cptcoveff;j++) \r
-       fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);\r
-     fprintf(ficresvpl,"******\n");\r
-     \r
-     varpl=matrix(1,nlstate,(int) bage, (int) fage);\r
-     oldm=oldms;savm=savms;\r
+  k=0;\r
+  for(cptcov=1;cptcov<=i1;cptcov++){\r
+    for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){\r
+      k=k+1;\r
+      fprintf(ficresvpl,"\n#****** ");\r
+      for(j=1;j<=cptcoveff;j++) \r
+       fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);\r
+      fprintf(ficresvpl,"******\n");\r
+      \r
+      varpl=matrix(1,nlstate,(int) bage, (int) fage);\r
+      oldm=oldms;savm=savms;\r
      varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k);\r
-   }\r
+    }\r
  }\r
 \r
   fclose(ficresvpl);\r