Preforecast routine cleaned. Change of notation, comments added.

Version 0.92

diff --git a/src/imach.c b/src/imach.c
index 0ec6e4e38b46fd99b5233d2efeeee9746e3941a0..5b7b7f3030a715e41ca9f8bd913658b2bfe09905 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -83,7 +83,7 @@
 #define ODIRSEPARATOR '\\'
 #endif
 
-char version[80]="Imach version 0.91, November 2002, INED-EUROREVES ";
+char version[80]="Imach version 0.92, February 2003, INED-EUROREVES ";
 int erreur; /* Error number */
 int nvar;
 int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;
@@ -1556,7 +1556,7 @@ void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev
 }
 
 /************ Prevalence ********************/
-void prevalence(int agemin, float agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, double calagedatem)
+void prevalence(int agemin, float agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass)
 {  
   /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people
      in each health status at the date of interview (if between dateprev1 and dateprev2).
@@ -1600,10 +1600,7 @@ void prevalence(int agemin, float agemax, int **s, double **agev, int nlstate, i
              if(agev[m][i]==0) agev[m][i]=agemax+1;
              if(agev[m][i]==1) agev[m][i]=agemax+2;
              if (m<lastpass) {
-               if (calagedatem>0) /* We compute prevalence at exact age, agev in fractional years */
-                 freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedatem %12)/12.)] += weight[i];
-               else
-                 freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
+               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]; 
              }
            }
@@ -2215,13 +2212,12 @@ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double
            vareij[i][j][(int)age] += doldm[i][j]*hf*hf;
       }
     }
-
+  
     /* pptj */
     matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov);
     matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp);
- 
-   for(j=nlstate+1;j<=nlstate+ndeath;j++)
-     for(i=nlstate+1;i<=nlstate+ndeath;i++){
+    for(j=nlstate+1;j<=nlstate+ndeath;j++)
+      for(i=nlstate+1;i<=nlstate+ndeath;i++)
        varppt[j][i]=doldmp[j][i];
     /* end ppptj */
     /*  x centered again */
@@ -2237,7 +2233,7 @@ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double
          prlim[i][i]=mobaverage[(int)age][i][ij];
       }
     }
-    
+             
     /* This for computing probability of death (h=1 means
        computed over hstepm (estepm) matrices product = hstepm*stepm months) 
        as a weighted average of prlim.
@@ -2298,7 +2294,7 @@ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double
   fclose(ficresprobmorprev);
   fclose(ficgp);
   fclose(fichtm);
-}
+}  
 
 /************ Variance of prevlim ******************/
 void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij)
@@ -3018,29 +3014,23 @@ int movingaverage(double ***probs, double bage,double fage, double ***mobaverage
 
 
 /************** Forecasting ******************/
-prevforecast(char fileres[], double anproj1,double mproj1,double jproj1,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anproj2,double p[], int i2){
+prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){
   /* proj1, year, month, day of starting projection 
      agemin, agemax range of age
      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;
   int *popage;
-  double calagedatem, agelim, kk1, kk2, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
+  double agec; /* generic age */
+  double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
   double *popeffectif,*popcount;
   double ***p3mat;
   double ***mobaverage;
   char fileresf[FILENAMELENGTH];
 
   agelim=AGESUP;
-  calagedatem=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM; /* offset
-      from the mean date of interviews between dateprev1 and dateprev2
-      (in fractional years converted to fractional months) */
-  /* Computing age-specific starting prevalence between exact 
-     dateprev1 and dateprev2 (in float years) and
- shifting ages from calagedatem.
-  */
-  prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedatem);
+  prevalence(ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
  
   strcpy(fileresf,"f"); 
   strcat(fileresf,fileres);
@@ -3064,8 +3054,6 @@ prevforecast(char fileres[], double anproj1,double mproj1,double jproj1,double a
   stepsize=(int) (stepm+YEARM-1)/YEARM;
   if (stepm<=12) stepsize=1;
   
-  agelim=AGESUP;
-  
   hstepm=1;
   hstepm=hstepm/stepm; 
   yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp  and
@@ -3078,50 +3066,54 @@ prevforecast(char fileres[], double anproj1,double mproj1,double jproj1,double a
   if(jprojmean==0) jprojmean=1;
   if(mprojmean==0) jprojmean=1;
   
-  fprintf(ficresf,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f ",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2); 
+  fprintf(ficresf,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2); 
   
-  for(cptcov=1, k=0;cptcov<=i2;cptcov++){
+  fprintf(ficresf,"#****** Routine prevforecast **\n");
+  for(cptcov=1, k=0;cptcov<=cptcoveff;cptcov++){
     for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
       k=k+1;
       fprintf(ficresf,"\n#******");
       for(j=1;j<=cptcoveff;j++) {
-       fprintf(ficresf," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+       fprintf(ficresf," V%d=%d, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
       }
       fprintf(ficresf,"******\n");
-      fprintf(ficresf,"# StartingAge FinalAge");
-      for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);
+      fprintf(ficresf,"# Covariate valuofcovar yearproj age");
+      for(j=1; j<=nlstate+ndeath;j++){ 
+       for(i=1; i<=nlstate;i++)              
+          fprintf(ficresf," p%d%d",i,j);
+       fprintf(ficresf," p.%d",j);
+      }
       for (yearp=0; yearp<=(anproj2-anproj1);yearp++) { 
        fprintf(ficresf,"\n");
-       fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);   
+       fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp);   
 
-       for (agedeb=(fage-((int)calagedatem %12/12.)); agedeb>=(ageminpar-((int)calagedatem %12)/12.); agedeb--){ 
-         nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); 
+       for (agec=fage; agec>=ageminpar; agec--){ 
+         nhstepm=(int) rint((agelim-agec)*YEARM/stepm); 
          nhstepm = nhstepm/hstepm; 
-         
          p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
          oldm=oldms;savm=savms;
-         hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  
+         hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);  
        
          for (h=0; h<=nhstepm; h++){
-           if (h==(int) (calagedatem+YEARM*yearp)) {
+           if (h==(int) (YEARM*yearp)) {
               fprintf(ficresf,"\n");
               for(j=1;j<=cptcoveff;j++) 
                 fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
-             fprintf(ficresf,"%.f %.f ",anproj1+yearp,agedeb+h*hstepm/YEARM*stepm);
+             fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm);
            } 
            for(j=1; j<=nlstate+ndeath;j++) {
-             kk1=0.;kk2=0;
+             ppij=0.;
              for(i=1; i<=nlstate;i++) {              
                if (mobilav==1) 
-                 kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod];
+                 ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec+1][i][cptcod];
                else {
-                 kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod];
+                 ppij=ppij+p3mat[i][j][h]*probs[(int)(agec+1)][i][cptcod];
                }
-               
+               if (h==(int)(YEARM*yearp))
+                 fprintf(ficresf," %.3f", p3mat[i][j][h]);
              }
-             if (h==(int)(calagedatem+12*yearp)){
-               fprintf(ficresf," %.3f", kk1);
-                       
+             if (h==(int)(YEARM*yearp)){
+               fprintf(ficresf," %.3f", ppij);
              }
            }
          }
@@ -3135,7 +3127,8 @@ prevforecast(char fileres[], double anproj1,double mproj1,double jproj1,double a
 
   fclose(ficresf);
 }
-/************** Forecasting ******************/
+
+/************** Forecasting *****not tested NB*************/
 populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){
   
   int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;
@@ -3151,7 +3144,7 @@ populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double
   agelim=AGESUP;
   calagedatem=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM;
   
-  prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedatem);
+  prevalence(ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
   
   
   strcpy(filerespop,"pop"); 
@@ -3324,7 +3317,7 @@ int main(int argc, char *argv[])
   int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */
   int mobilav=0,popforecast=0;
   int hstepm, nhstepm;
-  double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1, calagedatem;
+  double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1;
 
   double bage, fage, age, agelim, agebase;
   double ftolpl=FTOL;
@@ -3738,7 +3731,7 @@ int main(int argc, char *argv[])
                agev[m][i]=-1;
              }
            }
-       }<
+       }
        else if(s[m][i] !=9){ /* Standard case, age in fractional
                                 years but with the precision of a
                                 month */
@@ -3972,8 +3965,8 @@ int main(int argc, char *argv[])
   ungetc(c,ficpar);
  
 
-  dateprev1=anprev1+mprev1/12.+jprev1/365.;
-  dateprev2=anprev2+mprev2/12.+jprev2/365.;
+  dateprev1=anprev1+(mprev1-1)/12.+(jprev1-1)/365.;
+  dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/365.;
 
   fscanf(ficpar,"pop_based=%d\n",&popbased);
   fprintf(ficparo,"pop_based=%d\n",popbased);   
@@ -3988,10 +3981,10 @@ int main(int argc, char *argv[])
   ungetc(c,ficpar);
 
   fscanf(ficpar,"prevforecast=%d starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mobil_average=%d\n",&prevfcast,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilavproj);
-  fprintf(ficparo,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",&prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
+  fprintf(ficparo,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
   printf("prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobilaverage=%d\n",&prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
-  fprintf(ficlog,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobilaverage=%d\n",&prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
-  fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobilaverage=%d\n",&prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
+  fprintf(ficlog,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobilaverage=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
+  fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobilaverage=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
   /* day and month of proj2 are not used but only year anproj2.*/
 
   while((c=getc(ficpar))=='#' && c!= EOF){
@@ -4129,6 +4122,7 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
 
   /* hstepm=1;   aff par mois*/
 
+  fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");
   for(cptcov=1,k=0;cptcov<=i1;cptcov++){
     for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
       k=k+1;
@@ -4146,13 +4140,13 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
        p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
        oldm=oldms;savm=savms;
        hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  
-       fprintf(ficrespij,"# Age");
+       fprintf(ficrespij,"# Cov Agex agex+h hpijx with i,j=");
        for(i=1; i<=nlstate;i++)
          for(j=1; j<=nlstate+ndeath;j++)
            fprintf(ficrespij," %1d-%1d",i,j);
        fprintf(ficrespij,"\n");
        for (h=0; h<=nhstepm; h++){
-         fprintf(ficrespij,"%d %f %f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );
+         fprintf(ficrespij,"%d %3.f %3.f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );
          for(i=1; i<=nlstate;i++)
            for(j=1; j<=nlstate+ndeath;j++)
              fprintf(ficrespij," %.5f", p3mat[i][j][h]);
@@ -4173,7 +4167,7 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
   /*if((stepm == 1) && (strcmp(model,".")==0)){*/
   if(prevfcast==1){
     if(stepm ==1){
-      prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilavproj, agedeb, fage, popforecast, popfile, anproj2,p, i1);
+      prevforecast(fileres, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
       if (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);
     } 
     else{
@@ -4214,9 +4208,7 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
   printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
   fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
 
-  calagedatem=-1;
-
-  prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedatem);
+  prevalence(ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
 
   if (mobilav!=0) {
     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);