* imach.c (Repository): Replace "freqsummary" at a correct

place. It differs from routine "prevalence" which may be called
many times. Probs is memory consuming and must be used with
parcimony.
Version 0.95a2 (should output exactly the same maximization than 0.8a2)

diff --git a/src/imach.c b/src/imach.c
index 97f90e7f06a3ada13ff25cee58428d5856bbb204..7a88e88cce60fdc5d4719b6a5eecc23b91083ad0 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -1,6 +1,9 @@
 /* $Id$
   $State$
   $Log$
+  Revision 1.83  2003/06/10 13:39:11  lievre
+  *** empty log message ***
+
   Revision 1.82  2003/06/05 15:57:20  brouard
   Add log in  imach.c and  fullversion number is now printed.
 
@@ -65,6 +68,7 @@
   read parameterfile
   read datafile
   concatwav
+  freqsummary
   if (mle >= 1)
     mlikeli
   print results files
@@ -129,7 +133,7 @@
 /* $Id$ */
 /* $State$ */
 
-char version[]="Imach version 0.95a1, June 2003, INED-EUROREVES ";
+char version[]="Imach version 0.95a2, June 2003, INED-EUROREVES ";
 char fullversion[]="$Revision$ $Date$"; 
 int erreur; /* Error number */
 int nvar;
@@ -1174,7 +1178,7 @@ double func( double *x)
        ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
       } /* end of wave */
     } /* end of individual */
-  }else{  /* ml=4 no inter-extrapolation */
+  }else if (mle==4){  /* ml=4 no inter-extrapolation */
     for (i=1,ipmx=0, sw=0.; i<=imx; i++){
       for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
       for(mi=1; mi<= wav[i]-1; mi++){
@@ -1196,16 +1200,55 @@ double func( double *x)
          oldm=newm;
        } /* end mult */
       
+       s1=s[mw[mi][i]][i];
+       s2=s[mw[mi+1][i]][i];
+       if( s2 > nlstate){ 
+         lli=log(out[s1][s2] - savm[s1][s2]);
+       }else{
+         lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */
+       }
+       ipmx +=1;
+       sw += weight[i];
+       ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
+       /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]);*/
+      } /* end of wave */
+    } /* end of individual */
+  }else{  /* ml=5 no inter-extrapolation no jackson =0.8a */
+    for (i=1,ipmx=0, sw=0.; i<=imx; i++){
+      for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
+      for(mi=1; mi<= wav[i]-1; mi++){
+       for (ii=1;ii<=nlstate+ndeath;ii++)
+         for (j=1;j<=nlstate+ndeath;j++){
+           oldm[ii][j]=(ii==j ? 1.0 : 0.0);
+           savm[ii][j]=(ii==j ? 1.0 : 0.0);
+         }
+       for(d=0; d<dh[mi][i]; d++){
+         newm=savm;
+         cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
+         for (kk=1; kk<=cptcovage;kk++) {
+           cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
+         }
+       
+         out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
+                      1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
+         savm=oldm;
+         oldm=newm;
+       } /* end mult */
+      
+       s1=s[mw[mi][i]][i];
+       s2=s[mw[mi+1][i]][i];
        lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */
        ipmx +=1;
        sw += weight[i];
        ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
+       /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]);*/
       } /* end of wave */
     } /* end of individual */
   } /* End of if */
   for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
   /* printf("l1=%f l2=%f ",ll[1],ll[2]); */
   l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
+  /*exit(0); */
   return -l;
 }
 
@@ -1501,7 +1544,7 @@ void lubksb(double **a, int n, int *indx, double b[])
 } 
 
 /************ Frequencies ********************/
-void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2,double jprev1, double mprev1,double anprev1,double jprev2, double mprev2,double anprev2)
+void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint)
 {  /* Some frequencies */
   
   int i, m, jk, k1,i1, j1, bool, z1,z2,j;
@@ -1555,7 +1598,7 @@ void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **ag
        if (bool==1){
          for(m=firstpass; m<=lastpass; m++){
            k2=anint[m][i]+(mint[m][i]/12.);
-           if ((k2>=dateprev1) && (k2<=dateprev2)) {
+           /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/
              if(agev[m][i]==0) agev[m][i]=iagemax+1;
              if(agev[m][i]==1) agev[m][i]=iagemax+2;
              if (s[m][i]>0 && s[m][i]<=nlstate) prop[s[m][i]][(int)agev[m][i]] += weight[i];
@@ -1568,12 +1611,12 @@ void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **ag
                dateintsum=dateintsum+k2;
                k2cpt++;
              }
-           }
+             /*}*/
          }
        }
       }
        
-      fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
+      /*      fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
 
       if  (cptcovn>0) {
        fprintf(ficresp, "\n#********** Variable "); 
@@ -1634,7 +1677,7 @@ void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **ag
          if( i <= iagemax){
            if(pos>=1.e-5){
              fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop);
-             probs[i][jk][j1]= pp[jk]/pos;
+             /*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
@@ -1667,7 +1710,7 @@ void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **ag
 }
 
 /************ Prevalence ********************/
-void prevalence(double agemin, double 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)
+void prevalence(double ***probs, double agemin, double 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).
@@ -2401,7 +2444,7 @@ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double
   fclose(ficresprobmorprev);
   fclose(ficgp);
   fclose(fichtm);
-}  
+}  /* end varevsij */
 
 /************ 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)
@@ -3137,7 +3180,7 @@ prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, doubl
   char fileresf[FILENAMELENGTH];
 
   agelim=AGESUP;
-  prevalence(ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
+  prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
  
   strcpy(fileresf,"f"); 
   strcat(fileresf,fileres);
@@ -3260,7 +3303,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, firstpass, lastpass);
+  prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
   
   
   strcpy(filerespop,"pop"); 
@@ -3402,7 +3445,7 @@ populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double
   free_ma3x(tabpop,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
   free_ma3x(tabpopprev,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
   fclose(ficrespop);
-}
+} /* End of popforecast */
 
 /***********************************************/
 /**************** Main Program *****************/
@@ -3843,7 +3886,7 @@ int main(int argc, char *argv[])
       if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){
        printf("Error! Month of death of individual %d on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]); 
        fprintf(ficlog,"Error! Month of death of individual %d on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]); 
-       s[m][i]=-1;
+       s[m][i]=-1; /* We prefer to skip it (and to skip it in version 0.8a1 too */
       }
     }
   }
@@ -3965,6 +4008,7 @@ int main(int argc, char *argv[])
     
   /* Calculates basic frequencies. Computes observed prevalence at single age
      and prints on file fileres'p'. */
+  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint);
 
     pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
     oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
@@ -4008,7 +4052,7 @@ int main(int argc, char *argv[])
        }
     }
   }
-  if(mle==1){
+  if(mle!=0){
     /* Computing hessian and covariance matrix */
     ftolhess=ftol; /* Usually correct */
     hesscov(matcov, p, npar, delti, ftolhess, func);
@@ -4139,8 +4183,8 @@ int main(int argc, char *argv[])
   fprintf(ficparo,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);
   fprintf(ficres,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);
 
-  probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
-  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
+  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint);
+  /*,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
 
   /*------------ gnuplot -------------*/
   strcpy(optionfilegnuplot,optionfilefiname);
@@ -4304,6 +4348,7 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
 
   fclose(ficrespij);
 
+  probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
 
   /*---------- Forecasting ------------------*/
   /*if((stepm == 1) && (strcmp(model,".")==0)){*/
@@ -4351,7 +4396,7 @@ 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);
 
   /* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */
-  prevalence(agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
+  prevalence(probs, agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
   /*  printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d,  mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\
 ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass);
   */