--- imach/src/imach.c	2003/06/10 13:39:11	1.83
+++ imach/src/imach.c	2003/06/13 21:44:43	1.84
@@ -1,6 +1,13 @@
-/* $Id: imach.c,v 1.83 2003/06/10 13:39:11 lievre Exp $
+/* $Id: imach.c,v 1.84 2003/06/13 21:44:43 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.84  2003/06/13 21:44:43  brouard
+  * 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)
+
   Revision 1.83  2003/06/10 13:39:11  lievre
   *** empty log message ***
 
@@ -68,6 +75,7 @@
   read parameterfile
   read datafile
   concatwav
+  freqsummary
   if (mle >= 1)
     mlikeli
   print results files
@@ -129,11 +137,11 @@
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.83 2003/06/10 13:39:11 lievre Exp $ */
+/* $Id: imach.c,v 1.84 2003/06/13 21:44:43 brouard Exp $ */
 /* $State: Exp $ */
 
-char version[]="Imach version 0.95a1, June 2003, INED-EUROREVES ";
-char fullversion[]="$Revision: 1.83 $ $Date: 2003/06/10 13:39:11 $"; 
+char version[]="Imach version 0.95a2, June 2003, INED-EUROREVES ";
+char fullversion[]="$Revision: 1.84 $ $Date: 2003/06/13 21:44:43 $"; 
 int erreur; /* Error number */
 int nvar;
 int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;
@@ -1177,7 +1185,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++){
@@ -1199,16 +1207,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;
 }
 
@@ -1504,7 +1551,7 @@ void lubksb(double **a, int n, int *indx
 } 
 
 /************ 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;
@@ -1558,7 +1605,7 @@ void  freqsummary(char fileres[], int ia
 	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];
@@ -1571,12 +1618,12 @@ void  freqsummary(char fileres[], int ia
 		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 "); 
@@ -1637,7 +1684,7 @@ void  freqsummary(char fileres[], int ia
 	  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
@@ -1670,7 +1717,7 @@ void  freqsummary(char fileres[], int ia
 }
 
 /************ 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).
@@ -2404,7 +2451,7 @@ void varevsij(char optionfilefiname[], d
   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)
@@ -3140,7 +3187,7 @@ prevforecast(char fileres[], double anpr
   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);
@@ -3263,7 +3310,7 @@ populforecast(char fileres[], double anp
   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"); 
@@ -3405,7 +3452,7 @@ populforecast(char fileres[], double anp
   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 *****************/
@@ -3846,7 +3893,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 */
       }
     }
   }
@@ -3968,6 +4015,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 */
@@ -4011,7 +4059,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);
@@ -4142,8 +4190,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);
@@ -4307,6 +4355,7 @@ Interval (in months) between two waves:
 
   fclose(ficrespij);
 
+  probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
 
   /*---------- Forecasting ------------------*/
   /*if((stepm == 1) && (strcmp(model,".")==0)){*/
@@ -4354,7 +4403,7 @@ Interval (in months) between two waves:
   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);
   */