--- imach/src/imach.c	2014/09/02 09:24:05	1.160
+++ imach/src/imach.c	2014/12/22 11:40:47	1.166
@@ -1,6 +1,33 @@
-/* $Id: imach.c,v 1.160 2014/09/02 09:24:05 brouard Exp $
+/* $Id: imach.c,v 1.166 2014/12/22 11:40:47 brouard Exp $
   $State: Exp $
   $Log: imach.c,v $
+  Revision 1.166  2014/12/22 11:40:47  brouard
+  *** empty log message ***
+
+  Revision 1.165  2014/12/16 11:20:36  brouard
+  Summary: After compiling on Visual C
+
+  * imach.c (Module): Merging 1.61 to 1.162
+
+  Revision 1.164  2014/12/16 10:52:11  brouard
+  Summary: Merging with Visual C after suppressing some warnings for unused variables. Also fixing Saito's bug 0.98Xn
+
+  * imach.c (Module): Merging 1.61 to 1.162
+
+  Revision 1.163  2014/12/16 10:30:11  brouard
+  * imach.c (Module): Merging 1.61 to 1.162
+
+  Revision 1.162  2014/09/25 11:43:39  brouard
+  Summary: temporary backup 0.99!
+
+  Revision 1.1  2014/09/16 11:06:58  brouard
+  Summary: With some code (wrong) for nlopt
+
+  Author:
+
+  Revision 1.161  2014/09/15 20:41:41  brouard
+  Summary: Problem with macro SQR on Intel compiler
+
   Revision 1.160  2014/09/02 09:24:05  brouard
   *** empty log message ***
 
@@ -478,9 +505,8 @@
  end
 */
 
+#define POWELL /* Instead of NLOPT */
 
-
- 
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
@@ -512,6 +538,13 @@
 #include <gsl/gsl_multimin.h>
 #endif
 
+#ifdef NLOPT
+#include <nlopt.h>
+typedef struct {
+  double (* function)(double [] );
+} myfunc_data ;
+#endif
+
 /* #include <libintl.h> */
 /* #define _(String) gettext (String) */
 
@@ -536,7 +569,7 @@
 #define YEARM 12. /**< Number of months per year */
 #define AGESUP 130
 #define AGEBASE 40
-#define AGEGOMP 10. /**< Minimal age for Gompertz adjustment */
+#define AGEGOMP 10 /**< Minimal age for Gompertz adjustment */
 #ifdef _WIN32
 #define DIRSEPARATOR '\\'
 #define CHARSEPARATOR "\\"
@@ -547,11 +580,11 @@
 #define ODIRSEPARATOR '\\'
 #endif
 
-/* $Id: imach.c,v 1.160 2014/09/02 09:24:05 brouard Exp $ */
+/* $Id: imach.c,v 1.166 2014/12/22 11:40:47 brouard Exp $ */
 /* $State: Exp $ */
 
-char version[]="Imach version 0.98nX, August 2014,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121)";
-char fullversion[]="$Revision: 1.160 $ $Date: 2014/09/02 09:24:05 $"; 
+char version[]="Imach version 0.99, September 2014,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121)";
+char fullversion[]="$Revision: 1.166 $ $Date: 2014/12/22 11:40:47 $"; 
 char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];
 int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
@@ -581,6 +614,7 @@ int **mw; /* mw[mi][i] is number of the
 int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */
 int **bh; /* bh[mi][i] is the bias (+ or -) for this individual if the delay between
 	   * wave mi and wave mi+1 is not an exact multiple of stepm. */
+int countcallfunc=0;  /* Count the number of calls to func */
 double jmean=1; /* Mean space between 2 waves */
 double **matprod2(); /* test */
 double **oldm, **newm, **savm; /* Working pointers to matrices */
@@ -664,7 +698,12 @@ static double maxarg1,maxarg2;
   
 #define SIGN(a,b) ((b)>0.0 ? fabs(a) : -fabs(a))
 #define rint(a) floor(a+0.5)
-
+/* http://www.thphys.uni-heidelberg.de/~robbers/cmbeasy/doc/html/myutils_8h-source.html */
+/* #define mytinydouble 1.0e-16 */
+/* #define DEQUAL(a,b) (fabs((a)-(b))<mytinydouble) */
+/* http://www.thphys.uni-heidelberg.de/~robbers/cmbeasy/doc/html/mynrutils_8h-source.html */
+/* static double dsqrarg; */
+/* #define DSQR(a) (DEQUAL((dsqrarg=(a)),0.0) ? 0.0 : dsqrarg*dsqrarg) */
 static double sqrarg;
 #define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 :sqrarg*sqrarg)
 #define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;} 
@@ -1090,6 +1129,19 @@ char *subdirf3(char fileres[], char *pre
   return tmpout;
 }
 
+char *asc_diff_time(long time_sec, char ascdiff[])
+{
+  long sec_left, days, hours, minutes;
+  days = (time_sec) / (60*60*24);
+  sec_left = (time_sec) % (60*60*24);
+  hours = (sec_left) / (60*60) ;
+  sec_left = (sec_left) %(60*60);
+  minutes = (sec_left) /60;
+  sec_left = (sec_left) % (60);
+  sprintf(ascdiff,"%ld day(s) %ld hour(s) %ld minute(s) %ld second(s)",days, hours, minutes, sec_left);  
+  return ascdiff;
+}
+
 /***************** f1dim *************************/
 extern int ncom; 
 extern double *pcom,*xicom;
@@ -1114,7 +1166,7 @@ double brent(double ax, double bx, doubl
   int iter; 
   double a,b,d,etemp;
   double fu=0,fv,fw,fx;
-  double ftemp;
+  double ftemp=0.;
   double p,q,r,tol1,tol2,u,v,w,x,xm; 
   double e=0.0; 
  
@@ -1128,7 +1180,7 @@ double brent(double ax, double bx, doubl
     /*		if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret)))*/
     printf(".");fflush(stdout);
     fprintf(ficlog,".");fflush(ficlog);
-#ifdef DEBUG
+#ifdef DEBUGBRENT
     printf("br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol);
     fprintf(ficlog,"br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol);
     /*		if ((fabs(x-xm) <= (tol2-0.5*(b-a)))||(2.0*fabs(fu-ftemp) <= ftol*1.e-2*(fabs(fu)+fabs(ftemp)))) { */
@@ -1198,21 +1250,29 @@ void mnbrak(double *ax, double *bx, doub
       } 
   *cx=(*bx)+GOLD*(*bx-*ax); 
   *fc=(*func)(*cx); 
-  while (*fb > *fc) { 
+  while (*fb > *fc) { /* Declining fa, fb, fc */
     r=(*bx-*ax)*(*fb-*fc); 
     q=(*bx-*cx)*(*fb-*fa); 
     u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/ 
-      (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r)); 
-    ulim=(*bx)+GLIMIT*(*cx-*bx); 
-    if ((*bx-u)*(u-*cx) > 0.0) { 
+      (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r)); /* Minimum abscisse of a parabolic estimated from (a,fa), (b,fb) and (c,fc). */
+    ulim=(*bx)+GLIMIT*(*cx-*bx); /* Maximum abscisse where function can be evaluated */
+    if ((*bx-u)*(u-*cx) > 0.0) { /* if u between b and c */
       fu=(*func)(u); 
-    } else if ((*cx-u)*(u-ulim) > 0.0) { 
+#ifdef DEBUG
+      /* f(x)=A(x-u)**2+f(u) */
+      double A, fparabu; 
+      A= (*fb - *fa)/(*bx-*ax)/(*bx+*ax-2*u);
+      fparabu= *fa - A*(*ax-u)*(*ax-u);
+      printf("mnbrak (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf),  (*u=%.12f, fu=%.12lf, fparabu=%.12f)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu, fparabu);
+      fprintf(ficlog, "mnbrak (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf),  (*u=%.12f, fu=%.12lf, fparabu=%.12f)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu, fparabu);
+#endif 
+    } else if ((*cx-u)*(u-ulim) > 0.0) { /* u is after c but before ulim */
       fu=(*func)(u); 
       if (fu < *fc) { 
 	SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx)) 
 	  SHFT(*fb,*fc,fu,(*func)(u)) 
 	  } 
-    } else if ((u-ulim)*(ulim-*cx) >= 0.0) { 
+    } else if ((u-ulim)*(ulim-*cx) >= 0.0) { /* u outside ulim (verifying that ulim is beyond c) */
       u=ulim; 
       fu=(*func)(u); 
     } else { 
@@ -1225,7 +1285,11 @@ void mnbrak(double *ax, double *bx, doub
 } 
 
 /*************** linmin ************************/
-
+/* Given an n -dimensional point p[1..n] and an n -dimensional direction xi[1..n] , moves and
+resets p to where the function func(p) takes on a minimum along the direction xi from p ,
+and replaces xi by the actual vector displacement that p was moved. Also returns as fret
+the value of func at the returned location p . This is actually all accomplished by calling the
+routines mnbrak and brent .*/
 int ncom; 
 double *pcom,*xicom;
 double (*nrfunc)(double []); 
@@ -1251,8 +1315,8 @@ void linmin(double p[], double xi[], int
   } 
   ax=0.0; 
   xx=1.0; 
-  mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim); 
-  *fret=brent(ax,xx,bx,f1dim,TOL,&xmin); 
+  mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim); /* Find a bracket a,x,b in direction n=xi ie xicom */
+  *fret=brent(ax,xx,bx,f1dim,TOL,&xmin); /* Find a minimum P+lambda n in that direction (lambdamin), with TOL between abscisses */
 #ifdef DEBUG
   printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);
   fprintf(ficlog,"retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);
@@ -1265,20 +1329,16 @@ void linmin(double p[], double xi[], int
   free_vector(pcom,1,n); 
 } 
 
-char *asc_diff_time(long time_sec, char ascdiff[])
-{
-  long sec_left, days, hours, minutes;
-  days = (time_sec) / (60*60*24);
-  sec_left = (time_sec) % (60*60*24);
-  hours = (sec_left) / (60*60) ;
-  sec_left = (sec_left) %(60*60);
-  minutes = (sec_left) /60;
-  sec_left = (sec_left) % (60);
-  sprintf(ascdiff,"%ld day(s) %ld hour(s) %ld minute(s) %ld second(s)",days, hours, minutes, sec_left);  
-  return ascdiff;
-}
 
 /*************** powell ************************/
+/*
+Minimization of a function func of n variables. Input consists of an initial starting point
+p[1..n] ; an initial matrix xi[1..n][1..n] , whose columns contain the initial set of di-
+rections (usually the n unit vectors); and ftol , the fractional tolerance in the function value
+such that failure to decrease by more than this amount on one iteration signals doneness. On
+output, p is set to the best point found, xi is the then-current direction set, fret is the returned
+function value at p , and iter is the number of iterations taken. The routine linmin is used.
+ */
 void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret, 
 	    double (*func)(double [])) 
 { 
@@ -1319,17 +1379,15 @@ void powell(double p[], double **xi, int
     if(*iter <=3){
       tml = *localtime(&rcurr_time);
       strcpy(strcurr,asctime(&tml));
-/*       asctime_r(&tm,strcurr); */
       rforecast_time=rcurr_time; 
       itmp = strlen(strcurr);
       if(strcurr[itmp-1]=='\n')  /* Windows outputs with a new line */
 	strcurr[itmp-1]='\0';
-      printf("\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);
+      printf("\nConsidering the time needed for the last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);
       fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);
       for(niterf=10;niterf<=30;niterf+=10){
 	rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time);
 	forecast_time = *localtime(&rforecast_time);
-/* 	asctime_r(&tmf,strfor); */
 	strcpy(strfor,asctime(&forecast_time));
 	itmp = strlen(strfor);
 	if(strfor[itmp-1]=='\n')
@@ -1342,8 +1400,8 @@ void powell(double p[], double **xi, int
       for (j=1;j<=n;j++) xit[j]=xi[j][i]; 
       fptt=(*fret); 
 #ifdef DEBUG
-      printf("fret=%lf \n",*fret);
-      fprintf(ficlog,"fret=%lf \n",*fret);
+	  printf("fret=%lf, %lf, %lf \n", *fret, *fret, *fret);
+	  fprintf(ficlog, "fret=%lf, %lf, %lf \n", *fret, *fret, *fret);
 #endif
       printf("%d",i);fflush(stdout);
       fprintf(ficlog,"%d",i);fflush(ficlog);
@@ -1361,13 +1419,13 @@ void powell(double p[], double **xi, int
 	fprintf(ficlog," x(%d)=%.12e",j,xit[j]);
       }
       for(j=1;j<=n;j++) {
-	printf(" p=%.12e",p[j]);
-	fprintf(ficlog," p=%.12e",p[j]);
+	printf(" p(%d)=%.12e",j,p[j]);
+	fprintf(ficlog," p(%d)=%.12e",j,p[j]);
       }
       printf("\n");
       fprintf(ficlog,"\n");
 #endif
-    } 
+    } /* end i */
     if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) {
 #ifdef DEBUG
       int k[2],l;
@@ -1400,20 +1458,45 @@ void powell(double p[], double **xi, int
       return; 
     } 
     if (*iter == ITMAX) nrerror("powell exceeding maximum iterations."); 
-    for (j=1;j<=n;j++) { 
+    for (j=1;j<=n;j++) { /* Computes an extrapolated point */
       ptt[j]=2.0*p[j]-pt[j]; 
       xit[j]=p[j]-pt[j]; 
       pt[j]=p[j]; 
     } 
     fptt=(*func)(ptt); 
-    if (fptt < fp) { 
-      t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); 
-      if (t < 0.0) { 
-	linmin(p,xit,n,fret,func); 
+    if (fptt < fp) { /* If extrapolated point is better, decide if we keep that new direction or not */
+      /* (x1 f1=fp), (x2 f2=*fret), (x3 f3=fptt), (xm fm) */
+      /* From x1 (P0) distance of x2 is at h and x3 is 2h */
+      /* Let f"(x2) be the 2nd derivative equal everywhere.  */
+      /* Then the parabolic through (x1,f1), (x2,f2) and (x3,f3) */
+      /* will reach at f3 = fm + h^2/2 f"m  ; f" = (f1 -2f2 +f3 ) / h**2 */
+      /* f1-f3 = delta(2h) = 2 h**2 f'' = 2(f1- 2f2 +f3) */
+      /* Thus we compare delta(2h) with observed f1-f3 */
+      /* or best gain on one ancient line 'del' with total  */
+      /* gain f1-f2 = f1 - f2 - 'del' with del  */
+      /* t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); */
+
+      t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del);
+      t= t- del*SQR(fp-fptt);
+      printf("t1= %.12lf, t2= %.12lf, t=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t);
+      fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t);
+#ifdef DEBUG
+      printf("t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt),
+	     (fp-(*fret)-del)*(fp-(*fret)-del),(fp-fptt)*(fp-fptt));
+      fprintf(ficlog,"t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt),
+	     (fp-(*fret)-del)*(fp-(*fret)-del),(fp-fptt)*(fp-fptt));
+      printf("tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t);
+      fprintf(ficlog, "tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t);
+#endif
+      if (t < 0.0) { /* Then we use it for last direction */
+	linmin(p,xit,n,fret,func); /* computes mean on the extrapolated direction.*/
 	for (j=1;j<=n;j++) { 
-	  xi[j][ibig]=xi[j][n]; 
-	  xi[j][n]=xit[j]; 
+	  xi[j][ibig]=xi[j][n]; /* Replace the direction with biggest decrease by n */
+	  xi[j][n]=xit[j];      /* and nth direction by the extrapolated */
 	}
+	printf("Gaining to use average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
+	fprintf(ficlog,"Gaining to use average direction of P0 P%d instead of biggest increase direction :\n",n,ibig);
+
 #ifdef DEBUG
 	printf("Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
 	fprintf(ficlog,"Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
@@ -1424,8 +1507,8 @@ void powell(double p[], double **xi, int
 	printf("\n");
 	fprintf(ficlog,"\n");
 #endif
-      }
-    } 
+      } /* end of t negative */
+    } /* end if (fptt < fp)  */
   } 
 } 
 
@@ -1513,7 +1596,7 @@ double **pmij(double **ps, double *cov,
   */
   double s1, lnpijopii;
   /*double t34;*/
-  int i,j,j1, nc, ii, jj;
+  int i,j, nc, ii, jj;
 
     for(i=1; i<= nlstate; i++){
       for(j=1; j<i;j++){
@@ -1655,6 +1738,25 @@ double ***hpxij(double ***po, int nhstep
   return po;
 }
 
+#ifdef NLOPT
+  double  myfunc(unsigned n, const double *p1, double *grad, void *pd){
+  double fret;
+  double *xt;
+  int j;
+  myfunc_data *d2 = (myfunc_data *) pd;
+/* xt = (p1-1); */
+  xt=vector(1,n); 
+  for (j=1;j<=n;j++)   xt[j]=p1[j-1]; /* xt[1]=p1[0] */
+
+  fret=(d2->function)(xt); /*  p xt[1]@8 is fine */
+  /* fret=(*func)(xt); /\*  p xt[1]@8 is fine *\/ */
+  printf("Function = %.12lf ",fret);
+  for (j=1;j<=n;j++) printf(" %d %.8lf", j, xt[j]); 
+  printf("\n");
+ free_vector(xt,1,n);
+  return fret;
+}
+#endif
 
 /*************** log-likelihood *************/
 double func( double *x)
@@ -1673,6 +1775,9 @@ double func( double *x)
   /*for(i=1;i<imx;i++) 
     printf(" %d\n",s[4][i]);
   */
+
+  ++countcallfunc;
+
   cov[1]=1.;
 
   for(k=1; k<=nlstate; k++) ll[k]=0.;
@@ -2054,11 +2159,23 @@ void likelione(FILE *ficres,double p[],
 
 void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, double ftol, double (*func)(double []))
 {
-  int i,j, iter;
+  int i,j, iter=0;
   double **xi;
   double fret;
   double fretone; /* Only one call to likelihood */
   /*  char filerespow[FILENAMELENGTH];*/
+
+#ifdef NLOPT
+  int creturn;
+  nlopt_opt opt;
+  /* double lb[9] = { -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL }; /\* lower bounds *\/ */
+  double *lb;
+  double minf; /* the minimum objective value, upon return */
+  double * p1; /* Shifted parameters from 0 instead of 1 */
+  myfunc_data dinst, *d = &dinst;
+#endif
+
+
   xi=matrix(1,npar,1,npar);
   for (i=1;i<=npar;i++)
     for (j=1;j<=npar;j++)
@@ -2075,14 +2192,41 @@ void mlikeli(FILE *ficres,double p[], in
     for(j=1;j<=nlstate+ndeath;j++)
       if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);
   fprintf(ficrespow,"\n");
-
+#ifdef POWELL
   powell(p,xi,npar,ftol,&iter,&fret,func);
+#endif
 
+#ifdef NLOPT
+#ifdef NEWUOA
+  opt = nlopt_create(NLOPT_LN_NEWUOA,npar);
+#else
+  opt = nlopt_create(NLOPT_LN_BOBYQA,npar);
+#endif
+  lb=vector(0,npar-1);
+  for (i=0;i<npar;i++) lb[i]= -HUGE_VAL;
+  nlopt_set_lower_bounds(opt, lb);
+  nlopt_set_initial_step1(opt, 0.1);
+  
+  p1= (p+1); /*  p *(p+1)@8 and p *(p1)@8 are equal p1[0]=p[1] */
+  d->function = func;
+  printf(" Func %.12lf \n",myfunc(npar,p1,NULL,d));
+  nlopt_set_min_objective(opt, myfunc, d);
+  nlopt_set_xtol_rel(opt, ftol);
+  if ((creturn=nlopt_optimize(opt, p1, &minf)) < 0) {
+    printf("nlopt failed! %d\n",creturn); 
+  }
+  else {
+    printf("found minimum after %d evaluations (NLOPT=%d)\n", countcallfunc ,NLOPT);
+    printf("found minimum at f(%g,%g) = %0.10g\n", p[0], p[1], minf);
+    iter=1; /* not equal */
+  }
+  nlopt_destroy(opt);
+#endif
   free_matrix(xi,1,npar,1,npar);
   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));
+  printf("\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
+  fprintf(ficlog,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
+  fprintf(ficres,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
 
 }
 
@@ -2091,7 +2235,7 @@ void hesscov(double **matcov, double p[]
 {
   double  **a,**y,*x,pd;
   double **hess;
-  int i, j,jk;
+  int i, j;
   int *indx;
 
   double hessii(double p[], double delta, int theta, double delti[],double (*func)(double []),int npar);
@@ -2225,7 +2369,7 @@ double hessii(double x[], double delta,
 	k=kmax;
       }
       else if((k1 >khi/nkhif) || (k2 >khi/nkhif)){ /* Keeps lastvalue before 3.84/2 KHI2 5% 1d.f. */
-	k=kmax; l=lmax*10.;
+	k=kmax; l=lmax*10;
       }
       else if((k1 >khi/nkhi) || (k2 >khi/nkhi)){ 
 	delts=delt;
@@ -2240,7 +2384,7 @@ double hessii(double x[], double delta,
 double hessij( double x[], double delti[], int thetai,int thetaj,double (*func)(double []),int npar)
 {
   int i;
-  int l=1, l1, lmax=20;
+  int l=1, lmax=20;
   double k1,k2,k3,k4,res,fx;
   double p2[MAXPARM+1];
   int k;
@@ -2355,7 +2499,7 @@ void pstamp(FILE *fichier)
 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, char strstart[])
 {  /* Some frequencies */
   
-  int i, m, jk, k1,i1, j1, bool, z1,j;
+  int i, m, jk, j1, bool, z1,j;
   int first;
   double ***freq; /* Frequencies */
   double *pp, **prop;
@@ -2536,10 +2680,10 @@ void prevalence(double ***probs, double
      We still use firstpass and lastpass as another selection.
   */
  
-  int i, m, jk, k1, i1, j1, bool, z1,j;
-  double ***freq; /* Frequencies */
-  double *pp, **prop;
-  double pos,posprop; 
+  int i, m, jk, j1, bool, z1,j;
+
+  double **prop;
+  double posprop; 
   double  y2; /* in fractional years */
   int iagemin, iagemax;
   int first; /** to stop verbosity which is redirected to log file */
@@ -2630,7 +2774,7 @@ void  concatwav(int wav[], int **dh, int
   int j, k=0,jk, ju, jl;
   double sum=0.;
   first=0;
-  jmin=1e+5;
+  jmin=100000;
   jmax=-1;
   jmean=0.;
   for(i=1; i<=imx; i++){
@@ -2870,7 +3014,7 @@ void evsij(double ***eij, double x[], in
 
 {
   /* Health expectancies, no variances */
-  int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2;
+  int i, j, nhstepm, hstepm, h, nstepm;
   int nhstepma, nstepma; /* Decreasing with age */
   double age, agelim, hf;
   double ***p3mat;
@@ -3192,7 +3336,7 @@ void varevsij(char optionfilefiname[], d
   double **dnewm,**doldm;
   double **dnewmp,**doldmp;
   int i, j, nhstepm, hstepm, h, nstepm ;
-  int k, cptcode;
+  int k;
   double *xp;
   double **gp, **gm;  /* for var eij */
   double ***gradg, ***trgradg; /*for var eij */
@@ -3492,10 +3636,9 @@ void varprevlim(char fileres[], double *
 {
   /* Variance of prevalence limit */
   /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/
-  double **newm;
+
   double **dnewm,**doldm;
   int i, j, nhstepm, hstepm;
-  int k, cptcode;
   double *xp;
   double *gp, *gm;
   double **gradg, **trgradg;
@@ -3574,9 +3717,9 @@ void varprevlim(char fileres[], double *
 /************ Variance of one-step probabilities  ******************/
 void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax, char strstart[])
 {
-  int i, j=0,  i1, k1, l1, t, tj;
+  int i, j=0,  k1, l1, tj;
   int k2, l2, j1,  z1;
-  int k=0,l, cptcode;
+  int k=0, l;
   int first=1, first1, first2;
   double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp;
   double **dnewm,**doldm;
@@ -3584,7 +3727,7 @@ void varprob(char optionfilefiname[], do
   double *gp, *gm;
   double **gradg, **trgradg;
   double **mu;
-  double age,agelim, cov[NCOVMAX+1];
+  double age, cov[NCOVMAX+1];
   double std=2.0; /* Number of standard deviation wide of confidence ellipsoids */
   int theta;
   char fileresprob[FILENAMELENGTH];
@@ -3956,8 +4099,8 @@ fprintf(fichtm," \n<ul><li><b>Graphs</b>
 <img src=\"%s%d_2.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); 
        /* Period (stable) prevalence in each health state */
        for(cpt=1; cpt<=nlstate;cpt++){
-	 fprintf(fichtm,"<br>- Convergence from each state (1 to %d) to period (stable) prevalence in state %d <a href=\"%s%d_%d.png\">%s%d_%d.png</a><br> \
-<img src=\"%s%d_%d.png\">",nlstate, cpt, subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1);
+	 fprintf(fichtm,"<br>- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.png\">%s%d_%d.png</a><br> \
+<img src=\"%s%d_%d.png\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1);
        }
      for(cpt=1; cpt<=nlstate;cpt++) {
         fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) : <a href=\"%s%d%d.png\">%s%d%d.png</a> <br> \
@@ -4042,7 +4185,7 @@ true period expectancies (those weighted
 void printinggnuplot(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
 
   char dirfileres[132],optfileres[132];
-  int m0,cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0;
+  int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0;
   int ng=0;
 /*   if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */
 /*     printf("Problem with file %s",optionfilegnuplot); */
@@ -4228,7 +4371,7 @@ plot [%.f:%.f]  ", ageminpar, agemaxpar)
        } /* end k2 */
      } /* end jk */
    } /* end ng */
- avoid:
+ /* avoid: */
    fflush(ficgp); 
 }  /* end gnuplot */
 
@@ -4282,8 +4425,7 @@ prevforecast(char fileres[], double anpr
      dateprev1 dateprev2 range of dates during which prevalence is computed
      anproj2 year of en of projection (same day and month as proj1).
   */
-  int yearp, stepsize, hstepm, nhstepm, j, k, c, cptcod, i, h, i1;
-  int *popage;
+  int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1;
   double agec; /* generic age */
   double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
   double *popeffectif,*popcount;
@@ -4577,8 +4719,8 @@ void prwizard(int ncovmodel, int nlstate
 
   /* Wizard to print covariance matrix template */
 
-  char ca[32], cb[32], cc[32];
-  int i,j, k, l, li, lj, lk, ll, jj, npar, itimes;
+  char ca[32], cb[32];
+  int i,j, k, li, lj, lk, ll, jj, npar, itimes;
   int numlinepar;
 
   printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
@@ -4799,7 +4941,7 @@ fprintf(fichtm,"<ul><li><h4>Life table</
 void printinggnuplotmort(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){
 
   char dirfileres[132],optfileres[132];
-  int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;
+
   int ng;
 
 
@@ -4824,10 +4966,10 @@ int readdata(char datafile[], int firsto
   /*-------- data file ----------*/
   FILE *fic;
   char dummy[]="                         ";
-  int i, j, n;
+  int i=0, j=0, n=0;
   int linei, month, year,iout;
   char line[MAXLINE], linetmp[MAXLINE];
-  char stra[80], strb[80];
+  char stra[MAXLINE], strb[MAXLINE];
   char *stratrunc;
   int lstra;
 
@@ -4855,9 +4997,7 @@ int readdata(char datafile[], int firsto
       continue;
     }
     trimbb(linetmp,line); /* Trims multiple blanks in line */
-    for (j=0; line[j]!='\0';j++){
-      line[j]=linetmp[j];
-    }
+    strcpy(line, linetmp);
   
 
     for (j=maxwav;j>=1;j--){
@@ -4996,7 +5136,7 @@ int readdata(char datafile[], int firsto
   fclose(fic);
  
   return (0);
-  endread:
+  /* endread: */
     printf("Exiting readdata: ");
     fclose(fic);
     return (1);
@@ -5027,7 +5167,7 @@ int decodemodel ( char model[], int last
  */
 {
   int i, j, k, ks;
-  int i1, j1, k1, k2;
+  int  j1, k1, k2;
   char modelsav[80];
   char stra[80], strb[80], strc[80], strd[80],stre[80];
 
@@ -5176,7 +5316,7 @@ int decodemodel ( char model[], int last
 
 
   return (0); /* with covar[new additional covariate if product] and Tage if age */ 
-  endread:
+  /*endread:*/
     printf("Exiting decodemodel: ");
     return (1);
 }
@@ -5273,7 +5413,7 @@ calandcheckages(int imx, int maxwav, dou
   fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, *agemin, *agemax); 
 
   return (0);
-  endread:
+ /* endread:*/
     printf("Exiting calandcheckages: ");
     return (1);
 }
@@ -5293,38 +5433,35 @@ int main(int argc, char *argv[])
   double ssval;
 #endif
   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, month, year,iout;
-  int jj, ll, li, lj, lk, imk;
+  int i,j, k, n=MAXN,iter=0,m,size=100, cptcod;
+
+  int jj, ll, li, lj, lk;
   int numlinepar=0; /* Current linenumber of parameter file */
   int itimes;
   int NDIM=2;
   int vpopbased=0;
 
-  char ca[32], cb[32], cc[32];
+  char ca[32], cb[32];
   /*  FILE *fichtm; *//* Html File */
   /* FILE *ficgp;*/ /*Gnuplot File */
   struct stat info;
-  double agedeb, agefin,hf;
+  double agedeb;
   double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;
 
   double fret;
-  double **xi,tmp,delta;
-
   double dum; /* Dummy variable */
   double ***p3mat;
   double ***mobaverage;
-  int *indx;
-  char line[MAXLINE], linepar[MAXLINE];
+
+  char line[MAXLINE];
   char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE];
   char pathr[MAXLINE], pathimach[MAXLINE]; 
-  char **bp, *tok, *val; /* pathtot */
+  char *tok, *val; /* pathtot */
   int firstobs=1, lastobs=10;
-  int sdeb, sfin; /* Status at beginning and end */
-  int c,  h , cpt,l;
-  int ju,jl, mi;
-  int i1,j1, jk,aa,bb, stepsize, ij;
-  int jnais,jdc,jint4,jint1,jint2,jint3,*tab; 
+  int c,  h , cpt;
+  int jl;
+  int i1, j1, jk, stepsize;
+  int *tab; 
   int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */
   int mobilav=0,popforecast=0;
   int hstepm, nhstepm;
@@ -5333,7 +5470,7 @@ int main(int argc, char *argv[])
   double jprev1=1, mprev1=1,anprev1=2000,jprev2=1, mprev2=1,anprev2=2000;
   double jpyram=1, mpyram=1,anpyram=2000,jpyram1=1, mpyram1=1,anpyram1=2000;
 
-  double bage, fage, age, agelim, agebase;
+  double bage=0, fage=110, age, agelim, agebase;
   double ftolpl=FTOL;
   double **prlim;
   double ***param; /* Matrix of parameters */
@@ -5344,19 +5481,18 @@ int main(int argc, char *argv[])
   double ***eij, ***vareij;
   double **varpl; /* Variances of prevalence limits by age */
   double *epj, vepp;
-  double kk1, kk2;
+
   double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;
   double **ximort;
   char *alph[]={"a","a","b","c","d","e"}, str[4]="1234";
   int *dcwave;
 
-  char z[1]="c", occ;
+  char z[1]="c";
 
   /*char  *strt;*/
   char strtend[80];
 
-  long total_usecs;
- 
+
 /*   setlocale (LC_ALL, ""); */
 /*   bindtextdomain (PACKAGE, LOCALEDIR); */
 /*   textdomain (PACKAGE); */
@@ -5627,7 +5763,7 @@ run imach with mle=-1 to get a correct t
     for(i=1; i <=nlstate; i++){
       for(j=1; j <=nlstate+ndeath-1; j++){
 	fscanf(ficpar,"%1d%1d",&i1,&j1);
-	if ((i1-i)*(j1-j)!=0){
+	if ( (i1-i) * (j1-j) != 0){
 	  printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1);
 	  exit(1);
 	}
@@ -5996,7 +6132,7 @@ Interval (in months) between two waves:
     
 #ifdef GSL
     printf("GSL optimization\n");  fprintf(ficlog,"Powell\n");
-#elsedef
+#else
     printf("Powell\n");  fprintf(ficlog,"Powell\n");
 #endif
     strcpy(filerespow,"pow-mort"); 
@@ -6007,7 +6143,7 @@ Interval (in months) between two waves:
     }
 #ifdef GSL
     fprintf(ficrespow,"# GSL optimization\n# iter -2*LL");
-#elsedef
+#else
     fprintf(ficrespow,"# Powell\n# iter -2*LL");
 #endif
     /*  for (i=1;i<=nlstate;i++)
@@ -6569,7 +6705,8 @@ Interval (in months) between two waves:
 
 	for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
 	  oldm=oldms;savm=savms; /* Segmentation fault */
-	  varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart);
+	  cptcod= 0; /* To be deleted */
+	  varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */
 	  fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n#  (weighted average of eij where weights are ");
 	  if(vpopbased==1)
 	    fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
@@ -6661,7 +6798,7 @@ Interval (in months) between two waves:
     if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
     free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
   }  /* mle==-3 arrives here for freeing */
- endfree:
+ /* endfree:*/
     free_matrix(prlim,1,nlstate,1,nlstate); /*here or after loop ? */
     free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);
     free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);