In freqsummary, the observed prevalence is now computed from all the

waves and the bug in the forecasting routine is corrected. We can
forecast prevalence if ncov=0.

diff --git a/src/imach.c b/src/imach.c
index f2225f763479147475a0876ad3158b6da987428f..b3357c4094ba54e2e677516e8834196c9d69f705 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -83,7 +83,7 @@
 #define ODIRSEPARATOR '\\'
 #endif
 
-char version[80]="Imach version 0.93, February 2003, INED-EUROREVES ";
+char version[80]="Imach version 0.94, February 2003, INED-EUROREVES ";
 int erreur; /* Error number */
 int nvar;
 int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;
@@ -1429,12 +1429,13 @@ void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev
   int i, m, jk, k1,i1, j1, bool, z1,z2,j;
   int first;
   double ***freq; /* Frequencies */
-  double *pp;
-  double pos, k2, dateintsum=0,k2cpt=0;
+  double *pp, **prop;
+  double pos,posprop, k2, dateintsum=0,k2cpt=0;
   FILE *ficresp;
   char fileresp[FILENAMELENGTH];
   
   pp=vector(1,nlstate);
+  prop=matrix(1,nlstate,agemin,agemax+3);
   probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
   strcpy(fileresp,"p");
   strcat(fileresp,fileres);
@@ -1460,6 +1461,10 @@ void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev
        for (jk=-1; jk<=nlstate+ndeath; jk++)  
          for(m=agemin; m <= agemax+3; m++)
            freq[i][jk][m]=0;
+
+    for (i=1; i<=nlstate; i++)  
+      for(m=agemin; m <= agemax+3; m++)
+       prop[i][m]=0;
       
       dateintsum=0;
       k2cpt=0;
@@ -1476,6 +1481,7 @@ void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev
            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;
+             if (s[m][i]>0 && s[m][i]<=nlstate) prop[s[m][i]][(int)agev[m][i]] += weight[i];
              if (m<lastpass) {
                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];
@@ -1533,10 +1539,11 @@ void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev
        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++)
+       }       
+       for(jk=1,pos=0,posprop=0; jk <=nlstate ; jk++){
          pos += pp[jk];
+         posprop += prop[jk][i];
+       }
        for(jk=1; jk <=nlstate ; jk++){
          if(pos>=1.e-5){
            if(first==1)
@@ -1549,12 +1556,12 @@ void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev
          }
          if( i <= (int) agemax){
            if(pos>=1.e-5){
-             fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos);
+             fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop);
              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);
+             fprintf(ficresp," %d NaNq %.0f %.0f",i,prop[jk][i],posprop);
          }
        }
        
@@ -1578,7 +1585,7 @@ void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev
   fclose(ficresp);
   free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
   free_vector(pp,1,nlstate);
-  
+  free_matrix(prop,1,nlstate,(int) agemin,(int) agemax+3);
   /* End of Freq */
 }
 
@@ -1592,12 +1599,12 @@ void prevalence(int agemin, float agemax, int **s, double **agev, int nlstate, i
  
   int i, m, jk, k1, i1, j1, bool, z1,z2,j;
   double ***freq; /* Frequencies */
-  double *pp;
-  double pos; 
+  double *pp, **prop;
+  double pos,posprop; 
   double  y2; /* in fractional years */
 
   pp=vector(1,nlstate);
-  
+  prop=matrix(1,nlstate,agemin,agemax+3); 
   freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);
   j1=0;
   
@@ -1608,10 +1615,9 @@ void prevalence(int agemin, float agemax, int **s, double **agev, int nlstate, i
     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<=nlstate; i++)  
+       for(m=agemin; m <= agemax+3; m++)
+         prop[i][m]=0;
      
       for (i=1; i<=imx; i++) { /* Each individual */
        bool=1;
@@ -1626,35 +1632,24 @@ void prevalence(int agemin, float agemax, int **s, double **agev, int nlstate, i
            if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */
              if(agev[m][i]==0) agev[m][i]=agemax+1;
              if(agev[m][i]==1) agev[m][i]=agemax+2;
-             if (m<lastpass) {
-               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 (s[m][i]>0 && s[m][i]<=nlstate) {
+               prop[s[m][i]][(int)agev[m][i]] += weight[i];
+               prop[s[m][i]][(int)(agemax+3)] += weight[i];
              }
            }
          } /* end selection of waves */
        }
       }
       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,posprop=0; jk <=nlstate ; jk++) {
+         posprop += prop[jk][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;
+           if(posprop>=1.e-5){
+            probs[i][jk][j1]= prop[jk][i]/posprop;
            }
          }
        }/* end jk */
@@ -1665,7 +1660,7 @@ void prevalence(int agemin, float agemax, int **s, double **agev, int nlstate, i
   
   free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
   free_vector(pp,1,nlstate);
-  
+  free_matrix(prop,1,nlstate,(int) agemin,(int) agemax+3);
 }  /* End of Freq */
 
 /************* Waves Concatenation ***************/
@@ -1744,6 +1739,7 @@ void  concatwav(int wav[], int **dh, int **bh,  int **mw, int **s, double *agedc
          if (j >= jmax) jmax=j;
          else if (j <= jmin)jmin=j;
          /*        if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
+         /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/
          sum=sum+j;
        }
        jk= j/stepm;
@@ -3048,7 +3044,7 @@ prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, doubl
      dateprev1 dateprev2 range of dates during which prevalence is computed
      anproj2 year of en of projection (same day and month as proj1).
   */
-  int yearp, stepsize, hstepm, nhstepm, j, k, c, cptcod, i, h;
+  int yearp, stepsize, hstepm, nhstepm, j, k, c, cptcod, i, h, i1;
   int *popage;
   double agec; /* generic age */
   double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
@@ -3093,11 +3089,15 @@ prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, doubl
   jprojmean=yp;
   if(jprojmean==0) jprojmean=1;
   if(mprojmean==0) jprojmean=1;
+
+  i1=cptcoveff;
+  if (cptcovn < 1){i1=1;}
   
   fprintf(ficresf,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2); 
   
   fprintf(ficresf,"#****** Routine prevforecast **\n");
-  for(cptcov=1, k=0;cptcov<=cptcoveff;cptcov++){
+
+  for(cptcov=1, k=0;cptcov<=i1;cptcov++){
     for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
       k=k+1;
       fprintf(ficresf,"\n#******");