--- imach096d/src/imach.c	2003/05/02 18:51:41	1.74
+++ imach096d/src/imach.c	2003/05/16 10:44:42	1.76
@@ -1,4 +1,4 @@
-/* $Id: imach.c,v 1.74 2003/05/02 18:51:41 brouard Exp $
+/* $Id: imach.c,v 1.76 2003/05/16 10:44:42 brouard Exp $
    Interpolated Markov Chain
 
   Short summary of the programme:
@@ -119,7 +119,7 @@
 #define ODIRSEPARATOR '\\'
 #endif
 
-char version[80]="Imach version 0.95, February 2003, INED-EUROREVES ";
+char version[80]="Imach version 0.95a, May 2003, INED-EUROREVES ";
 int erreur; /* Error number */
 int nvar;
 int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;
@@ -141,7 +141,7 @@ double jmean; /* Mean space between 2 wa
 double **oldm, **newm, **savm; /* Working pointers to matrices */
 double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
 FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
-FILE *ficlog;
+FILE *ficlog, *ficrespow;
 FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor;
 FILE *ficresprobmorprev;
 FILE *fichtm; /* Html File */
@@ -327,6 +327,21 @@ void free_vector(double*v, int nl, int n
 }
 
 /************************ivector *******************************/
+char *cvector(long nl,long nh)
+{
+  char *v;
+  v=(char *) malloc((size_t)((nh-nl+1+NR_END)*sizeof(char)));
+  if (!v) nrerror("allocation failure in cvector");
+  return v-nl+NR_END;
+}
+
+/******************free ivector **************************/
+void free_cvector(char *v, long nl, long nh)
+{
+  free((FREE_ARG)(v+nl-NR_END));
+}
+
+/************************ivector *******************************/
 int *ivector(long nl,long nh)
 {
   int *v;
@@ -647,11 +662,15 @@ void powell(double p[], double **xi, int
     del=0.0; 
     printf("\nPowell iter=%d -2*LL=%.12f",*iter,*fret);
     fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f",*iter,*fret);
-    for (i=1;i<=n;i++) 
+    fprintf(ficrespow,"%d %.12f",*iter,*fret);
+    for (i=1;i<=n;i++) {
       printf(" %d %.12f",i, p[i]);
-    fprintf(ficlog," %d %.12f",i, p[i]);
+      fprintf(ficlog," %d %.12lf",i, p[i]);
+      fprintf(ficrespow," %.12lf", p[i]);
+    }
     printf("\n");
     fprintf(ficlog,"\n");
+    fprintf(ficrespow,"\n");
     for (i=1;i<=n;i++) { 
       for (j=1;j<=n;j++) xit[j]=xi[j][i]; 
       fptt=(*fret); 
@@ -1187,14 +1206,27 @@ void mlikeli(FILE *ficres,double p[], in
   int i,j, iter;
   double **xi;
   double fret;
+  char filerespow[FILENAMELENGTH];
   xi=matrix(1,npar,1,npar);
   for (i=1;i<=npar;i++)
     for (j=1;j<=npar;j++)
       xi[i][j]=(i==j ? 1.0 : 0.0);
   printf("Powell\n");  fprintf(ficlog,"Powell\n");
+  strcpy(filerespow,"pow"); 
+  strcat(filerespow,fileres);
+  if((ficrespow=fopen(filerespow,"w"))==NULL) {
+    printf("Problem with resultfile: %s\n", filerespow);
+    fprintf(ficlog,"Problem with resultfile: %s\n", filerespow);
+  }
+  fprintf(ficrespow,"# Powell\n# iter -2*LL");
+  for (i=1;i<=nlstate;i++)
+    for(j=1;j<=nlstate+ndeath;j++)
+      if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);
+  fprintf(ficrespow,"\n");
   powell(p,xi,npar,ftol,&iter,&fret,func);
 
-   printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));
+  fclose(ficrespow);
+  printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));
   fprintf(ficlog,"\n#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
   fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
 
@@ -2792,7 +2824,7 @@ fprintf(fichtm," \n<ul><li><b>Graphs</b>
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
      }
      /* Pij */
-     fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before: pe%s%d1.png<br>
+     fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: pe%s%d1.png<br>
 <img src=\"pe%s%d1.png\">",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);     
      /* Quasi-incidences */
      fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: pe%s%d2.png<br>
@@ -2823,12 +2855,12 @@ health expectancies in states (1) and (2
  - Health expectancies with their variances (no covariance): <a href=\"t%s\">t%s</a> <br>\n
  - Standard deviation of stable prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);
 
- if(popforecast==1) fprintf(fichtm,"\n
- - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n
- - Population forecasting (if popforecast=1): <a href=\"pop%s\">pop%s</a> <br>\n
-	<br>",fileres,fileres,fileres,fileres);
- else 
-   fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)<br><br></li>\n",popforecast, stepm, model);
+/*  if(popforecast==1) fprintf(fichtm,"\n */
+/*  - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n */
+/*  - Population forecasting (if popforecast=1): <a href=\"pop%s\">pop%s</a> <br>\n */
+/* 	<br>",fileres,fileres,fileres,fileres); */
+/*  else  */
+/*    fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)<br><br></li>\n",popforecast, stepm, model); */
 fprintf(fichtm," <ul><li><b>Graphs</b></li><p>");
 
  m=cptcoveff;
@@ -2845,7 +2877,7 @@ fprintf(fichtm," <ul><li><b>Graphs</b></
        fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
      }
      for(cpt=1; cpt<=nlstate;cpt++) {
-       fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident
+       fprintf(fichtm,"<br>- Observed and period prevalence (with confident
 interval) in state (%d): v%s%d%d.png <br>
 <img src=\"v%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);  
      }
@@ -2946,14 +2978,14 @@ m=pow(2,cptcoveff);
     }
   }
   
-  /* CV preval stat */
+  /* CV preval stable (period) */
   for (k1=1; k1<= m ; k1 ++) { 
-    for (cpt=1; cpt<nlstate ; cpt ++) {
+    for (cpt=1; cpt<=nlstate ; cpt ++) {
       k=3;
       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);
       
@@ -3138,6 +3170,7 @@ prevforecast(char fileres[], double anpr
   
   fprintf(ficresf,"#****** Routine prevforecast **\n");
 
+/* 	      if (h==(int)(YEARM*yearp)){ */
   for(cptcov=1, k=0;cptcov<=i1;cptcov++){
     for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
       k=k+1;
@@ -3164,7 +3197,7 @@ prevforecast(char fileres[], double anpr
 	  hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);  
 	
 	  for (h=0; h<=nhstepm; h++){
-	    if (h==(int) (YEARM*yearp)) {
+	    if (h*hstepm/YEARM*stepm ==yearp) {
               fprintf(ficresf,"\n");
               for(j=1;j<=cptcoveff;j++) 
                 fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
@@ -3178,19 +3211,20 @@ prevforecast(char fileres[], double anpr
 		else {
 		  ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod];
 		}
-		if (h==(int)(YEARM*yearp))
+		if (h*hstepm/YEARM*stepm== yearp) {
 		  fprintf(ficresf," %.3f", p3mat[i][j][h]);
-	      }
-	      if (h==(int)(YEARM*yearp)){
+		}
+	      } /* end i */
+	      if (h*hstepm/YEARM*stepm==yearp) {
 		fprintf(ficresf," %.3f", ppij);
 	      }
-	    }
-	  }
+	    }/* end j */
+	  } /* end h */
 	  free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
-	}
-      }
-    }
-  }
+	} /* end agec */
+      } /* end yearp */
+    } /* end cptcod */
+  } /* end  cptcov */
        
   if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
 
@@ -3784,11 +3818,21 @@ int main(int argc, char *argv[])
 
   for (i=1; i<=imx; i++) {
     for(m=2; (m<= maxwav); m++) {
-      if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){
+      if (((int)mint[m][i]== 99) && (s[m][i] <= nlstate)){
 	anint[m][i]=9999;
 	s[m][i]=-1;
       }
-      if(moisdc[i]==99 && andc[i]==9999 & s[m][i]>nlstate) s[m][i]=-1;
+      if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){
+	printf("Error! Date of death (month %2d and year %4d) of individual %d on line %d was unknown %d, set an arbitrary year of death\n",(int)moisdc[i],(int)andc[i],num[i],i);
+	fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %d on line %d was unknown %d, set an arbitrary year of death\n",(int)moisdc[i],(int)andc[i],num[i],i);
+	s[m][i]=-1;
+      }
+      if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){
+	printf("Warning! Month of death of individual %d on line %d was unknown %2d, we set it to June\n",num[i],i,(int)moisdc[i]); 
+	fprintf(ficlog,"Warning! Month of death of individual %d on line %d was unknown %f, we set it to June\n",num[i],i,moisdc[i]); 
+	moisdc[i]=6;
+	s[m][i]=-1;
+      }
     }
   }
 
@@ -3798,11 +3842,11 @@ int main(int argc, char *argv[])
       if(s[m][i] >0){
 	if (s[m][i] >= nlstate+1) {
 	  if(agedc[i]>0)
-	    if(moisdc[i]!=99 && andc[i]!=9999)
+	    if((int)moisdc[i]!=99 && (int)andc[i]!=9999)
 	      agev[m][i]=agedc[i];
 	  /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/
 	    else {
-	      if (andc[i]!=9999){
+	      if ((int)andc[i]!=9999){
 		printf("Warning negative age at death: %d line:%d\n",num[i],i);
 		fprintf(ficlog,"Warning negative age at death: %d line:%d\n",num[i],i);
 		agev[m][i]=-1;
@@ -3813,7 +3857,7 @@ int main(int argc, char *argv[])
 				 years but with the precision of a
 				 month */
 	  agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]);
-	  if(mint[m][i]==99 || anint[m][i]==9999)
+	  if((int)mint[m][i]==99 || (int)anint[m][i]==9999)
 	    agev[m][i]=1;
 	  else if(agev[m][i] <agemin){ 
 	    agemin=agev[m][i];
@@ -4111,12 +4155,13 @@ int main(int argc, char *argv[])
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n
 \n
 Total number of observations=%d <br>\n
+Youngest age at first pass %.2f, oldest age %.2f<br>\n
 Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
 <hr  size=\"2\" color=\"#EC5E5E\">
  <ul><li><h4>Parameter files</h4>\n
  - Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n
  - Log file of the run: <a href=\"%s\">%s</a><br>\n
- - Gnuplot file name: <a href=\"%s\">%s</a></ul>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,filelog,filelog,optionfilegnuplot,optionfilegnuplot);
+ - Gnuplot file name: <a href=\"%s\">%s</a></ul>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,agemin,agemax,jmin,jmax,jmean,fileres,fileres,filelog,filelog,optionfilegnuplot,optionfilegnuplot);
    fclose(fichtm);
 
   printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);
@@ -4468,8 +4513,9 @@ ageminpar, agemax, s[lastpass][imx], age
   strcpy(plotcmd,GNUPLOTPROGRAM);
   strcat(plotcmd," ");
   strcat(plotcmd,optionfilegnuplot);
-  printf("Starting: %s\n",plotcmd);fflush(stdout);
+  printf("Starting graphs with: %s",plotcmd);fflush(stdout);
   system(plotcmd);
+  printf(" Wait...");
 
  /*#ifdef windows*/
   while (z[0] != 'q') {