Version 0.7

hPijx are now computed for h=1 year instead of 2 years

diff --git a/src/imach.c b/src/imach.c
index e9b172fc5967a6a952316d201c41942e432df9b6..97f01bd3bbcad712842d0324e5ffb1b82b55d594 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -85,7 +85,7 @@ double jmean; /* Mean space between 2 waves */
 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 *ficgp, *fichtm,*ficresprob;\r
+FILE *ficgp, *fichtm,*ficresprob,*ficpop;\r
 FILE *ficreseij;\r
   char filerese[FILENAMELENGTH];\r
  FILE  *ficresvij;\r
@@ -1860,6 +1860,7 @@ int main()
   char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], fileresf[FILENAMELENGTH];\r
   char filerest[FILENAMELENGTH];\r
   char fileregp[FILENAMELENGTH];\r
+  char popfile[FILENAMELENGTH];\r
   char path[80],pathc[80],pathcd[80],pathtot[80],model[20];\r
   int firstobs=1, lastobs=10;\r
   int sdeb, sfin; /* Status at beginning and end */\r
@@ -1867,8 +1868,9 @@ int main()
   int ju,jl, mi;\r
   int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;\r
   int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab; \r
-  int mobilav=0, fprev, lprev ,fprevfore=1, lprevfore=1,nforecast;\r
+  int mobilav=0, fprev, lprev ,fprevfore=1, lprevfore=1,nforecast,popforecast=0;\r
   int hstepm, nhstepm;\r
+  int *popage;\r
 \r
   double bage, fage, age, agelim, agebase;\r
   double ftolpl=FTOL;\r
@@ -1882,7 +1884,8 @@ int main()
   double ***eij, ***vareij;\r
   double **varpl; /* Variances of prevalence limits by age */\r
   double *epj, vepp;\r
-  double kk1;\r
+  double kk1, kk2;\r
+  double *popeffectif,*popcount;\r
 \r
   char version[80]="Imach version 64b, May 2001, INED-EUROREVES ";\r
   char *alph[]={"a","a","b","c","d","e"}, str[4];\r
@@ -1898,7 +1901,7 @@ int main()
   gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */\r
 \r
 \r
-  printf("\nIMACH, Version 0.64b");\r
+  printf("\nIMACH, Version 0.7");\r
   printf("\nEnter the parameter file name: ");\r
 \r
 #ifdef windows\r
@@ -1955,6 +1958,8 @@ while((c=getc(ficpar))=='#' && c!= EOF){
   ungetc(c,ficpar);\r
   \r
   fscanf(ficpar,"fprevalence=%d lprevalence=%d pop_based=%d\n",&fprev,&lprev,&popbased);\r
+  fprintf(ficparo,"fprevalence=%d lprevalence=%d pop_based=%d\n",fprev,lprev,popbased);\r
+ \r
  while((c=getc(ficpar))=='#' && c!= EOF){\r
     ungetc(c,ficpar);\r
     fgets(line, MAXLINE, ficpar);\r
@@ -1964,6 +1969,18 @@ while((c=getc(ficpar))=='#' && c!= EOF){
   ungetc(c,ficpar);\r
   \r
   fscanf(ficpar,"fprevalence=%d lprevalence=%d nforecast=%d mob_average=%d\n",&fprevfore,&lprevfore,&nforecast,&mobilav);\r
+  fprintf(ficparo,"fprevalence=%d lprevalence=%d nforecast=%d mob_average=%d\n",fprevfore,lprevfore,nforecast,mobilav);\r
+     \r
+  \r
+while((c=getc(ficpar))=='#' && c!= EOF){\r
+    ungetc(c,ficpar);\r
+    fgets(line, MAXLINE, ficpar);\r
+    puts(line);\r
+    fputs(line,ficparo);\r
+  }\r
+  ungetc(c,ficpar);\r
+  \r
+  fscanf(ficpar,"popforecast=%d popfile=%s\n",&popforecast,popfile);\r
  \r
   covar=matrix(0,NCOVMAX,1,n); \r
   cptcovn=0; \r
@@ -2335,8 +2352,10 @@ printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx,
     }\r
     \r
     /*--------- results files --------------*/\r
-    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);\r
-    \r
+    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);\r
+   fprintf(ficres,"fprevalence=%d lprevalence=%d pop_based=%d\n",fprev,lprev,popbased); \r
+   fprintf(ficres,"fprevalence=%d lprevalence=%d nforecast=%d mob_average=%d\n",fprevfore,lprevfore,nforecast,mobilav);\r
+\r
    jk=1;\r
    fprintf(ficres,"# Parameters\n");\r
    printf("# Parameters\n");\r
@@ -2618,7 +2637,7 @@ chdir(path);
     printf("Problem with %s \n",optionfilehtm);goto end;\r
   }\r
 \r
- fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.64b </font> <hr size=\"2\" color=\"#EC5E5E\"> \r
+ fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.7 </font> <hr size=\"2\" color=\"#EC5E5E\"> \r
 Titre=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\r
 Total number of observations=%d <br>\r
 Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\r
@@ -2813,8 +2832,25 @@ fclose(fichtm);
   agelim=AGESUP;\r
   hstepm=stepsize*YEARM; /* Every year of age */\r
   hstepm=hstepm/stepm; /* Typically 2 years, = 2 years/6 months = 4 */ \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
-   k=0;\r
   for(cptcov=1;cptcov<=i1;cptcov++){\r
     for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){\r
       k=k+1;\r
@@ -2825,22 +2861,26 @@ fclose(fichtm);
       fprintf(ficresf,"******\n");\r
       fprintf(ficresf,"# StartingAge FinalAge Horizon(in years)");\r
       for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);\r
+      if (popforecast==1)  fprintf(ficresf," [Population]");\r
 \r
       for (agedeb=fage; agedeb>=bage; agedeb--){ \r
-       fprintf(ficresf,"\n%d %.f %.f 0 ",k,agedeb, agedeb);\r
+       fprintf(ficresf,"\n%.f %.f 0",agedeb, agedeb);\r
        if (mobilav==1) {\r
        for(j=1; j<=nlstate;j++) \r
-         fprintf(ficresf," %.5f ",mobaverage[(int)agedeb][j][cptcod]);\r
+         fprintf(ficresf," %.3f",mobaverage[(int)agedeb][j][cptcod]);\r
        }\r
        else {\r
          for(j=1; j<=nlstate;j++) \r
-         fprintf(ficresf," %.5f ",probs[(int)agedeb][j][cptcod]);\r
+         fprintf(ficresf," %.3f",probs[(int)agedeb][j][cptcod]);\r
        }  \r
-      for(j=1; j<=ndeath;j++) fprintf(ficresf," 0.00000");\r
+\r
+       for(j=1; j<=ndeath;j++) fprintf(ficresf," 0.00000");\r
+       if (popforecast==1) fprintf(ficresf," [%.f] ",popeffectif[(int)agedeb]);\r
       }\r
-      for (cpt=1; cpt<=nforecast;cpt++)   \r
+      \r
+      for (cpt=1; cpt<=nforecast;cpt++) { \r
+       fprintf(ficresf,"\n");\r
       for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */\r
-       \r
        nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); \r
        nhstepm = nhstepm/hstepm; \r
        /*printf("agedeb=%.lf stepm=%d hstepm=%d nhstepm=%d \n",agedeb,stepm,hstepm,nhstepm);*/\r
@@ -2852,24 +2892,35 @@ fclose(fichtm);
        for (h=0; h<=nhstepm; h++){\r
        \r
         if (h*hstepm/YEARM*stepm==cpt)\r
-           fprintf(ficresf,"\n%d %.f %.f %.f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm, h*hstepm/YEARM*stepm);\r
+           fprintf(ficresf,"\n%.f %.f %.f",agedeb, agedeb+ h*hstepm/YEARM*stepm, h*hstepm/YEARM*stepm);\r
         \r
          \r
-         for(j=1; j<=nlstate+ndeath;j++) {\r
-           kk1=0.;\r
-           for(i=1; i<=nlstate;i++) {        \r
-             if (mobilav==1)\r
-             kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb][i][cptcod];\r
-             else kk1=kk1+p3mat[i][j][h]*probs[(int)agedeb][i][cptcod];\r
-           }    \r
-         if (h*hstepm/YEARM*stepm==cpt) fprintf(ficresf," %.5f ", kk1);\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][i][cptcod];\r
+            else kk1=kk1+p3mat[i][j][h]*probs[(int)agedeb][i][cptcod];\r
+            if (popforecast==1) kk2=kk1*popeffectif[(int)agedeb];\r
+           }\r
+          if (h*hstepm/YEARM*stepm==cpt) {\r
+            fprintf(ficresf," %.3f", kk1);\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
+  if (popforecast==1) {\r
+    free_ivector(popage,0,AGESUP);\r
+    free_vector(popeffectif,0,AGESUP);\r
+    free_vector(popcount,0,AGESUP);\r
+  }\r
   free_imatrix(s,1,maxwav+1,1,n);\r
   free_vector(weight,1,n);\r
   fclose(ficresf);\r