--- imach/src/imach.c	2006/01/24 19:37:15	1.109
+++ imach/src/imach.c	2006/01/25 00:51:50	1.110
@@ -1,6 +1,9 @@
-/* $Id: imach.c,v 1.109 2006/01/24 19:37:15 brouard Exp $
+/* $Id: imach.c,v 1.110 2006/01/25 00:51:50 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.110  2006/01/25 00:51:50  brouard
+  (Module): Lots of cleaning and bugs added (Gompertz)
+
   Revision 1.109  2006/01/24 19:37:15  brouard
   (Module): Comments (lines starting with a #) are allowed in data.
 
@@ -281,11 +284,11 @@ extern int errno;
 #define ODIRSEPARATOR '/'
 #endif
 
-/* $Id: imach.c,v 1.109 2006/01/24 19:37:15 brouard Exp $ */
+/* $Id: imach.c,v 1.110 2006/01/25 00:51:50 brouard Exp $ */
 /* $State: Exp $ */
 
 char version[]="Imach version 0.98a, January 2006, INED-EUROREVES ";
-char fullversion[]="$Revision: 1.109 $ $Date: 2006/01/24 19:37:15 $"; 
+char fullversion[]="$Revision: 1.110 $ $Date: 2006/01/25 00:51:50 $"; 
 int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
 int nvar;
 int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;
@@ -298,6 +301,7 @@ int popbased=0;
 int *wav; /* Number of waves for this individuual 0 is possible */
 int maxwav; /* Maxim number of waves */
 int jmin, jmax; /* min, max spacing between 2 waves */
+int ijmin, ijmax; /* Individuals having jmin and jmax */ 
 int gipmx, gsw; /* Global variables on the number of contributions 
 		   to the likelihood and the sum of weights (done by funcone)*/
 int mle, weightopt;
@@ -2228,8 +2232,14 @@ void  concatwav(int wav[], int **dh, int
 	      fprintf(ficlog,"   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
 	    }
 	    k=k+1;
-	    if (j >= jmax) jmax=j;
-	    if (j <= jmin) jmin=j;
+	    if (j >= jmax){
+	      jmax=j;
+	      ijmax=i;
+	    }
+	    if (j <= jmin){
+	      jmin=j;
+	      ijmin=i;
+	    }
 	    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);*/
@@ -2240,8 +2250,14 @@ void  concatwav(int wav[], int **dh, int
 /* 	  if (j<0) printf("%d %lf %lf %d %d %d\n", i,agev[mw[mi+1][i]][i], agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); */
 
 	  k=k+1;
-	  if (j >= jmax) jmax=j;
-	  else if (j <= jmin)jmin=j;
+	  if (j >= jmax) {
+	    jmax=j;
+	    ijmax=i;
+	  }
+	  else if (j <= jmin){
+	    jmin=j;
+	    ijmin=i;
+	  }
 	  /*	    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){
@@ -2284,8 +2300,8 @@ void  concatwav(int wav[], int **dh, int
     } /* end wave */
   }
   jmean=sum/k;
-  printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);
-  fprintf(ficlog,"Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);
+  printf("Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, num[ijmin], jmax, num[ijmax], jmean);
+  fprintf(ficlog,"Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, ijmin, jmax, ijmax, jmean);
  }
 
 /*********** Tricode ****************************/
@@ -3996,6 +4012,7 @@ double gompertz(double x[])
 { 
   double A,B,L=0.0,sump=0.,num=0.;
   int i,n=0; /* n is the size of the sample */
+
   for (i=0;i<=imx-1 ; i++) {
     sump=sump+weight[i];
     /*    sump=sump+1;*/
@@ -4008,14 +4025,15 @@ double gompertz(double x[])
 
   for (i=1;i<=imx ; i++)
     {
-      if (cens[i]==1 & wav[i]>1)
+      if (cens[i] == 1 && wav[i]>1)
 	A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)));
       
-      if (cens[i]==0 & wav[i]>1)
+      if (cens[i] == 0 && wav[i]>1)
 	A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)))
 	     +log(x[1]/YEARM)+x[2]*(agedc[i]-agegomp)+log(YEARM);  
       
-      if (wav[i]>1 & agecens[i]>15) {
+      /*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */
+      if (wav[i] > 1 ) { /* ??? */
 	L=L+A*weight[i];
 	/* 	printf("\ni=%d A=%f L=%lf x[1]=%lf x[2]=%lf ageex=%lf agecens=%lf cens=%d agedc=%lf weight=%lf\n",i,A,L,x[1],x[2],ageexmed[i]*12,agecens[i]*12,cens[i],agedc[i]*12,weight[i]);*/
       }
@@ -4084,7 +4102,7 @@ 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=100,cptcode, cptcod;
-  int linei;
+  int linei, month, year,iout;
   int jj, ll, li, lj, lk, imk;
   int numlinepar=0; /* Current linenumber of parameter file */
   int itimes;
@@ -4492,14 +4510,16 @@ int main(int argc, char *argv[])
 
   i=1;
   linei=0;
-  while ((fgets(line, MAXLINE, fic) != NULL) ||((i >= firstobs) && (i <=lastobs)))    {
+  while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs)))    {
     linei=linei+1;
-    printf("IIIII= %d linei=%d\n",i,linei);
     for(j=strlen(line); j>=0;j--){  /* Untabifies line */
 	if(line[j] == '\t')
 	  line[j] = ' ';
       }
-      for(j=strlen(line)-1; (line[j]==' ')||(line[j]==10);j--){;};line[j+1]=0;  /* Trims blanks at end of line */
+      for(j=strlen(line)-1; (line[j]==' ')||(line[j]==10)||(line[j]==13);j--){
+	;
+      };
+      line[j+1]=0;  /* Trims blanks at end of line */
       if(line[0]=='#'){
 	fprintf(ficlog,"Comment line\n%s\n",line);
 	printf("Comment line\n%s\n",line);
@@ -4517,65 +4537,47 @@ int main(int argc, char *argv[])
 	s[j][i]=lval;
 
 	strcpy(line,stra);
-	cutv(stra, strb,line,'/');
-        errno=0;
-        lval=strtol(strb,&endptr,10); 
-	if( strb[0]=='\0' || (*endptr != '\0')){
-	  printf("Error reading data around '%d'.at line number %ld %s for individual %d\nShould be a year of exam at wave %d.  Exiting.\n",lval, i,line,linei,j);
-	  exit(1);
-	}
-	anint[j][i]=(double)(lval); 
-
-	strcpy(line,stra);
 	cutv(stra, strb,line,' ');
-        errno=0;
-        lval=strtol(strb,&endptr,10); 
-	if( strb[0]=='\0' || (*endptr != '\0')){
-	  printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a month of exam at wave %d.  Exiting.\n",lval, i,line, linei,j);
+	if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
+	}
+	else  if(iout=sscanf(strb,".") != 0){
+	  month=99;
+	  year=9999;
+	}else{
+	  printf("Error reading data around '%s'.at line number %ld %s for individual %d\nShould be a year of exam at wave %d.  Exiting.\n",strb, i,line,linei,j);
 	  exit(1);
 	}
-	mint[j][i]=(double)(lval); 
+	anint[j][i]= (double) year; 
+	mint[j][i]= (double)month; 
 	strcpy(line,stra);
-      }
+      } /* ENd Waves */
 	
-      cutv(stra, strb,line,'/'); 
-      errno=0;
-      lval=strtol(strb,&endptr,10); 
-      if( strb[0]=='\0' || (*endptr != '\0')){
-	printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a year of death.  Exiting.\n",lval, i,line,linei);
+      cutv(stra, strb,line,' '); 
+      if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
+      }
+      else  if(iout=sscanf(strb,".") != 0){
+	month=99;
+	year=9999;
+      }else{
+	printf("Error reading data around '%s'.at line number %ld %s for individual %d\nShould be a year of exam at wave %d.  Exiting.\n",strb, i,line,linei,j);
 	exit(1);
       }
-      andc[i]=(double)(lval); 
+      andc[i]=(double) year; 
+      moisdc[i]=(double) month; 
       strcpy(line,stra);
 
       cutv(stra, strb,line,' '); 
-      errno=0;
-      lval=strtol(strb,&endptr,10); 
-      if( strb[0]=='\0' || (*endptr != '\0')){
-	printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a month of death.  Exiting.\n",lval,i,line, linei);
-	exit(1);
-      }
-      moisdc[i]=(double)(lval); 
-
-      strcpy(line,stra);
-      cutv(stra, strb,line,'/'); 
-      errno=0;
-      lval=strtol(strb,&endptr,10); 
-      if( strb[0]=='\0' || (*endptr != '\0')){
-	printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a year of birth.  Exiting.\n",lval, i,line, linei);
-	exit(1);
+      if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
       }
-      annais[i]=(double)(lval);
-
-      strcpy(line,stra);
-      cutv(stra, strb,line,' ');
-      errno=0;
-      lval=strtol(strb,&endptr,10); 
-      if( strb[0]=='\0' || (*endptr != '\0')){
-	printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a month of birth.  Exiting.\n",lval,i,line,linei);
+      else  if(iout=sscanf(strb,".") != 0){
+	month=99;
+	year=9999;
+      }else{
+	printf("Error reading data around '%s'.at line number %ld %s for individual %d\nShould be a year of exam at wave %d.  Exiting.\n",strb, i,line,linei,j);
 	exit(1);
       }
-      moisnais[i]=(double)(lval); 
+      annais[i]=(double)(year);
+      moisnais[i]=(double)(month); 
       strcpy(line,stra);
 
       cutv(stra, strb,line,' '); 
@@ -4596,7 +4598,7 @@ int main(int argc, char *argv[])
 	  printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a covar (meaning 0 for the reference or 1).  Exiting.\n",lval, i,line,linei);
 	  exit(1);
 	}
-	if(lval <0 || lval >1){
+	if(lval <-1 || lval >1){
 	  printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a value of the %d covar (meaning 0 for the reference or 1. IMaCh does not build design variables, do it your self).  Exiting.\n",lval,i,line,linei,j);
 	  exit(1);
 	}
@@ -4611,7 +4613,6 @@ int main(int argc, char *argv[])
       }
       else
 	num[i]=atol(stra);
-      printf ("num [i] %ld %d\n",i, num[i]);fflush(stdout);
       /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){
 	printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/
 
@@ -4940,6 +4941,7 @@ Interval (in months) between two waves:
   p=param[1][1]; /* *(*(*(param +1)+1)+0) */
 
   globpr=0; /* To get the number ipmx of contributions and the sum of weights*/
+
   if (mle==-3){
     ximort=matrix(1,NDIM,1,NDIM);
     cens=ivector(1,n);
@@ -4949,9 +4951,9 @@ Interval (in months) between two waves:
  
     for (i=1; i<=imx; i++){
       dcwave[i]=-1;
-      for (j=1; j<=lastpass; j++)
-	if (s[j][i]>nlstate) {
-	  dcwave[i]=j;
+      for (m=firstpass; m<=lastpass; m++)
+	if (s[m][i]>nlstate) {
+	  dcwave[i]=m;
 	  /*	printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/
 	  break;
 	}
@@ -4960,12 +4962,16 @@ Interval (in months) between two waves:
     for (i=1; i<=imx; i++) {
       if (wav[i]>0){
 	ageexmed[i]=agev[mw[1][i]][i];
-	j=wav[i];agecens[i]=1.; 
-	if (ageexmed[i]>1 & wav[i]>0) agecens[i]=agev[mw[j][i]][i];
-	cens[i]=1;
-	
-	if (ageexmed[i]<1) cens[i]=-1;
-	if (agedc[i]< AGESUP & agedc[i]>1 & dcwave[i]>firstpass & dcwave[i]<=lastpass) cens[i]=0 ;
+	j=wav[i];
+	agecens[i]=1.; 
+
+	if (ageexmed[i]> 1 && wav[i] > 0){
+	  agecens[i]=agev[mw[j][i]][i];
+	  cens[i]= 1;
+	}else if (ageexmed[i]< 1) 
+	  cens[i]= -1;
+	if (agedc[i]< AGESUP && agedc[i]>1 && dcwave[i]>firstpass && dcwave[i]<=lastpass)
+	  cens[i]=0 ;
       }
       else cens[i]=-1;
     }
@@ -4974,29 +4980,29 @@ Interval (in months) between two waves:
       for (j=1;j<=NDIM;j++)
 	ximort[i][j]=(i == j ? 1.0 : 0.0);
     }
-
-    p[1]=0.1; p[2]=0.1;
+    
+    p[1]=0.1; p[NDIM]=0.1;
     /*printf("%lf %lf", p[1], p[2]);*/
     
     
-  printf("Powell\n");  fprintf(ficlog,"Powell\n");
-  strcpy(filerespow,"pow-mort"); 
-  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");
-
+    printf("Powell\n");  fprintf(ficlog,"Powell\n");
+    strcpy(filerespow,"pow-mort"); 
+    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,ximort,NDIM,ftol,&iter,&fret,gompertz);
     fclose(ficrespow);
     
-    hesscov(matcov, p, NDIM,delti, 1e-4, gompertz); 
+    hesscov(matcov, p, NDIM, delti, 1e-4, gompertz); 
 
     for(i=1; i <=NDIM; i++)
       for(j=i+1;j<=NDIM;j++)
@@ -5014,48 +5020,48 @@ Interval (in months) between two waves:
     for (i=1;i<=NDIM;i++) 
       printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i]));
 
-lsurv=vector(1,AGESUP);
+    lsurv=vector(1,AGESUP);
     lpop=vector(1,AGESUP);
     tpop=vector(1,AGESUP);
     lsurv[agegomp]=100000;
-   
-     for (k=agegomp;k<=AGESUP;k++) {
+    
+    for (k=agegomp;k<=AGESUP;k++) {
       agemortsup=k;
       if (p[1]*exp(p[2]*(k-agegomp))>1) break;
     }
-   
-      for (k=agegomp;k<agemortsup;k++)
+    
+    for (k=agegomp;k<agemortsup;k++)
       lsurv[k+1]=lsurv[k]-lsurv[k]*(p[1]*exp(p[2]*(k-agegomp)));
-
+    
     for (k=agegomp;k<agemortsup;k++){
       lpop[k]=(lsurv[k]+lsurv[k+1])/2.;
       sumlpop=sumlpop+lpop[k];
     }
-
- tpop[agegomp]=sumlpop;
+    
+    tpop[agegomp]=sumlpop;
     for (k=agegomp;k<(agemortsup-3);k++){
       /*  tpop[k+1]=2;*/
       tpop[k+1]=tpop[k]-lpop[k];
-       }
-   
-   
-       printf("\nAge   lx     qx    dx    Lx     Tx     e(x)\n");
+    }
+    
+    
+    printf("\nAge   lx     qx    dx    Lx     Tx     e(x)\n");
     for (k=agegomp;k<(agemortsup-2);k++) 
       printf("%d %.0lf %lf %.0lf %.0lf %.0lf %lf\n",k,lsurv[k],p[1]*exp(p[2]*(k-agegomp)),(p[1]*exp(p[2]*(k-agegomp)))*lsurv[k],lpop[k],tpop[k],tpop[k]/lsurv[k]);
-
-
+    
+    
     replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */
     printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
     
     printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \
 		     stepm, weightopt,\
 		     model,imx,p,matcov,agemortsup);
-
+    
     free_vector(lsurv,1,AGESUP);
     free_vector(lpop,1,AGESUP);
     free_vector(tpop,1,AGESUP);
   } /* Endof if mle==-3 */
-
+  
   else{ /* For mle >=1 */
   
     likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */