--- imach/src/imach.c	2002/11/19 00:17:38	1.60
+++ imach/src/imach.c	2002/11/19 14:08:13	1.61
@@ -1,4 +1,4 @@
-/* $Id: imach.c,v 1.60 2002/11/19 00:17:38 brouard Exp $
+/* $Id: imach.c,v 1.61 2002/11/19 14:08:13 brouard Exp $
    Interpolated Markov Chain
 
   Short summary of the programme:
@@ -928,56 +928,86 @@ double func( double *x)
   cov[1]=1.;
 
   for(k=1; k<=nlstate; k++) ll[k]=0.;
-  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;
+
+  if(mle==1){
+    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 */
+	
+	} /* end mult */
       
-      /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
-      /* But now since version 0.9 we anticipate for bias and large stepm.
-       * If stepm is larger than one month (smallest stepm) and if the exact delay 
-       * (in months) between two waves is not a multiple of stepm, we rounded to 
-       * the nearest (and in case of equal distance, to the lowest) interval but now
-       * we keep into memory the bias bh[mi][i] and also the previous matrix product
-       * (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the
-       * probability in order to take into account the bias as a fraction of the way
-       * from savm to out if bh is neagtive or even beyond if bh is positive. bh varies
-       * -stepm/2 to stepm/2 .
-       * For stepm=1 the results are the same as for previous versions of Imach.
-       * For stepm > 1 the results are less biased than in previous versions. 
-       */
-      s1=s[mw[mi][i]][i];
-      s2=s[mw[mi+1][i]][i];
-      bbh=(double)bh[mi][i]/(double)stepm;
-      lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.-bbh)*out[s1][s2]));
-      /*lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.-bbh)*out[s1][s2]));*/
-      /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/
-      /*if(lli ==000.0)*/
-      /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */
-      ipmx +=1;
-      sw += weight[i];
-      ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
-    } /* end of wave */
-  } /* end of individual */
-
+	/*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
+	/* But now since version 0.9 we anticipate for bias and large stepm.
+	 * If stepm is larger than one month (smallest stepm) and if the exact delay 
+	 * (in months) between two waves is not a multiple of stepm, we rounded to 
+	 * the nearest (and in case of equal distance, to the lowest) interval but now
+	 * we keep into memory the bias bh[mi][i] and also the previous matrix product
+	 * (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the
+	 * probability in order to take into account the bias as a fraction of the way
+	 * from savm to out if bh is neagtive or even beyond if bh is positive. bh varies
+	 * -stepm/2 to stepm/2 .
+	 * For stepm=1 the results are the same as for previous versions of Imach.
+	 * For stepm > 1 the results are less biased than in previous versions. 
+	 */
+	s1=s[mw[mi][i]][i];
+	s2=s[mw[mi+1][i]][i];
+	bbh=(double)bh[mi][i]/(double)stepm;
+	lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.-bbh)*out[s1][s2]));
+	/*lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.-bbh)*out[s1][s2]));*/
+	/*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/
+	/*if(lli ==000.0)*/
+	/*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */
+	ipmx +=1;
+	sw += weight[i];
+	ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
+      } /* end of wave */
+    } /* end of individual */
+  }  else{ 
+    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 */
+      
+	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;
+      } /* 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 */
@@ -3122,7 +3152,7 @@ populforecast(char fileres[], double anp
 
 int main(int argc, char *argv[])
 {
-
+  int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav);
   int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod;
   double agedeb, agefin,hf;
   double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;
@@ -3160,7 +3190,6 @@ int main(int argc, char *argv[])
   double *epj, vepp;
   double kk1, kk2;
   double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;
-  /*int *movingaverage; */
 
   char *alph[]={"a","a","b","c","d","e"}, str[4];
 
@@ -3663,7 +3692,7 @@ int main(int argc, char *argv[])
      so we point p on param[1][1] so that p[1] maps on param[1][1][1] */
   p=param[1][1]; /* *(*(*(param +1)+1)+0) */
 
-  if(mle==1){
+  if(mle>=1){ /* Could be 1 or 2 */
     mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
   }