--- imach/src/imach.c	2016/09/02 07:25:01	1.245
+++ imach/src/imach.c	2016/09/07 14:10:18	1.248
@@ -1,6 +1,15 @@
-/* $Id: imach.c,v 1.245 2016/09/02 07:25:01 brouard Exp $
+/* $Id: imach.c,v 1.248 2016/09/07 14:10:18 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.248  2016/09/07 14:10:18  brouard
+  *** empty log message ***
+
+  Revision 1.247  2016/09/02 11:11:21  brouard
+  *** empty log message ***
+
+  Revision 1.246  2016/09/02 08:49:22  brouard
+  *** empty log message ***
+
   Revision 1.245  2016/09/02 07:25:01  brouard
   *** empty log message ***
 
@@ -937,12 +946,12 @@ typedef struct {
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.245 2016/09/02 07:25:01 brouard Exp $ */
+/* $Id: imach.c,v 1.248 2016/09/07 14:10:18 brouard Exp $ */
 /* $State: Exp $ */
 #include "version.h"
 char version[]=__IMACH_VERSION__;
 char copyright[]="February 2016,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018";
-char fullversion[]="$Revision: 1.245 $ $Date: 2016/09/02 07:25:01 $"; 
+char fullversion[]="$Revision: 1.248 $ $Date: 2016/09/07 14:10:18 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -2581,6 +2590,7 @@ Earliest age to start was %d-%d=%d, ncvl
   /* If we start from prlim again, prlim tends to a constant matrix */
 
   int i, ii,j,k;
+  int first=0;
   double *min, *max, *meandiff, maxmax,sumnew=0.;
   /* double **matprod2(); */ /* test */
   double **out, cov[NCOVMAX+1], **bmij();
@@ -2706,7 +2716,12 @@ Earliest age to start was %d-%d=%d, ncvl
     }
   } /* age loop */
     /* After some age loop it doesn't converge */
-  printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\
+  if(first){
+    first=1;
+    printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. Others in log file only...\n\
+Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);
+  }
+  fprintf(ficlog,"Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\
 Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);
   /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */
   free_vector(min,1,nlstate);
@@ -3593,10 +3608,11 @@ double funcone( double *x)
       agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */
       ageend=agev[mw[mi][i]][i] + (dh[mi][i])*stepm/YEARM; /* Age at end of effective wave and at the end of transition */
       for(d=0; d<dh[mi][i]; d++){  /* Delay between two effective waves */
+      /* for(d=0; d<=0; d++){  /\* Delay between two effective waves Only one matrix to speed up*\/ */
 	/*dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]
 	  and mw[mi+1][i]. dh depends on stepm.*/
 	newm=savm;
-	agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
+	agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;  /* Here d is needed */
 	cov[2]=agexact;
 	if(nagesqr==1)
 	  cov[3]= agexact*agexact;
@@ -3649,7 +3665,7 @@ double funcone( double *x)
       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]); */
       if(globpr){
-	fprintf(ficresilk,"%9ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\
+	fprintf(ficresilk,"%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\
  %11.6f %11.6f %11.6f ", \
 		num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,
 		2*weight[i]*lli,out[s1][s2],savm[s1][s2]);
@@ -4007,8 +4023,8 @@ double hessij( double x[], double **hess
       kmax=kmax+10;
     }
     if(kmax >=10 || firstime ==1){
-      printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you may increase ftol=%.2e\n",thetai,thetaj, ftol);
-      fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you may increase ftol=%.2e\n",thetai,thetaj, ftol);
+      printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol);
+      fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol);
       printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e  res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
       fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e  res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
     }
@@ -6401,7 +6417,7 @@ void printinggnuplot(char fileresu[], ch
 	if(TKresult[nres]!= k1)
 	  continue;
 	/* We are interested in selected combination by the resultline */
-	printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt);
+	/* printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); */
 	fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files  and live state =%d ", cpt);
 	for (k=1; k<=cptcoveff; k++){    /* For each covariate k get corresponding value lv for combination k1 */
 	  lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */
@@ -6410,14 +6426,14 @@ void printinggnuplot(char fileresu[], ch
 	  /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
 	  vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */
 	  /* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */
-	  printf(" V%d=%d ",Tvaraff[k],vlv);
+	  /* printf(" V%d=%d ",Tvaraff[k],vlv); */
 	  fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
 	}
 	for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
-	  printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+	  /* printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); */
 	  fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
 	}	
-	printf("\n#\n");
+	/* printf("\n#\n"); */
 	fprintf(ficgp,"\n#\n");
 	if(invalidvarcomb[k1]){
 	  fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
@@ -8830,7 +8846,7 @@ Dummy[k] 0=dummy (0 1), 1 quantitative (
 }
 
 int calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn )
-{
+{/* Check ages at death */
   int i, m;
   int firstone=0;
   
@@ -9874,6 +9890,12 @@ int main(int argc, char *argv[])
   delti=delti3[1][1];
   /*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/
   if(mle==-1){ /* Print a wizard for help writing covariance matrix */
+/* We could also provide initial parameters values giving by simple logistic regression 
+ * only one way, that is without matrix product. We will have nlstate maximizations */
+      /* for(i=1;i<nlstate;i++){ */
+      /* 	/\*reducing xi for 1 to npar to 1 to ncovmodel; *\/ */
+      /*    mlikeli(ficres,p, ncovmodel, ncovmodel, nlstate, ftol, funcnoprod); */
+      /* } */
     prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
     printf(" You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso);
     fprintf(ficlog," You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso);
@@ -9882,6 +9904,11 @@ int main(int argc, char *argv[])
     fclose (ficlog);
     goto end;
     exit(0);
+  } else if(mle==-2) { /* Guessing from means */
+    prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
+    printf(" You chose mle=-2, look at file %s for a template of covariance matrix \n",filereso);
+    fprintf(ficlog," You chose mle=-2, look at file %s for a template of covariance matrix \n",filereso);
+   
   }  else if(mle==-5) { /* Main Wizard */
     prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
     printf(" You chose mle=-3, look at file %s for a template of covariance matrix \n",filereso);
@@ -10205,7 +10232,7 @@ Please run with mle=-1 to get a correct
   */
 
   concatwav(wav, dh, bh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);
-  /* */
+  /* Concatenates waves */
  
   free_vector(moisdc,1,n);
   free_vector(andc,1,n);
@@ -10408,9 +10435,9 @@ Interval (in months) between two waves:
   /* For mortality only */
   if (mle==-3){
     ximort=matrix(1,NDIM,1,NDIM); 
-		for(i=1;i<=NDIM;i++)
-			for(j=1;j<=NDIM;j++)
-				ximort[i][j]=0.;
+    for(i=1;i<=NDIM;i++)
+      for(j=1;j<=NDIM;j++)
+	ximort[i][j]=0.;
     /*     ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */
     cens=ivector(1,n);
     ageexmed=vector(1,n);
@@ -10646,6 +10673,10 @@ Please run with mle=-1 to get a correct
     printf("\n");
     if(mle>=1){ /* Could be 1 or 2, Real Maximization */
       /* mlikeli uses func not funcone */
+      /* for(i=1;i<nlstate;i++){ */
+      /* 	/\*reducing xi for 1 to npar to 1 to ncovmodel; *\/ */
+      /*    mlikeli(ficres,p, ncovmodel, ncovmodel, nlstate, ftol, funcnoprod); */
+      /* } */
       mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
     }
     if(mle==0) {/* No optimization, will print the likelihoods for the datafile */