--- imach096d/src/imach.c	2003/05/16 16:49:47	1.77
+++ imach096d/src/imach.c	2003/06/13 21:44:43	1.84
@@ -1,4 +1,21 @@
-/* $Id: imach.c,v 1.77 2003/05/16 16:49:47 brouard 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 ***
+
+  Revision 1.82  2003/06/05 15:57:20  brouard
+  Add log in  imach.c and  fullversion number is now printed.
+
+*/
+/*
    Interpolated Markov Chain
 
   Short summary of the programme:
@@ -58,6 +75,7 @@
   read parameterfile
   read datafile
   concatwav
+  freqsummary
   if (mle >= 1)
     mlikeli
   print results files
@@ -119,7 +137,11 @@
 #define ODIRSEPARATOR '\\'
 #endif
 
-char version[80]="Imach version 0.95a, May 2003, INED-EUROREVES ";
+/* $Id: imach.c,v 1.84 2003/06/13 21:44:43 brouard Exp $ */
+/* $State: Exp $ */
+
+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;
@@ -1163,7 +1185,42 @@ 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++){
+	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];
+	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++){
@@ -1185,16 +1242,20 @@ double func( double *x)
 	  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;
 }
 
@@ -1490,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;
@@ -1544,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];
@@ -1557,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 "); 
@@ -1623,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
@@ -1656,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).
@@ -1799,7 +1860,7 @@ void  concatwav(int wav[], int **dh, int
 	  sum=sum+j;
 	  /*if (j<0) printf("j=%d num=%d \n",j,i);*/
 	  /*	  printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/
-	  if(j<0)printf("Warning! Negative delay (%d) between waves %d and %d of individual at line %d who is aged %.1f with statuses %d %d\n ",j,mw[mi][i],mw[mi+1][i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+	  if(j<0)printf("Error! Negative delay (%d to death) between waves %d and %d of individual %d at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
 	  }
 	}
 	else{
@@ -1810,7 +1871,7 @@ void  concatwav(int wav[], int **dh, int
 	  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]);*/
-	  if(j<0)printf("Warning! Negative delay (%d to death) between waves %d and %d of individual at line %d who is aged %.1f with statuses %d %d\n ",j,mw[mi][i],mw[mi+1][i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+	  if(j<0)printf("Error! Negative delay (%d) between waves %d and %d of individual %d at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
 	  sum=sum+j;
 	}
 	jk= j/stepm;
@@ -2390,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)
@@ -2987,7 +3048,7 @@ m=pow(2,cptcoveff);
       fprintf(ficgp,"\nset out \"p%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,fileres,k1,k+cpt+1,k+1);
       
-      for (i=1; i<= nlstate ; i ++)
+      for (i=1; i< nlstate ; i ++)
 	fprintf(ficgp,"+$%d",k+i+1);
       fprintf(ficgp,")) t\"prev(%d,%d)\" w l",cpt,cpt+1);
       
@@ -3126,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);
@@ -3249,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"); 
@@ -3391,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 *****************/
@@ -3454,7 +3515,7 @@ int main(int argc, char *argv[])
      gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
   getcwd(pathcd, size);
 
-  printf("\n%s",version);
+  printf("\n%s\n%s",version,fullversion);
   if(argc <=1){
     printf("\nEnter the parameter file name: ");
     scanf("%s",pathtot);
@@ -3832,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 */
       }
     }
   }
@@ -3954,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 */
@@ -3997,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);
@@ -4128,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);
@@ -4293,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)){*/
@@ -4340,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);
   */