mobilav can be 3, 5 or 7 instead of only 1 (which defaults to 5).

diff --git a/src/imach.c b/src/imach.c
index 2d8cd1e22323a4321f08c634cb9696f3b2a9a129..78fa65f2074aa95ca97df728abe656c7c50b3673 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -39,7 +39,7 @@
   hPijx.
 
   Also this programme outputs the covariance matrix of the parameters but also
-  of the life expectancies. It also computes the prevalence limits. 
+  of the life expectancies. It also computes the stable prevalence. 
   
   Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr).
            Institut national d'études démographiques, Paris.
@@ -83,7 +83,7 @@
 #define ODIRSEPARATOR '\\'
 #endif
 
-char version[80]="Imach version 0.8j, July 2002, INED-EUROREVES ";
+char version[80]="Imach version 0.8k, July 2002, INED-EUROREVES ";
 int erreur; /* Error number */
 int nvar;
 int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;
@@ -701,12 +701,12 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret
        printf("\n");
        fprintf(ficlog,"\n");
 #endif
-      } 
+      }
     } 
   } 
 } 
 
-/**** Prevalence limit ****************/
+/**** Prevalence limit (stable prevalence)  ****************/
 
 double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij)
 {
@@ -1827,13 +1827,16 @@ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double
     strcpy(digitp,"-populbased-");
   else
     strcpy(digitp,"-stablbased-");
-  if(mobilav==1)
+  if(mobilav!=0)
     strcat(digitp,"mobilav-");
   else
     strcat(digitp,"nomobil-");
-  if (mobilav==1) {
+  if (mobilav!=0) {
     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
-    movingaverage(probs, bage, fage, mobaverage);
+    if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){
+      fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
+      printf(" Error in movingaverage mobilav=%d\n",mobilav);
+    }
   }
 
   strcpy(fileresprobmorprev,"prmorprev"); 
@@ -1928,10 +1931,10 @@ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
 
       if (popbased==1) {
-       if(mobilav !=1){
+       if(mobilav ==0){
          for(i=1; i<=nlstate;i++)
            prlim[i][i]=probs[(int)age][i][ij];
-       }else{ /* mobilav=1 */ 
+       }else{ /* mobilav */ 
          for(i=1; i<=nlstate;i++)
            prlim[i][i]=mobaverage[(int)age][i][ij];
        }
@@ -1956,10 +1959,10 @@ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
  
       if (popbased==1) {
-       if(mobilav !=1){
+       if(mobilav ==0){
          for(i=1; i<=nlstate;i++)
            prlim[i][i]=probs[(int)age][i][ij];
-       }else{ /* mobilav=1 */ 
+       }else{ /* mobilav */ 
          for(i=1; i<=nlstate;i++)
            prlim[i][i]=mobaverage[(int)age][i][ij];
        }
@@ -2025,10 +2028,10 @@ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double
     prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij);
  
     if (popbased==1) {
-      if(mobilav !=1){
+      if(mobilav ==0){
        for(i=1; i<=nlstate;i++)
          prlim[i][i]=probs[(int)age][i][ij];
-      }else{ /* mobilav=1 */ 
+      }else{ /* mobilav */ 
        for(i=1; i<=nlstate;i++)
          prlim[i][i]=mobaverage[(int)age][i][ij];
       }
@@ -2106,7 +2109,7 @@ void varprevlim(char fileres[], double **varpl, double **matcov, double x[], dou
   double age,agelim;
   int theta;
    
-  fprintf(ficresvpl,"# Standard deviation of prevalence's limit\n");
+  fprintf(ficresvpl,"# Standard deviation of stable prevalences \n");
   fprintf(ficresvpl,"# Age");
   for(i=1; i<=nlstate;i++)
       fprintf(ficresvpl," %1d-%1d",i,i);
@@ -2599,9 +2602,9 @@ void printinggnuplot(char fileres[], double ageminpar, double agemaxpar, double
     fprintf(ficlog,"Problem with file %s",optionfilegnuplot);
   }
 
-#ifdef windows
+  /*#ifdef windows */
     fprintf(ficgp,"cd \"%s\" \n",pathc);
-#endif
+    /*#endif */
 m=pow(2,cptcoveff);
   
  /* 1eme*/
@@ -2614,7 +2617,7 @@ m=pow(2,cptcoveff);
        if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
        else fprintf(ficgp," \%%*lf (\%%*lf)");
      }
-     fprintf(ficgp,"\" t\"Stationary prevalence\" w l 0,\"vpl%s\" every :::%d::%d u 1:($2+2*$3) \"\%%lf",fileres,k1-1,k1-1);
+     fprintf(ficgp,"\" t\"Stable prevalence\" w l 0,\"vpl%s\" every :::%d::%d u 1:($2+2*$3) \"\%%lf",fileres,k1-1,k1-1);
      for (i=1; i<= nlstate ; i ++) {
        if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
        else fprintf(ficgp," \%%*lf (\%%*lf)");
@@ -2768,27 +2771,39 @@ m=pow(2,cptcoveff);
 
 
 /*************** Moving average **************/
-void movingaverage(double ***probs, double bage,double fage, double ***mobaverage){
+int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav){
 
   int i, cpt, cptcod;
+  int mobilavrange, mob;
   double age;
-  for (age=bage; age<=fage; age++)
-    for (i=1; i<=nlstate;i++)
-      for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)
-       mobaverage[(int)age][i][cptcod]=0.;
-  
-  for (age=bage+4; age<=fage; age++){
-    for (i=1; i<=nlstate;i++){
-      for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
-       for (cpt=0;cpt<=4;cpt++){
-         mobaverage[(int)age-2][i][cptcod]=mobaverage[(int)age-2][i][cptcod]+probs[(int)age-cpt][i][cptcod];
+  if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){
+    if(mobilav==1) mobilavrange=5; /* default */
+    else mobilavrange=mobilav;
+    for (age=bage; age<=fage; age++)
+      for (i=1; i<=nlstate;i++)
+       for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)
+         mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod];
+    /* We keep the original values on the extreme ages bage, fage and for 
+       fage+1 and bage-1 we use a 3 terms moving average; for fage+2 bage+2
+       we use a 5 terms etc. until the borders are no more concerned. 
+    */ 
+    for (mob=3;mob <=mobilavrange;mob=mob+2){
+      for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){
+       for (i=1; i<=nlstate;i++){
+         for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
+           mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod];
+             for (cpt=1;cpt<=(mob-1)/2;cpt++){
+               mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod];
+               mobaverage[(int)age][i][cptcod] +=probs[(int)age+cpt][i][cptcod];
+             }
+           mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/mob;
+         }
        }
-       mobaverage[(int)age-2][i][cptcod]=mobaverage[(int)age-2][i][cptcod]/5;
-      }
-    }
-  }
-  
-}
+      }/* end age */
+    }/* end mob */
+  }else return -1;
+  return 0;
+}/* End movingaverage */
 
 
 /************** Forecasting ******************/
@@ -2818,9 +2833,12 @@ calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;
 
   if (cptcoveff==0) ncodemax[cptcoveff]=1;
 
-  if (mobilav==1) {
+  if (mobilav!=0) {
     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
-    movingaverage(probs, ageminpar,fage, mobaverage);
+    if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){
+      fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
+      printf(" Error in movingaverage mobilav=%d\n",mobilav);
+    }
   }
 
   stepsize=(int) (stepm+YEARM-1)/YEARM;
@@ -2891,7 +2909,7 @@ calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;
     }
   }
        
-  if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+  if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
 
   fclose(ficresf);
 }
@@ -2924,9 +2942,12 @@ populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double
 
   if (cptcoveff==0) ncodemax[cptcoveff]=1;
 
-  if (mobilav==1) {
+  if (mobilav!=0) {
     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
-    movingaverage(probs, ageminpar, fage, mobaverage);
+    if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){
+      fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
+      printf(" Error in movingaverage mobilav=%d\n",mobilav);
+    }
   }
 
   stepsize=(int) (stepm+YEARM-1)/YEARM;
@@ -3039,7 +3060,7 @@ populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double
    } 
   }
  
-  if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+  if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
 
   if (popforecast==1) {
     free_ivector(popage,0,AGESUP);
@@ -3750,12 +3771,17 @@ while((c=getc(ficpar))=='#' && c!= EOF){
 
  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
 /*------------ gnuplot -------------*/
-  strcpy(optionfilegnuplot,optionfilefiname);
-  strcat(optionfilegnuplot,".gp");
-  if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
-    printf("Problem with file %s",optionfilegnuplot);
-  }
-  fclose(ficgp);
+ strcpy(optionfilegnuplot,optionfilefiname);
+ strcat(optionfilegnuplot,".gp");
+ if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
+   printf("Problem with file %s",optionfilegnuplot);
+ }
+ else{
+   fprintf(ficgp,"\n# %s\n", version); 
+   fprintf(ficgp,"# %s\n", optionfilegnuplot); 
+   fprintf(ficgp,"set missing 'NaNq'\n");
+}
+ fclose(ficgp);
  printinggnuplot(fileres, ageminpar,agemaxpar,fage, pathc,p);
 /*--------- index.htm --------*/
 
@@ -3792,17 +3818,17 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
  fclose(ficres);
 
 
-  /*--------------- Prevalence limit --------------*/
+  /*--------------- Prevalence limit  (stable prevalence) --------------*/
   
   strcpy(filerespl,"pl");
   strcat(filerespl,fileres);
   if((ficrespl=fopen(filerespl,"w"))==NULL) {
-    printf("Problem with Prev limit resultfile: %s\n", filerespl);goto end;
-    fprintf(ficlog,"Problem with Prev limit resultfile: %s\n", filerespl);goto end;
+    printf("Problem with stable prevalence resultfile: %s\n", filerespl);goto end;
+    fprintf(ficlog,"Problem with stable prevalence resultfile: %s\n", filerespl);goto end;
   }
-  printf("Computing prevalence limit: result on file '%s' \n", filerespl);
-  fprintf(ficlog,"Computing prevalence limit: result on file '%s' \n", filerespl);
-  fprintf(ficrespl,"#Prevalence limit\n");
+  printf("Computing stable prevalence: result on file '%s' \n", filerespl);
+  fprintf(ficlog,"Computing stable prevalence: result on file '%s' \n", filerespl);
+  fprintf(ficrespl,"#Stable prevalence \n");
   fprintf(ficrespl,"#Age ");
   for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
   fprintf(ficrespl,"\n");
@@ -3950,9 +3976,12 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
   fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
   calagedate=-1;
   prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
-  if (mobilav==1) {
+  if (mobilav!=0) {
     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
-    movingaverage(probs, ageminpar, fage, mobaverage);
+    if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){
+      fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
+      printf(" Error in movingaverage mobilav=%d\n",mobilav);
+    }
   }
 
   k=0;
@@ -3994,10 +4023,10 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
       for(age=bage; age <=fage ;age++){
        prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
        if (popbased==1) {
-         if(mobilav !=1){
+         if(mobilav ==0){
            for(i=1; i<=nlstate;i++)
              prlim[i][i]=probs[(int)age][i][k];
-         }else{ /* mobilav=1 */ 
+         }else{ /* mobilav */ 
            for(i=1; i<=nlstate;i++)
              prlim[i][i]=mobaverage[(int)age][i][k];
          }
@@ -4032,15 +4061,15 @@ free_matrix(mint,1,maxwav,1,n);
   fclose(ficpar);
   free_vector(epj,1,nlstate+1);
   
-  /*------- Variance limit prevalence------*/   
+  /*------- Variance of stable prevalence------*/   
 
   strcpy(fileresvpl,"vpl");
   strcat(fileresvpl,fileres);
   if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
-    printf("Problem with variance prev lim resultfile: %s\n", fileresvpl);
+    printf("Problem with variance of stable prevalence  resultfile: %s\n", fileresvpl);
     exit(0);
   }
-  printf("Computing Variance-covariance of Prevalence limit: file '%s' \n", fileresvpl);
+  printf("Computing Variance-covariance of stable prevalence: file '%s' \n", fileresvpl);
 
   k=0;
   for(cptcov=1;cptcov<=i1;cptcov++){
@@ -4075,7 +4104,7 @@ free_matrix(mint,1,maxwav,1,n);
   free_vector(delti,1,npar);
   free_matrix(agev,1,maxwav,1,imx);
   free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
-  if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
+  if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
 
   fprintf(fichtm,"\n</body>");
   fclose(fichtm);
@@ -4109,9 +4138,10 @@ free_matrix(mint,1,maxwav,1,n);
  strcpy(plotcmd,GNUPLOTPROGRAM);
  strcat(plotcmd," ");
  strcat(plotcmd,optionfilegnuplot);
+ printf("Starting: %s\n",plotcmd);fflush(stdout);
  system(plotcmd);
 
-#ifdef windows
+ /*#ifdef windows*/
   while (z[0] != 'q') {
     /* chdir(path); */
     printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: ");
@@ -4121,7 +4151,7 @@ free_matrix(mint,1,maxwav,1,n);
     else if (z[0] == 'g') system(plotcmd);
     else if (z[0] == 'q') exit(0);
   }
-#endif 
+  /*#endif */
 }