imach/src/imach.c - diff

Return to imach.c CVS log

Up to [Local Repository] / imach / src

Diff for /imach/src/imach.c between versions 1.135 and 1.166

-version 1.135, 2009/10/29 15:33:14
+version 1.166, 2014/12/22 11:40:47
  Line 1
  /* $Id$
    $State$
    $Log$
+   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 ***
+   Revision 1.159  2014/09/01 10:34:10  brouard
+   Summary: WIN32
+   Author: Brouard
+   Revision 1.158  2014/08/27 17:11:51  brouard
+   *** empty log message ***
+   Revision 1.157  2014/08/27 16:26:55  brouard
+   Summary: Preparing windows Visual studio version
+   Author: Brouard
+   In order to compile on Visual studio, time.h is now correct and time_t
+   and tm struct should be used. difftime should be used but sometimes I
+   just make the differences in raw time format (time(&now).
+   Trying to suppress #ifdef LINUX
+   Add xdg-open for __linux in order to open default browser.
+   Revision 1.156  2014/08/25 20:10:10  brouard
+   *** empty log message ***
+   Revision 1.155  2014/08/25 18:32:34  brouard
+   Summary: New compile, minor changes
+   Author: Brouard
+   Revision 1.154  2014/06/20 17:32:08  brouard
+   Summary: Outputs now all graphs of convergence to period prevalence
+   Revision 1.153  2014/06/20 16:45:46  brouard
+   Summary: If 3 live state, convergence to period prevalence on same graph
+   Author: Brouard
+   Revision 1.152  2014/06/18 17:54:09  brouard
+   Summary: open browser, use gnuplot on same dir than imach if not found in the path
+   Revision 1.151  2014/06/18 16:43:30  brouard
+   *** empty log message ***
+   Revision 1.150  2014/06/18 16:42:35  brouard
+   Summary: If gnuplot is not in the path try on same directory than imach binary (OSX)
+   Author: brouard
+   Revision 1.149  2014/06/18 15:51:14  brouard
+   Summary: Some fixes in parameter files errors
+   Author: Nicolas Brouard
+   Revision 1.148  2014/06/17 17:38:48  brouard
+   Summary: Nothing new
+   Author: Brouard
+   Just a new packaging for OS/X version 0.98nS
+   Revision 1.147  2014/06/16 10:33:11  brouard
+   *** empty log message ***
+   Revision 1.146  2014/06/16 10:20:28  brouard
+   Summary: Merge
+   Author: Brouard
+   Merge, before building revised version.
+   Revision 1.145  2014/06/10 21:23:15  brouard
+   Summary: Debugging with valgrind
+   Author: Nicolas Brouard
+   Lot of changes in order to output the results with some covariates
+   After the Edimburgh REVES conference 2014, it seems mandatory to
+   improve the code.
+   No more memory valgrind error but a lot has to be done in order to
+   continue the work of splitting the code into subroutines.
+   Also, decodemodel has been improved. Tricode is still not
+   optimal. nbcode should be improved. Documentation has been added in
+   the source code.
+   Revision 1.143  2014/01/26 09:45:38  brouard
+   Summary: Version 0.98nR (to be improved, but gives same optimization results as 0.98k. Nice, promising
+   * imach.c (Module): Trying to merge old staffs together while being at Tokyo. Not tested...
+   (Module): Version 0.98nR Running ok, but output format still only works for three covariates.
+   Revision 1.142  2014/01/26 03:57:36  brouard
+   Summary: gnuplot changed plot w l 1 has to be changed to plot w l lt 2
+   * imach.c (Module): Trying to merge old staffs together while being at Tokyo. Not tested...
+   Revision 1.141  2014/01/26 02:42:01  brouard
+   * imach.c (Module): Trying to merge old staffs together while being at Tokyo. Not tested...
+   Revision 1.140  2011/09/02 10:37:54  brouard
+   Summary: times.h is ok with mingw32 now.
+   Revision 1.139  2010/06/14 07:50:17  brouard
+   After the theft of my laptop, I probably lost some lines of codes which were not uploaded to the CVS tree.
+   I remember having already fixed agemin agemax which are pointers now but not cvs saved.
+   Revision 1.138  2010/04/30 18:19:40  brouard
+   *** empty log message ***
+   Revision 1.137  2010/04/29 18:11:38  brouard
+   (Module): Checking covariates for more complex models
+   than V1+V2. A lot of change to be done. Unstable.
+   Revision 1.136  2010/04/26 20:30:53  brouard
+   (Module): merging some libgsl code. Fixing computation
+   of likelione (using inter/intrapolation if mle = 0) in order to
+   get same likelihood as if mle=1.
+   Some cleaning of code and comments added.
    Revision 1.135  2009/10/29 15:33:14  brouard
    (Module): Now imach stops if date of birth, at least year of birth, is not given. Some cleaning of the code.
- Line 340
+ Line 477
        begin-prev-date,...
    open gnuplot file
    open html file
-   period (stable) prevalence
+   period (stable) prevalence      | pl_nom    1-1 2-2 etc by covariate
-    for age prevalim()
+    for age prevalim()             | #****** V1=0  V2=1  V3=1  V4=0 ******
-   h Pij x
+                                   | 65 1 0 2 1 3 1 4 0  0.96326 0.03674
-   variance of p varprob
+     freexexit2 possible for memory heap.
+   h Pij x                         | pij_nom  ficrestpij
+    # Cov Agex agex+h hpijx with i,j= 1-1 1-2     1-3     2-1     2-2     2-3
+ 85   85    1.00000             0.00000 0.00000 0.00000 1.00000 0.00000
+ 85   86    0.68299             0.22291 0.09410 0.71093 0.00000 0.28907
+ 65   99    0.00364             0.00322 0.99314 0.00350 0.00310 0.99340
+ 65  100    0.00214             0.00204 0.99581 0.00206 0.00196 0.99597
+   variance of p one-step probabilities varprob  | prob_nom   ficresprob #One-step probabilities and stand. devi in ()
+    Standard deviation of one-step probabilities | probcor_nom   ficresprobcor #One-step probabilities and correlation matrix
+    Matrix of variance covariance of one-step probabilities |  probcov_nom ficresprobcov #One-step probabilities and covariance matrix
    forecasting if prevfcast==1 prevforecast call prevalence()
    health expectancies
    Variance-covariance of DFLE
- Line 356
+ Line 505
   end
  */
+ #define POWELL /* Instead of NLOPT */
  #include <math.h>
  #include <stdio.h>
  #include <stdlib.h>
  #include <string.h>
+ #ifdef _WIN32
+ #include <io.h>
+ #else
  #include <unistd.h>
+ #endif
  #include <limits.h>
  #include <sys/types.h>
  #include <sys/stat.h>
  #include <errno.h>
- extern int errno;
+ /* extern int errno; */
+ /* #ifdef LINUX */
+ /* #include <time.h> */
+ /* #include "timeval.h" */
+ /* #else */
  /* #include <sys/time.h> */
+ /* #endif */
  #include <time.h>
- #include "timeval.h"
+ #ifdef GSL
+ #include <gsl/gsl_errno.h>
+ #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) */
- #define MAXLINE 256
+ #define MAXLINE 1024 /* Was 256. Overflow with 312 with 2 states and 4 covariates. Should be ok */
  #define GNUPLOTPROGRAM "gnuplot"
  /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
- Line 387  extern int errno;
+ Line 557  extern int errno;
  #define GLOCK_ERROR_NOPATH              -1      /* empty path */
  #define GLOCK_ERROR_GETCWD              -2      /* cannot get cwd */
- #define MAXPARM 128 /* Maximum number of parameters for the optimization */
+ #define MAXPARM 128 /**< Maximum number of parameters for the optimization */
- #define NPARMAX 64 /* (nlstate+ndeath-1)*nlstate*ncovmodel */
+ #define NPARMAX 64 /**< (nlstate+ndeath-1)*nlstate*ncovmodel */
  #define NINTERVMAX 8
- #define NLSTATEMAX 8 /* Maximum number of live states (for func) */
+ #define NLSTATEMAX 8 /**< Maximum number of live states (for func) */
- #define NDEATHMAX 8 /* Maximum number of dead states (for func) */
+ #define NDEATHMAX 8 /**< Maximum number of dead states (for func) */
- #define NCOVMAX 20 /* Maximum number of covariates */
+ #define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */
+ #define codtabm(h,k)  1 & (h-1) >> (k-1) ;
  #define MAXN 20000
- #define YEARM 12. /* Number of months per year */
+ #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 UNIX
+ #ifdef _WIN32
- #define DIRSEPARATOR '/'
- #define CHARSEPARATOR "/"
- #define ODIRSEPARATOR '\\'
- #else
  #define DIRSEPARATOR '\\'
  #define CHARSEPARATOR "\\"
  #define ODIRSEPARATOR '/'
+ #else
+ #define DIRSEPARATOR '/'
+ #define CHARSEPARATOR "/"
+ #define ODIRSEPARATOR '\\'
  #endif
  /* $Id$ */
  /* $State$ */
- char version[]="Imach version 0.98l, October 2009, INED-EUROREVES-Institut de longevite ";
+ 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$ $Date$";
  char strstart[80];
  char optionfilext[10], optionfilefiname[FILENAMELENGTH];
  int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings  */
  int nvar=0, nforce=0; /* Number of variables, number of forces */
- int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov=0; /* Number of covariates, of covariates with '*age' */
+ /* Number of covariates model=V2+V1+ V3*age+V2*V4 */
+ int cptcovn=0; /**< cptcovn number of covariates added in the model (excepting constant and age and age*product) */
+ int cptcovt=0; /**< cptcovt number of covariates added in the model (excepting constant and age) */
+ int cptcovs=0; /**< cptcovs number of simple covariates V2+V1 =2 */
+ int cptcovage=0; /**< Number of covariates with age: V3*age only =1 */
+ int cptcovprodnoage=0; /**< Number of covariate products without age */
+ int cptcoveff=0; /* Total number of covariates to vary for printing results */
+ int cptcov=0; /* Working variable */
  int npar=NPARMAX;
  int nlstate=2; /* Number of live states */
  int ndeath=1; /* Number of dead states */
- Line 436  int **mw; /* mw[mi][i] is number of the
+ Line 614  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 */
  double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
- FILE *fic,*ficpar, *ficparo,*ficres, *ficresp, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
+ /*FILE *fic ; */ /* Used in readdata only */
+ FILE *ficpar, *ficparo,*ficres, *ficresp, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;
  FILE *ficlog, *ficrespow;
  int globpr=0; /* Global variable for printing or not */
  double fretone; /* Only one call to likelihood */
- Line 477  char popfile[FILENAMELENGTH];
+ Line 658  char popfile[FILENAMELENGTH];
  char optionfilegnuplot[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH], optionfilehtmcov[FILENAMELENGTH] ;
- struct timeval start_time, end_time, curr_time, last_time, forecast_time;
+ /* struct timeval start_time, end_time, curr_time, last_time, forecast_time; */
- struct timezone tzp;
+ /* struct timezone tzp; */
- extern int gettimeofday();
+ /* extern int gettimeofday(); */
- struct tm tmg, tm, tmf, *gmtime(), *localtime();
+ struct tm tml, *gmtime(), *localtime();
- long time_value;
- extern long time();
+ extern time_t time();
+ struct tm start_time, end_time, curr_time, last_time, forecast_time;
+ time_t  rstart_time, rend_time, rcurr_time, rlast_time, rforecast_time; /* raw time */
+ struct tm tm;
  char strcurr[80], strfor[80];
  char *endptr;
- Line 512  static double maxarg1,maxarg2;
+ Line 698  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;}
- Line 535  double dateintmean=0;
+ Line 726  double dateintmean=0;
  double *weight;
  int **s; /* Status */
- double *agedc, **covar, idx;
+ double *agedc;
- int **nbcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff;
+ double  **covar; /**< covar[j,i], value of jth covariate for individual i,
+                   * covar=matrix(0,NCOVMAX,1,n);
+                   * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; */
+ double  idx;
+ int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
+ int *Ndum; /** Freq of modality (tricode */
+ int **codtab; /**< codtab=imatrix(1,100,1,10); */
+ int **Tvard, *Tprod, cptcovprod, *Tvaraff;
  double *lsurv, *lpop, *tpop;
- double ftol=FTOL; /* Tolerance for computing Max Likelihood */
+ double ftol=FTOL; /**< Tolerance for computing Max Likelihood */
- double ftolhess; /* Tolerance for computing hessian */
+ double ftolhess; /**< Tolerance for computing hessian */
  /**************** split *************************/
  static  int split( char *path, char *dirc, char *name, char *ext, char *finame )
- Line 610  void replace_back_to_slash(char *s, char
+ Line 808  void replace_back_to_slash(char *s, char
  }
  char *trimbb(char *out, char *in)
- { /* Trim multiple blanks in line */
+ { /* Trim multiple blanks in line but keeps first blanks if line starts with blanks */
    char *s;
    s=out;
    while (*in != '\0'){
-     while( *in == ' ' && *(in+1) == ' ' && *(in+1) != '\0'){
+     while( *in == ' ' && *(in+1) == ' '){ /* && *(in+1) != '\0'){*/
        in++;
      }
      *out++ = *in++;
- Line 623  char *trimbb(char *out, char *in)
+ Line 821  char *trimbb(char *out, char *in)
    return s;
  }
+ char *cutl(char *blocc, char *alocc, char *in, char occ)
+ {
+   /* cuts string in into blocc and alocc where blocc ends before first occurence of char 'occ'
+      and alocc starts after first occurence of char 'occ' : ex cutv(blocc,alocc,"abcdef2ghi2j",'2')
+      gives blocc="abcdef2ghi" and alocc="j".
+      If occ is not found blocc is null and alocc is equal to in. Returns blocc
+   */
+   char *s, *t;
+   t=in;s=in;
+   while ((*in != occ) && (*in != '\0')){
+     *alocc++ = *in++;
+   }
+   if( *in == occ){
+     *(alocc)='\0';
+     s=++in;
+   }
+   if (s == t) {/* occ not found */
+     *(alocc-(in-s))='\0';
+     in=s;
+   }
+   while ( *in != '\0'){
+     *blocc++ = *in++;
+   }
+   *blocc='\0';
+   return t;
+ }
+ char *cutv(char *blocc, char *alocc, char *in, char occ)
+ {
+   /* cuts string in into blocc and alocc where blocc ends before last occurence of char 'occ'
+      and alocc starts after last occurence of char 'occ' : ex cutv(blocc,alocc,"abcdef2ghi2j",'2')
+      gives blocc="abcdef2ghi" and alocc="j".
+      If occ is not found blocc is null and alocc is equal to in. Returns alocc
+   */
+   char *s, *t;
+   t=in;s=in;
+   while (*in != '\0'){
+     while( *in == occ){
+       *blocc++ = *in++;
+       s=in;
+     }
+     *blocc++ = *in++;
+   }
+   if (s == t) /* occ not found */
+     *(blocc-(in-s))='\0';
+   else
+     *(blocc-(in-s)-1)='\0';
+   in=s;
+   while ( *in != '\0'){
+     *alocc++ = *in++;
+   }
+   *alocc='\0';
+   return s;
+ }
  int nbocc(char *s, char occ)
  {
    int i,j=0;
- Line 635  int nbocc(char *s, char occ)
+ Line 890  int nbocc(char *s, char occ)
    return j;
  }
- void cutv(char *u,char *v, char*t, char occ)
+ /* void cutv(char *u,char *v, char*t, char occ) */
- {
+ /* { */
-   /* cuts string t into u and v where u ends before first occurence of char 'occ'
+ /*   /\* cuts string t into u and v where u ends before last occurence of char 'occ'  */
-      and v starts after first occurence of char 'occ' : ex cutv(u,v,"abcdef2ghi2j",'2')
+ /*      and v starts after last occurence of char 'occ' : ex cutv(u,v,"abcdef2ghi2j",'2') */
-      gives u="abcedf" and v="ghi2j" */
+ /*      gives u="abcdef2ghi" and v="j" *\/ */
-   int i,lg,j,p=0;
+ /*   int i,lg,j,p=0; */
-   i=0;
+ /*   i=0; */
-   for(j=0; j<=strlen(t)-1; j++) {
+ /*   lg=strlen(t); */
-     if((t[j]!= occ) && (t[j+1]== occ)) p=j+1;
+ /*   for(j=0; j<=lg-1; j++) { */
-   }
+ /*     if((t[j]!= occ) && (t[j+1]== occ)) p=j+1; */
+ /*   } */
-   lg=strlen(t);
+ /*   for(j=0; j<p; j++) { */
-   for(j=0; j<p; j++) {
+ /*     (u[j] = t[j]); */
-     (u[j] = t[j]);
+ /*   } */
-   }
+ /*      u[p]='\0'; */
-      u[p]='\0';
-    for(j=0; j<= lg; j++) {
+ /*    for(j=0; j<= lg; j++) { */
-     if (j>=(p+1))(v[j-p-1] = t[j]);
+ /*     if (j>=(p+1))(v[j-p-1] = t[j]); */
+ /*   } */
+ /* } */
+ #ifdef _WIN32
+ char * strsep(char **pp, const char *delim)
+ {
+   char *p, *q;
+   if ((p = *pp) == NULL)
+     return 0;
+   if ((q = strpbrk (p, delim)) != NULL)
+   {
+     *pp = q + 1;
+     *q = '\0';
    }
+   else
+     *pp = 0;
+   return p;
  }
+ #endif
  /********************** nrerror ********************/
- Line 764  double **matrix(long nrl, long nrh, long
+ Line 1037  double **matrix(long nrl, long nrh, long
    for (i=nrl+1; i<=nrh; i++) m[i]=m[i-1]+ncol;
    return m;
-   /* print *(*(m+1)+70) or print m[1][70]; print m+1 or print &(m[1])
+   /* print *(*(m+1)+70) or print m[1][70]; print m+1 or print &(m[1]) or &(m[1][0])
+ m[i] = address of ith row of the table. &(m[i]) is its value which is another adress
+ that of m[i][0]. In order to get the value p m[i][0] but it is unitialized.
     */
  }
- Line 854  char *subdirf3(char fileres[], char *pre
+ Line 1129  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;
- Line 877  double brent(double ax, double bx, doubl
+ Line 1165  double brent(double ax, double bx, doubl
  {
    int iter;
    double a,b,d,etemp;
-   double fu,fv,fw,fx;
+   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;
- Line 892  double brent(double ax, double bx, doubl
+ Line 1180  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)))) { */
- Line 962  void mnbrak(double *ax, double *bx, doub
+ Line 1250  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));
+       (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);
+     ulim=(*bx)+GLIMIT*(*cx-*bx); /* Maximum abscisse where function can be evaluated */
-     if ((*bx-u)*(u-*cx) > 0.0) {
+     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 {
- Line 989  void mnbrak(double *ax, double *bx, doub
+ Line 1285  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 []);
- Line 1015  void linmin(double p[], double xi[], int
+ Line 1315  void linmin(double p[], double xi[], int
    }
    ax=0.0;
    xx=1.0;
-   mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim);
+   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);
+   *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);
- Line 1029  void linmin(double p[], double xi[], int
+ Line 1329  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,"%d day(s) %d hour(s) %d minute(s) %d 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 []))
  {
- Line 1060  void powell(double p[], double **xi, int
+ Line 1356  void powell(double p[], double **xi, int
    xits=vector(1,n);
    *fret=(*func)(p);
    for (j=1;j<=n;j++) pt[j]=p[j];
+     rcurr_time = time(NULL);
    for (*iter=1;;++(*iter)) {
      fp=(*fret);
      ibig=0;
      del=0.0;
-     last_time=curr_time;
+     rlast_time=rcurr_time;
-     (void) gettimeofday(&curr_time,&tzp);
+     /* (void) gettimeofday(&curr_time,&tzp); */
-     printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec);fflush(stdout);
+     rcurr_time = time(NULL);
-     fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec); fflush(ficlog);
+     curr_time = *localtime(&rcurr_time);
- /*     fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tv_sec-start_time.tv_sec); */
+     printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout);
+     fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog);
+ /*     fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */
     for (i=1;i<=n;i++) {
        printf(" %d %.12f",i, p[i]);
        fprintf(ficlog," %d %.12lf",i, p[i]);
- Line 1078  void powell(double p[], double **xi, int
+ Line 1377  void powell(double p[], double **xi, int
      fprintf(ficlog,"\n");
      fprintf(ficrespow,"\n");fflush(ficrespow);
      if(*iter <=3){
-       tm = *localtime(&curr_time.tv_sec);
+       tml = *localtime(&rcurr_time);
-       strcpy(strcurr,asctime(&tm));
+       strcpy(strcurr,asctime(&tml));
- /*       asctime_r(&tm,strcurr); */
+       rforecast_time=rcurr_time;
-       forecast_time=curr_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,curr_time.tv_sec-last_time.tv_sec);
+       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,curr_time.tv_sec-last_time.tv_sec);
+       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){
-         forecast_time.tv_sec=curr_time.tv_sec+(niterf-*iter)*(curr_time.tv_sec-last_time.tv_sec);
+         rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time);
-         tmf = *localtime(&forecast_time.tv_sec);
+         forecast_time = *localtime(&rforecast_time);
- /*      asctime_r(&tmf,strfor); */
+         strcpy(strfor,asctime(&forecast_time));
-         strcpy(strfor,asctime(&tmf));
          itmp = strlen(strfor);
          if(strfor[itmp-1]=='\n')
          strfor[itmp-1]='\0';
-         printf("   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr);
+         printf("   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
-         fprintf(ficlog,"   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr);
+         fprintf(ficlog,"   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
        }
      }
      for (i=1;i<=n;i++) {
        for (j=1;j<=n;j++) xit[j]=xi[j][i];
        fptt=(*fret);
  #ifdef DEBUG
-       printf("fret=%lf \n",*fret);
+           printf("fret=%lf, %lf, %lf \n", *fret, *fret, *fret);
-       fprintf(ficlog,"fret=%lf \n",*fret);
+           fprintf(ficlog, "fret=%lf, %lf, %lf \n", *fret, *fret, *fret);
  #endif
        printf("%d",i);fflush(stdout);
        fprintf(ficlog,"%d",i);fflush(ficlog);
- Line 1122  void powell(double p[], double **xi, int
+ Line 1419  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]);
+         printf(" p(%d)=%.12e",j,p[j]);
-         fprintf(ficlog," p=%.12e",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;
- Line 1161  void powell(double p[], double **xi, int
+ Line 1458  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) {
+     if (fptt < fp) { /* If extrapolated point is better, decide if we keep that new direction or not */
-       t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt);
+       /* (x1 f1=fp), (x2 f2=*fret), (x3 f3=fptt), (xm fm) */
-       if (t < 0.0) {
+       /* From x1 (P0) distance of x2 is at h and x3 is 2h */
-         linmin(p,xit,n,fret,func);
+       /* 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][ibig]=xi[j][n]; /* Replace the direction with biggest decrease by n */
-           xi[j][n]=xit[j];
+           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);
- Line 1185  void powell(double p[], double **xi, int
+ Line 1507  void powell(double p[], double **xi, int
          printf("\n");
          fprintf(ficlog,"\n");
  #endif
-       }
+       } /* end of t negative */
-     }
+     } /* end if (fptt < fp)  */
    }
  }
- Line 1199  double **prevalim(double **prlim, int nl
+ Line 1521  double **prevalim(double **prlim, int nl
    int i, ii,j,k;
    double min, max, maxmin, maxmax,sumnew=0.;
-   double **matprod2();
+   /* double **matprod2(); */ /* test */
    double **out, cov[NCOVMAX+1], **pmij();
    double **newm;
    double agefin, delaymax=50 ; /* Max number of years to converge */
- Line 1215  double **prevalim(double **prlim, int nl
+ Line 1537  double **prevalim(double **prlim, int nl
    for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){
      newm=savm;
      /* Covariates have to be included here again */
-      cov[2]=agefin;
+     cov[2]=agefin;
-       for (k=1; k<=cptcovn;k++) {
+     for (k=1; k<=cptcovn;k++) {
-         cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
+       cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
-         /*      printf("ij=%d k=%d Tvar[k]=%d nbcode=%d cov=%lf codtab[ij][Tvar[k]]=%d \n",ij,k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], codtab[ij][Tvar[k]]);*/
+       /*printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtab[%d][Tvar[%d]]=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], ij, k, codtab[ij][Tvar[k]]);*/
-       }
+     }
-       for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
+     /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
-       for (k=1; k<=cptcovprod;k++)
+     /* for (k=1; k<=cptcovprod;k++) /\* Useless *\/ */
-         cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
+     /*   cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]] * nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; */
-       /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
+     /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
-       /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
+     /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
-       /*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/
+     /*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/
-     out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);
+     /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
+     /* out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /\* Bug Valgrind *\/ */
+     out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /* Bug Valgrind */
      savm=oldm;
      oldm=newm;
      maxmax=0.;
- Line 1240  double **prevalim(double **prlim, int nl
+ Line 1564  double **prevalim(double **prlim, int nl
          sumnew=0;
          for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k];
          prlim[i][j]= newm[i][j]/(1-sumnew);
+         /*printf(" prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d \n", i, j, i, j, prlim[i][j],(int)agefin);*/
          max=FMAX(max,prlim[i][j]);
          min=FMIN(min,prlim[i][j]);
        }
- Line 1256  double **prevalim(double **prlim, int nl
+ Line 1581  double **prevalim(double **prlim, int nl
  double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
  {
-   double s1, s2;
+   /* According to parameters values stored in x and the covariate's values stored in cov,
+      computes the probability to be observed in state j being in state i by appying the
+      model to the ncovmodel covariates (including constant and age).
+      lnpijopii=ln(pij/pii)= aij+bij*age+cij*v1+dij*v2+... = sum_nc=1^ncovmodel xij(nc)*cov[nc]
+      and, according on how parameters are entered, the position of the coefficient xij(nc) of the
+      ncth covariate in the global vector x is given by the formula:
+      j<i nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel
+      j>=i nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel
+      Computes ln(pij/pii) (lnpijopii), deduces pij/pii by exponentiation,
+      sums on j different of i to get 1-pii/pii, deduces pii, and then all pij.
+      Outputs ps[i][j] the probability to be observed in j being in j according to
+      the values of the covariates cov[nc] and corresponding parameter values x[nc+shiftij]
+   */
+   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++){
-         for (nc=1, s2=0.;nc <=ncovmodel; nc++){
+         for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
-           /*s2 += param[i][j][nc]*cov[nc];*/
+           /*lnpijopii += param[i][j][nc]*cov[nc];*/
-           s2 += x[(i-1)*nlstate*ncovmodel+(j-1)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
+           lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];
- /*       printf("Int j<i s1=%.17e, s2=%.17e\n",s1,s2); */
+ /*       printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
          }
-         ps[i][j]=s2;
+         ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
- /*      printf("s1=%.17e, s2=%.17e\n",s1,s2); */
+ /*      printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
        }
        for(j=i+1; j<=nlstate+ndeath;j++){
-         for (nc=1, s2=0.;nc <=ncovmodel; nc++){
+         for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
-           s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
+           /*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/
- /*        printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2); */
+           lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];
+ /*        printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */
          }
-         ps[i][j]=s2;
+         ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
        }
      }
-     /*ps[3][2]=1;*/
      for(i=1; i<= nlstate; i++){
        s1=0;
        for(j=1; j<i; j++){
-         s1+=exp(ps[i][j]);
+         s1+=exp(ps[i][j]); /* In fact sums pij/pii */
          /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
        }
        for(j=i+1; j<=nlstate+ndeath; j++){
-         s1+=exp(ps[i][j]);
+         s1+=exp(ps[i][j]); /* In fact sums pij/pii */
          /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
        }
+       /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
        ps[i][i]=1./(s1+1.);
+       /* Computing other pijs */
        for(j=1; j<i; j++)
          ps[i][j]= exp(ps[i][j])*ps[i][i];
        for(j=i+1; j<=nlstate+ndeath; j++)
- Line 1305  double **pmij(double **ps, double *cov,
+ Line 1645  double **pmij(double **ps, double *cov,
        }
      }
- /*        for(ii=1; ii<= nlstate+ndeath; ii++){ */
+     /* for(ii=1; ii<= nlstate+ndeath; ii++){ */
- /*       for(jj=1; jj<= nlstate+ndeath; jj++){ */
+     /*   for(jj=1; jj<= nlstate+ndeath; jj++){ */
- /*         printf("ddd %lf ",ps[ii][jj]); */
+     /*  printf(" pmij  ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */
- /*       } */
+     /*   } */
- /*       printf("\n "); */
+     /*   printf("\n "); */
- /*        } */
+     /* } */
- /*        printf("\n ");printf("%lf ",cov[2]); */
+     /* printf("\n ");printf("%lf ",cov[2]);*/
-        /*
+     /*
        for(i=1; i<= npar; i++) printf("%f ",x[i]);
        goto end;*/
      return ps;
- Line 1321  double **pmij(double **ps, double *cov,
+ Line 1661  double **pmij(double **ps, double *cov,
  /**************** Product of 2 matrices ******************/
- double **matprod2(double **out, double **in,long nrl, long nrh, long ncl, long nch, long ncolol, long ncoloh, double **b)
+ double **matprod2(double **out, double **in,int nrl, int nrh, int ncl, int nch, int ncolol, int ncoloh, double **b)
  {
    /* Computes the matrix product of in(1,nrh-nrl+1)(1,nch-ncl+1) times
       b(1,nch-ncl+1)(1,ncoloh-ncolol+1) into out(...) */
    /* in, b, out are matrice of pointers which should have been initialized
       before: only the contents of out is modified. The function returns
       a pointer to pointers identical to out */
-   long i, j, k;
+   int i, j, k;
    for(i=nrl; i<= nrh; i++)
-     for(k=ncolol; k<=ncoloh; k++)
+     for(k=ncolol; k<=ncoloh; k++){
-       for(j=ncl,out[i][k]=0.; j<=nch; j++)
+       out[i][k]=0.;
-         out[i][k] +=in[i][j]*b[j][k];
+       for(j=ncl; j<=nch; j++)
+         out[i][k] +=in[i][j]*b[j][k];
+     }
    return out;
  }
- Line 1375  double ***hpxij(double ***po, int nhstep
+ Line 1716  double ***hpxij(double ***po, int nhstep
          cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
        for (k=1; k<=cptcovage;k++)
          cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
-       for (k=1; k<=cptcovprod;k++)
+       for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */
          cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
- Line 1397  double ***hpxij(double ***po, int nhstep
+ Line 1738  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)
- Line 1415  double func( double *x)
+ Line 1775  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.;
    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];
+       /* Computes the values of the ncovmodel covariates of the model
+          depending if the covariates are fixed or variying (age dependent) and stores them in cov[]
+          Then computes with function pmij which return a matrix p[i][j] giving the elementary probability
+          to be observed in j being in i according to the model.
+        */
+       for (k=1; k<=cptcovn;k++){ /* Simple and product covariates without age* products */
+         cov[2+k]=covar[Tvar[k]][i];
+       }
+       /* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4]
+          is 6, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2]
+          has been calculated etc */
        for(mi=1; mi<= wav[i]-1; mi++){
          for (ii=1;ii<=nlstate+ndeath;ii++)
            for (j=1;j<=nlstate+ndeath;j++){
- Line 1432  double func( double *x)
+ Line 1805  double func( double *x)
            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];
+             cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; /* Tage[kk] gives the data-covariate associated with age */
            }
            out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
 ,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
- Line 1686  double funcone( double *x)
+ Line 2059  double funcone( double *x)
          for (kk=1; kk<=cptcovage;kk++) {
            cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
          }
+         /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
          out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
 ,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
+         /* 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 */
- Line 1712  double funcone( double *x)
+ Line 2088  double funcone( double *x)
          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])); /* exponential inter-extrapolation */
        } else if (mle==4){  /* mle=4 no inter-extrapolation */
          lli=log(out[s1][s2]); /* Original formula */
-       } else{  /* ml>=5 no inter-extrapolation no jackson =0.8a */
+       } else{  /* mle=0 back to 1 */
-         lli=log(out[s1][s2]); /* Original formula */
+         lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
+         /*lli=log(out[s1][s2]); */ /* Original formula */
        } /* End of if */
        ipmx +=1;
        sw += weight[i];
        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,"%9d %6d %2d %2d %1d %1d %3d %11.6f %8.4f\
+         fprintf(ficresilk,"%9ld %6d %2d %2d %1d %1d %3d %11.6f %8.4f\
   %11.6f %11.6f %11.6f ", \
                  num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],
 *weight[i]*lli,out[s1][s2],savm[s1][s2]);
- Line 1782  void likelione(FILE *ficres,double p[],
+ Line 2159  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++)
- Line 1803  void mlikeli(FILE *ficres,double p[], in
+ Line 2192  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));
+   printf("\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
-   fprintf(ficlog,"\n#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
+   fprintf(ficlog,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
-   fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
+   fprintf(ficres,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
  }
- Line 1819  void hesscov(double **matcov, double p[]
+ Line 2235  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);
- Line 1932  double hessii(double x[], double delta,
+ Line 2348  double hessii(double x[], double delta,
    fx=func(x);
    for (i=1;i<=npar;i++) p2[i]=x[i];
-   for(l=0 ; l <=lmax; l++){
+   for(l=0 ; l <=lmax; l++){  /* Enlarging the zone around the Maximum */
      l1=pow(10,l);
      delts=delt;
      for(k=1 ; k <kmax; k=k+1){
        delt = delta*(l1*k);
        p2[theta]=x[theta] +delt;
-       k1=func(p2)-fx;
+       k1=func(p2)-fx;   /* Might be negative if too close to the theoretical maximum */
        p2[theta]=x[theta]-delt;
        k2=func(p2)-fx;
        /*res= (k1-2.0*fx+k2)/delt/delt; */
- Line 1953  double hessii(double x[], double delta,
+ Line 2369  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;
- Line 1968  double hessii(double x[], double delta,
+ Line 2384  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;
- Line 2083  void pstamp(FILE *fichier)
+ Line 2499  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;
- Line 2107  void  freqsummary(char fileres[], int ia
+ Line 2523  void  freqsummary(char fileres[], int ia
    first=1;
-   for(k1=1; k1<=j;k1++){
+   /* for(k1=1; k1<=j ; k1++){   /* Loop on covariates */
-     for(i1=1; i1<=ncodemax[k1];i1++){
+   /*  for(i1=1; i1<=ncodemax[k1];i1++){ /* Now it is 2 */
-       j1++;
+   /*    j1++;
+ */
+   for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){
        /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
          scanf("%d", i);*/
        for (i=-5; i<=nlstate+ndeath; i++)
          for (jk=-5; jk<=nlstate+ndeath; jk++)
            for(m=iagemin; m <= iagemax+3; m++)
              freq[i][jk][m]=0;
-     for (i=1; i<=nlstate; i++)
+       for (i=1; i<=nlstate; i++)
-       for(m=iagemin; m <= iagemax+3; m++)
+         for(m=iagemin; m <= iagemax+3; m++)
-         prop[i][m]=0;
+           prop[i][m]=0;
        dateintsum=0;
        k2cpt=0;
        for (i=1; i<=imx; i++) {
          bool=1;
-         if  (cptcovn>0) {
+         if  (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
            for (z1=1; z1<=cptcoveff; z1++)
-             if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]])
+             if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]){
-               bool=0;
+                 /* Tests if the value of each of the covariates of i is equal to filter j1 */
+               bool=0;
+               /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtab[%d][%d]=%d, nbcode[Tvaraff][codtab[%d][%d]=%d, j1=%d\n",
+                 bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtab[j1][z1],
+                 j1,z1,nbcode[Tvaraff[z1]][codtab[j1][z1]],j1);*/
+               /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtab[7][3]=1 and nbcde[3][?]=1*/
+             }
          }
          if (bool==1){
            for(m=firstpass; m<=lastpass; m++){
              k2=anint[m][i]+(mint[m][i]/12.);
- Line 2149  void  freqsummary(char fileres[], int ia
+ Line 2574  void  freqsummary(char fileres[], int ia
                /*}*/
            }
          }
-       }
+       } /* end i */
        /*      fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
        pstamp(ficresp);
- Line 2157  void  freqsummary(char fileres[], int ia
+ Line 2582  void  freqsummary(char fileres[], int ia
          fprintf(ficresp, "\n#********** Variable ");
          for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
          fprintf(ficresp, "**********\n#");
+         fprintf(ficlog, "\n#********** Variable ");
+         for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
+         fprintf(ficlog, "**********\n#");
        }
        for(i=1; i<=nlstate;i++)
          fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);
- Line 2233  void  freqsummary(char fileres[], int ia
+ Line 2661  void  freqsummary(char fileres[], int ia
            printf("Others in log...\n");
          fprintf(ficlog,"\n");
        }
-     }
+       /*}*/
    }
    dateintmean=dateintsum/k2cpt;
- Line 2252  void prevalence(double ***probs, double
+ Line 2680  void prevalence(double ***probs, double
       We still use firstpass and lastpass as another selection.
    */
-   int i, m, jk, k1, i1, j1, bool, z1,j;
+   int i, m, jk, j1, bool, z1,j;
-   double ***freq; /* Frequencies */
-   double *pp, **prop;
+   double **prop;
-   double pos,posprop;
+   double posprop;
    double  y2; /* in fractional years */
    int iagemin, iagemax;
+   int first; /** to stop verbosity which is redirected to log file */
    iagemin= (int) agemin;
    iagemax= (int) agemax;
- Line 2266  void prevalence(double ***probs, double
+ Line 2695  void prevalence(double ***probs, double
    /*  freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/
    j1=0;
-   j=cptcoveff;
+   /*j=cptcoveff;*/
    if (cptcovn<1) {j=1;ncodemax[1]=1;}
-   for(k1=1; k1<=j;k1++){
+   first=1;
-     for(i1=1; i1<=ncodemax[k1];i1++){
+   for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){
-       j1++;
+     /*for(i1=1; i1<=ncodemax[k1];i1++){
+       j1++;*/
        for (i=1; i<=nlstate; i++)
          for(m=iagemin; m <= iagemax+3; m++)
- Line 2301  void prevalence(double ***probs, double
+ Line 2731  void prevalence(double ***probs, double
          }
        }
        for(i=iagemin; i <= iagemax+3; i++){
          for(jk=1,posprop=0; jk <=nlstate ; jk++) {
            posprop += prop[jk][i];
          }
          for(jk=1; jk <=nlstate ; jk++){
            if( i <=  iagemax){
              if(posprop>=1.e-5){
                probs[i][jk][j1]= prop[jk][i]/posprop;
-             } else
+             } else{
-               printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\n",jk,i,j1,probs[i][jk][j1]);
+               if(first==1){
+                 first=0;
+                 printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]);
+               }
+             }
            }
          }/* end jk */
        }/* end i */
-     } /* end i1 */
+     /*} *//* end i1 */
-   } /* end k1 */
+   } /* end j1 */
    /*  free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);*/
    /*free_vector(pp,1,nlstate);*/
- Line 2341  void  concatwav(int wav[], int **dh, int
+ Line 2774  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++){
- Line 2436  void  concatwav(int wav[], int **dh, int
+ Line 2869  void  concatwav(int wav[], int **dh, int
              dh[mi][i]=jk;
              bh[mi][i]=0;
            }else{ /* We want a negative bias in order to only have interpolation ie
-                   * at the price of an extra matrix product in likelihood */
+                   * to avoid the price of an extra matrix product in likelihood */
              dh[mi][i]=jk+1;
              bh[mi][i]=ju;
            }
- Line 2462  void  concatwav(int wav[], int **dh, int
+ Line 2895  void  concatwav(int wav[], int **dh, int
    }
    jmean=sum/k;
    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);
+   fprintf(ficlog,"Delay (in months) between two waves Min=%d (for indiviudal %d) Max=%d (%d) Mean=%f\n\n ",jmin, ijmin, jmax, ijmax, jmean);
   }
  /*********** Tricode ****************************/
- void tricode(int *Tvar, int **nbcode, int imx)
+ void tricode(int *Tvar, int **nbcode, int imx, int *Ndum)
  {
+   /**< Uses cptcovn+2*cptcovprod as the number of covariates */
-   /*      Tvar[i]=atoi(stre); /* find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 */
+   /*      Tvar[i]=atoi(stre);  find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1
+   /* Boring subroutine which should only output nbcode[Tvar[j]][k]
+    * Tvar[5] in V2+V1+V3*age+V2*V4 is 2 (V2)
+   /* nbcode[Tvar[j]][1]=
+   */
+   int ij=1, k=0, j=0, i=0, maxncov=NCOVMAX;
+   int modmaxcovj=0; /* Modality max of covariates j */
+   int cptcode=0; /* Modality max of covariates j */
+   int modmincovj=0; /* Modality min of covariates j */
-   int Ndum[20],ij=1, k=0, j=0, i=0, maxncov=NCOVMAX;
-   int cptcode=0;
    cptcoveff=0;
-   for (k=0; k<maxncov; k++) Ndum[k]=0;
+   for (k=-1; k < maxncov; k++) Ndum[k]=0;
-   for (k=1; k<=7; k++) ncodemax[k]=0; /* Horrible constant again */
+   for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */
-   for (j=1; j<=(cptcovn+2*cptcovprod); j++) { /* For each covariate */
+   /* Loop on covariates without age and products */
-     for (i=1; i<=imx; i++) { /*reads the data file to get the maximum
+   for (j=1; j<=(cptcovs); j++) { /* model V1 + V2*age+ V3 + V3*V4 : V1 + V3 = 2 only */
-                                modality*/
+     for (i=1; i<=imx; i++) { /* Lopp on individuals: reads the data file to get the maximum value of the
-       ij=(int)(covar[Tvar[j]][i]); /* ij is the modality of this individual, might be -1*/
+                                modality of this covariate Vj*/
-       Ndum[ij]++; /*counts the occurence of this modality */
+       ij=(int)(covar[Tvar[j]][i]); /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i
+                                     * If product of Vn*Vm, still boolean *:
+                                     * If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables
+                                     * 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0   */
+       /* Finds for covariate j, n=Tvar[j] of Vn . ij is the
+                                       modality of the nth covariate of individual i. */
+       if (ij > modmaxcovj)
+         modmaxcovj=ij;
+       else if (ij < modmincovj)
+         modmincovj=ij;
+       if ((ij < -1) && (ij > NCOVMAX)){
+         printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX );
+         exit(1);
+       }else
+       Ndum[ij]++; /*counts and stores the occurence of this modality 0, 1, -1*/
+       /*  If coded 1, 2, 3 , counts the number of 1 Ndum[1], number of 2, Ndum[2], etc */
        /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/
-       if (ij > cptcode) cptcode=ij; /* getting the maximum value of the modality of the covariate  (should be 0 or 1 now)
+       /* getting the maximum value of the modality of the covariate
-                                        Tvar[j]. If V=sex and male is 0 and
+          (should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and
-                                        female is 1, then  cptcode=1.*/
+          female is 1, then modmaxcovj=1.*/
      }
+     printf(" Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", j, Tvar[j], modmincovj, modmaxcovj);
-     for (i=0; i<=cptcode; i++) { /* i=-1 ?*/
+     cptcode=modmaxcovj;
-       if(Ndum[i]!=0) ncodemax[j]++; /* Nomber of modalities of the j
+     /* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */
-                                        th covariate. In fact
+    /*for (i=0; i<=cptcode; i++) {*/
-                                        ncodemax[j]=2
+     for (i=modmincovj;  i<=modmaxcovj; i++) { /* i=-1 ? 0 and 1*//* For each value of the modality of model-cov j */
-                                        (dichotom. variables only) but
+       printf("Frequencies of covariates %d V%d %d\n", j, Tvar[j], Ndum[i]);
-                                        it can be more */
+       if( Ndum[i] != 0 ){ /* Counts if nobody answered, empty modality */
+         ncodemax[j]++;  /* ncodemax[j]= Number of non-null modalities of the j th covariate. */
+       }
+       /* In fact  ncodemax[j]=2 (dichotom. variables only) but it could be more for
+          historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */
      } /* Ndum[-1] number of undefined modalities */
-     ij=1;
+     /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */
-     for (i=1; i<=ncodemax[j]; i++) { /* i= 1 to 2 */
+     /* For covariate j, modalities could be 1, 2, 3, 4. If Ndum[2]=0 ncodemax[j] is not 4 but 3 */
-       for (k=0; k<= maxncov; k++) { /* k=-1 ?*/
+     /* If Ndum[3}= 635; Ndum[4]=0; Ndum[5]=0; Ndum[6]=27; Ndum[7]=125;
+        modmincovj=3; modmaxcovj = 7;
+        There are only 3 modalities non empty (or 2 if 27 is too few) : ncodemax[j]=3;
+        which will be coded 0, 1, 2 which in binary on 3-1 digits are 0=00 1=01, 2=10; defining two dummy
+        variables V1_1 and V1_2.
+        nbcode[Tvar[j]][ij]=k;
+        nbcode[Tvar[j]][1]=0;
+        nbcode[Tvar[j]][2]=1;
+        nbcode[Tvar[j]][3]=2;
+     */
+     ij=1; /* ij is similar to i but can jumps over null modalities */
+     for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 */
+       for (k=0; k<= cptcode; k++) { /* k=-1 ? k=0 to 1 *//* Could be 1 to 4 */
+         /*recode from 0 */
          if (Ndum[k] != 0) { /* If at least one individual responded to this modality k */
            nbcode[Tvar[j]][ij]=k;  /* stores the modality in an array nbcode.
                                       k is a modality. If we have model=V1+V1*sex
- Line 2507  void tricode(int *Tvar, int **nbcode, in
+ Line 2980  void tricode(int *Tvar, int **nbcode, in
            ij++;
          }
          if (ij > ncodemax[j]) break;
-       }
+       }  /* end of loop on */
-     }
+     } /* end of loop on modality */
-   }
+   } /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/
-  for (k=0; k< maxncov; k++) Ndum[k]=0;
+  for (k=-1; k< maxncov; k++) Ndum[k]=0;
-  for (i=1; i<=ncovmodel-2; i++) { /* -2, cste and age */
+   for (i=1; i<=ncovmodel-2; i++) { /* -2, cste and age */
     /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/
     ij=Tvar[i]; /* Tvar might be -1 if status was unknown */
     Ndum[ij]++;
   }
   ij=1;
-  for (i=1; i<= maxncov; i++) {
+  for (i=0; i<=  maxncov-1; i++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */
+    /*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/
     if((Ndum[i]!=0) && (i<=ncovcol)){
-      Tvaraff[ij]=i; /*For printing */
+      /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/
+      Tvaraff[ij]=i; /*For printing (unclear) */
       ij++;
-    }
+    }else
+        Tvaraff[ij]=0;
   }
   ij--;
-  cptcoveff=ij; /*Number of simple covariates*/
+  cptcoveff=ij; /*Number of total covariates*/
  }
  /*********** Health Expectancies ****************/
  void evsij(double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int cij, int estepm,char strstart[] )
  {
    /* 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;
- Line 2858  void varevsij(char optionfilefiname[], d
+ Line 3336  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 */
- Line 3125  void varevsij(char optionfilefiname[], d
+ Line 3603  void varevsij(char optionfilefiname[], d
    free_vector(gmp,nlstate+1,nlstate+ndeath);
    free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath);
    free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/
-   fprintf(ficgp,"\nunset parametric;unset label; set ter png small;set size 0.65, 0.65");
+   fprintf(ficgp,"\nunset parametric;unset label; set ter png small size 320, 240");
    /* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */
    fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";");
  /*   fprintf(ficgp,"\n plot \"%s\"  u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */
  /*   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */
  /*   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */
-   fprintf(ficgp,"\n plot \"%s\"  u 1:($3) not w l 1 ",subdirf(fileresprobmorprev));
+   fprintf(ficgp,"\n plot \"%s\"  u 1:($3) not w l lt 1 ",subdirf(fileresprobmorprev));
-   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)) t \"95\%% interval\" w l 2 ",subdirf(fileresprobmorprev));
+   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)) t \"95\%% interval\" w l lt 2 ",subdirf(fileresprobmorprev));
-   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)) not w l 2 ",subdirf(fileresprobmorprev));
+   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)) not w l lt 2 ",subdirf(fileresprobmorprev));
    fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev));
    fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"%s%s.png\"> <br>\n", estepm,subdirf3(optionfilefiname,"varmuptjgr",digitp),digit);
    /*  fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", stepm,YEARM,digitp,digit);
- Line 3158  void varprevlim(char fileres[], double *
+ Line 3636  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;
- Line 3240  void varprevlim(char fileres[], double *
+ Line 3717  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;
+   int first=1, first1, first2;
    double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp;
    double **dnewm,**doldm;
    double *xp;
    double *gp, *gm;
    double **gradg, **trgradg;
    double **mu;
-   double age,agelim, cov[NCOVMAX];
+   double age, cov[NCOVMAX+1];
    double std=2.0; /* Number of standard deviation wide of confidence ellipsoids */
    int theta;
    char fileresprob[FILENAMELENGTH];
    char fileresprobcov[FILENAMELENGTH];
    char fileresprobcor[FILENAMELENGTH];
    double ***varpij;
    strcpy(fileresprob,"prob");
- Line 3329  standard deviations wide on each axis. <
+ Line 3805  standard deviations wide on each axis. <
  To be simple, these graphs help to understand the significativity of each parameter in relation to a second other one.<br> \n");
    cov[1]=1;
-   tj=cptcoveff;
+   /* tj=cptcoveff; */
+   tj = (int) pow(2,cptcoveff);
    if (cptcovn<1) {tj=1;ncodemax[1]=1;}
    j1=0;
-   for(t=1; t<=tj;t++){
+   for(j1=1; j1<=tj;j1++){
-     for(i1=1; i1<=ncodemax[t];i1++){
+     /*for(i1=1; i1<=ncodemax[t];i1++){ */
-       j1++;
+     /*j1++;*/
        if  (cptcovn>0) {
          fprintf(ficresprob, "\n#********** Variable ");
          for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
- Line 3357  To be simple, these graphs help to under
+ Line 3834  To be simple, these graphs help to under
          fprintf(ficresprobcor, "**********\n#");
        }
+       gradg=matrix(1,npar,1,(nlstate)*(nlstate+ndeath));
+       trgradg=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);
+       gp=vector(1,(nlstate)*(nlstate+ndeath));
+       gm=vector(1,(nlstate)*(nlstate+ndeath));
        for (age=bage; age<=fage; age ++){
          cov[2]=age;
          for (k=1; k<=cptcovn;k++) {
-           cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];
+           cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];/* j1 1 2 3 4
+                                                          * 1  1 1 1 1
+                                                          * 2  2 1 1 1
+                                                          * 3  1 2 1 1
+                                                          */
+           /* nbcode[1][1]=0 nbcode[1][2]=1;*/
          }
          for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
          for (k=1; k<=cptcovprod;k++)
            cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
-         gradg=matrix(1,npar,1,(nlstate)*(nlstate+ndeath));
-         trgradg=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);
-         gp=vector(1,(nlstate)*(nlstate+ndeath));
-         gm=vector(1,(nlstate)*(nlstate+ndeath));
          for(theta=1; theta <=npar; theta++){
            for(i=1; i<=npar; i++)
- Line 3407  To be simple, these graphs help to under
+ Line 3889  To be simple, these graphs help to under
          matprod2(dnewm,trgradg,1,(nlstate)*(nlstate+ndeath),1,npar,1,npar,matcov);
          matprod2(doldm,dnewm,1,(nlstate)*(nlstate+ndeath),1,npar,1,(nlstate)*(nlstate+ndeath),gradg);
-         free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
-         free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));
-         free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
-         free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
          pmij(pmmij,cov,ncovmodel,x,nlstate);
- Line 3444  To be simple, these graphs help to under
+ Line 3922  To be simple, these graphs help to under
          i=0;
          for (k=1; k<=(nlstate);k++){
            for (l=1; l<=(nlstate+ndeath);l++){
-             i=i++;
+             i++;
              fprintf(ficresprobcov,"\n%d %d-%d",(int)age,k,l);
              fprintf(ficresprobcor,"\n%d %d-%d",(int)age,k,l);
              for (j=1; j<=i;j++){
+               /* printf(" k=%d l=%d i=%d j=%d\n",k,l,i,j);fflush(stdout); */
                fprintf(ficresprobcov," %11.3e",varpij[i][j][(int)age]);
                fprintf(ficresprobcor," %11.3e",varpij[i][j][(int) age]/sqrt(varpij[i][i][(int) age])/sqrt(varpij[j][j][(int)age]));
              }
            }
          }/* end of loop for state */
        } /* end of loop for age */
+       free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
+       free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));
+       free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
+       free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
        /* Confidence intervalle of pij  */
        /*
          fprintf(ficgp,"\nunset parametric;unset label");
- Line 3467  To be simple, these graphs help to under
+ Line 3950  To be simple, these graphs help to under
        */
        /* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/
-       first1=1;
+       first1=1;first2=2;
        for (k2=1; k2<=(nlstate);k2++){
          for (l2=1; l2<=(nlstate+ndeath);l2++){
            if(l2==k2) continue;
- Line 3489  To be simple, these graphs help to under
+ Line 3972  To be simple, these graphs help to under
                    lc1=((v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.;
                    lc2=((v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.;
                    if ((lc2 <0) || (lc1 <0) ){
-                     printf("Error: One eigen value of 2x2 matrix of covariance is negative, lc1=%11.3e, lc2=%11.3e, v1=%11.3e, v2=%11.3e, cv12=%11.3e.\n It means that the matrix was not well estimated (varpij), for i=%2d, j=%2d, age=%4d .\n See files %s and %s. Continuing by making them positive: WRONG RESULTS.\n", lc1, lc2, v1, v2, cv12, i, j, (int)age,fileresprobcov, fileresprobcor);
+                     if(first2==1){
-                     fprintf(ficlog,"Error: One eigen value of 2x2 matrix of covariance is negative, lc1=%11.3e, lc2=%11.3e, v1=%11.3e, v2=%11.3e, cv12=%11.3e\n", lc1, lc2, v1, v2, cv12);fflush(ficlog);
+                       first1=0;
-                     lc1=fabs(lc1);
+                     printf("Strange: j1=%d One eigen value of 2x2 matrix of covariance is negative, lc1=%11.3e, lc2=%11.3e, v1=%11.3e, v2=%11.3e, cv12=%11.3e.\n It means that the matrix was not well estimated (varpij), for i=%2d, j=%2d, age=%4d .\n See files %s and %s. Probably WRONG RESULTS. See log file for details...\n", j1, lc1, lc2, v1, v2, cv12, i, j, (int)age,fileresprobcov, fileresprobcor);
-                     lc2=fabs(lc2);
+                     }
+                     fprintf(ficlog,"Strange: j1=%d One eigen value of 2x2 matrix of covariance is negative, lc1=%11.3e, lc2=%11.3e, v1=%11.3e, v2=%11.3e, cv12=%11.3e.\n It means that the matrix was not well estimated (varpij), for i=%2d, j=%2d, age=%4d .\n See files %s and %s. Probably WRONG RESULTS.\n", j1, lc1, lc2, v1, v2, cv12, i, j, (int)age,fileresprobcov, fileresprobcor);fflush(ficlog);
+                     /* lc1=fabs(lc1); */ /* If we want to have them positive */
+                     /* lc2=fabs(lc2); */
                    }
                    /* Eigen vectors */
- Line 3514  To be simple, these graphs help to under
+ Line 4000  To be simple, these graphs help to under
                      first=0;
                      fprintf(ficgp,"\nset parametric;unset label");
                      fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2);
-                     fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");
+                     fprintf(ficgp,"\nset ter png small size 320, 240");
                      fprintf(fichtmcov,"\n<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup>\
   :<a href=\"%s%d%1d%1d-%1d%1d.png\">\
  %s%d%1d%1d-%1d%1d.png</A>, ",k1,l1,k2,l2,\
- Line 3545  To be simple, these graphs help to under
+ Line 4031  To be simple, these graphs help to under
            } /* k12 */
          } /*l1 */
        }/* k1 */
-     } /* loop covariates */
+       /* } /* loop covariates */
    }
    free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);
    free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);
- Line 3591  void printinghtml(char fileres[], char t
+ Line 4077  void printinghtml(char fileres[], char t
  fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");
-  m=cptcoveff;
+  m=pow(2,cptcoveff);
   if (cptcovn < 1) {m=1;ncodemax[1]=1;}
   jj1=0;
- Line 3605  fprintf(fichtm," \n<ul><li><b>Graphs</b>
+ Line 4091  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: <a href=\"%s%d1.png\">%s%d1.png</a><br> \
+      fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s%d_1.png\">%s%d_1.png</a><br> \
- <img src=\"%s%d1.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);
+ <img src=\"%s%d_1.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),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: <a href=\"%s%d2.png\">%s%d2.png</a><br> \
+  before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: <a href=\"%s%d_2.png\">%s%d_2.png</a><br> \
- <img src=\"%s%d2.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);
+ <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++){
+        for(cpt=1; cpt<=nlstate;cpt++){
-          fprintf(fichtm,"<br>- Period (stable) prevalence in each health state : <a href=\"%s%d%d.png\">%s%d%d.png</a><br> \
+          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\">",subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1);
+ <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 : <a href=\"%s%d%d.png\">%s%d%d.png</a> <br> \
+         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> \
- <img src=\"%s%d%d.png\">",cpt,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1);
+ <img src=\"%s%d%d.png\">",cpt,nlstate,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1);
       }
     } /* end i1 */
   }/* End k1 */
- Line 3665  fprintf(fichtm," \n<ul><li><b>Graphs</b>
+ Line 4151  fprintf(fichtm," \n<ul><li><b>Graphs</b>
   fflush(fichtm);
   fprintf(fichtm," <ul><li><b>Graphs</b></li><p>");
-  m=cptcoveff;
+  m=pow(2,cptcoveff);
   if (cptcovn < 1) {m=1;ncodemax[1]=1;}
   jj1=0;
- Line 3680  fprintf(fichtm," \n<ul><li><b>Graphs</b>
+ Line 4166  fprintf(fichtm," \n<ul><li><b>Graphs</b>
       }
       for(cpt=1; cpt<=nlstate;cpt++) {
         fprintf(fichtm,"<br>- Observed (cross-sectional) and period (incidence based) \
- prevalence (with 95%% confidence interval) in state (%d): %s%d%d.png <br>\
+ prevalence (with 95%% confidence interval) in state (%d): %s%d_%d.png <br>\
- <img src=\"%s%d%d.png\">",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1);
+ <img src=\"%s%d_%d.png\">",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1);
       }
       fprintf(fichtm,"\n<br>- Total life expectancy by age and \
  health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \
- Line 3699  true period expectancies (those weighted
+ Line 4185  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); */
- Line 3714  void printinggnuplot(char fileres[], cha
+ Line 4200  void printinggnuplot(char fileres[], cha
    strcpy(dirfileres,optionfilefiname);
    strcpy(optfileres,"vpl");
   /* 1eme*/
+   fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'vpl' files\n");
    for (cpt=1; cpt<= nlstate ; cpt ++) {
-    for (k1=1; k1<= m ; k1 ++) {
+     for (k1=1; k1<= m ; k1 ++) { /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */
-      fprintf(ficgp,"\nset out \"%s%d%d.png\" \n",subdirf2(optionfilefiname,"v"),cpt,k1);
+      fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"v"),cpt,k1);
-      fprintf(ficgp,"\n#set out \"v%s%d%d.png\" \n",optionfilefiname,cpt,k1);
+      fprintf(ficgp,"\n#set out \"v%s%d_%d.png\" \n",optionfilefiname,cpt,k1);
       fprintf(ficgp,"set xlabel \"Age\" \n\
  set ylabel \"Probability\" \n\
- set ter png small\n\
+ set ter png small size 320, 240\n\
- set size 0.65,0.65\n\
  plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,subdirf2(fileres,"vpl"),k1-1,k1-1);
       for (i=1; i<= nlstate ; i ++) {
         if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
         else        fprintf(ficgp," \%%*lf (\%%*lf)");
       }
-      fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);
+      fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);
       for (i=1; i<= nlstate ; i ++) {
         if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
         else fprintf(ficgp," \%%*lf (\%%*lf)");
       }
-      fprintf(ficgp,"\" t\"95\%% CI\" w l 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);
+      fprintf(ficgp,"\" t\"95\%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);
       for (i=1; i<= nlstate ; i ++) {
         if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
         else fprintf(ficgp," \%%*lf (\%%*lf)");
       }
-      fprintf(ficgp,"\" t\"\" w l 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l 2",subdirf2(fileres,"p"),k1-1,k1-1,2+4*(cpt-1));
+      fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l lt 2",subdirf2(fileres,"p"),k1-1,k1-1,2+4*(cpt-1));
     }
    }
    /*2 eme*/
+   fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files\n");
    for (k1=1; k1<= m ; k1 ++) {
      fprintf(ficgp,"\nset out \"%s%d.png\" \n",subdirf2(optionfilefiname,"e"),k1);
-     fprintf(ficgp,"set ylabel \"Years\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] ",ageminpar,fage);
+     fprintf(ficgp,"set ylabel \"Years\" \nset ter png small size 320, 240\nplot [%.f:%.f] ",ageminpar,fage);
      for (i=1; i<= nlstate+1 ; i ++) {
        k=2*i;
- Line 3761  plot [%.f:%.f] \"%s\" every :::%d::%d u
+ Line 4247  plot [%.f:%.f] \"%s\" every :::%d::%d u
          if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
          else fprintf(ficgp," \%%*lf (\%%*lf)");
        }
-       fprintf(ficgp,"\" t\"\" w l 0,");
+       fprintf(ficgp,"\" t\"\" w l lt 0,");
        fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2+$3*2) \"\%%lf",subdirf2(fileres,"t"),k1-1,k1-1);
        for (j=1; j<= nlstate+1 ; j ++) {
          if (j==i) fprintf(ficgp," \%%lf (\%%lf)");
          else fprintf(ficgp," \%%*lf (\%%*lf)");
        }
-       if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0");
+       if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
-       else fprintf(ficgp,"\" t\"\" w l 0,");
+       else fprintf(ficgp,"\" t\"\" w l lt 0,");
      }
    }
- Line 3779  plot [%.f:%.f] \"%s\" every :::%d::%d u
+ Line 4265  plot [%.f:%.f] \"%s\" every :::%d::%d u
        /*       k=2+nlstate*(2*cpt-2); */
        k=2+(nlstate+1)*(cpt-1);
        fprintf(ficgp,"\nset out \"%s%d%d.png\" \n",subdirf2(optionfilefiname,"exp"),cpt,k1);
-       fprintf(ficgp,"set ter png small\n\
+       fprintf(ficgp,"set ter png small size 320, 240\n\
- set size 0.65,0.65\n\
  plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileres,"e"),k1-1,k1-1,k,cpt);
        /*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
          for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) ");
- Line 3800  plot [%.f:%.f] \"%s\" every :::%d::%d u
+ Line 4285  plot [%.f:%.f] \"%s\" every :::%d::%d u
    }
    /* CV preval stable (period) */
-   for (k1=1; k1<= m ; k1 ++) {
+   for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */
-     for (cpt=1; cpt<=nlstate ; cpt ++) {
+     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
        k=3;
-       fprintf(ficgp,"\nset out \"%s%d%d.png\" \n",subdirf2(optionfilefiname,"p"),cpt,k1);
+       fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, cov=%d state=%d",k1, cpt);
+       fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"p"),cpt,k1);
        fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
- set ter png small\nset size 0.65,0.65\n\
+ set ter png small size 320, 240\n\
  unset log y\n\
- plot [%.f:%.f] \"%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,subdirf2(fileres,"pij"),k1,k+cpt+1,k+1);
+ plot [%.f:%.f]  ", ageminpar, agemaxpar);
+       for (i=1; i<= nlstate ; i ++){
-       for (i=1; i< nlstate ; i ++)
+         if(i==1)
-         fprintf(ficgp,"+$%d",k+i+1);
+           fprintf(ficgp,"\"%s\"",subdirf2(fileres,"pij"));
-       fprintf(ficgp,")) t\"prev(%d,%d)\" w l",cpt,cpt+1);
+         else
+           fprintf(ficgp,", '' ");
-       l=3+(nlstate+ndeath)*cpt;
+         l=(nlstate+ndeath)*(i-1)+1;
-       fprintf(ficgp,",\"%s\" u ($1==%d ? ($3):1/0):($%d/($%d",subdirf2(fileres,"pij"),k1,l+cpt+1,l+1);
+         fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
-       for (i=1; i< nlstate ; i ++) {
+         for (j=1; j<= (nlstate-1) ; j ++)
-         l=3+(nlstate+ndeath)*cpt;
+           fprintf(ficgp,"+$%d",k+l+j);
-         fprintf(ficgp,"+$%d",l+i+1);
+         fprintf(ficgp,")) t \"prev(%d,%d)\" w l",i,cpt);
-       }
+       } /* nlstate */
-       fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);
+       fprintf(ficgp,"\n");
-     }
+     } /* end cpt state*/
-   }
+   } /* end covariate */
    /* proba elementaires */
    for(i=1,jk=1; i <=nlstate; i++){
- Line 3835  plot [%.f:%.f] \"%s\" u ($1==%d ? ($3):1
+ Line 4321  plot [%.f:%.f] \"%s\" u ($1==%d ? ($3):1
        }
      }
     }
+   /*goto avoid;*/
     for(ng=1; ng<=2;ng++){ /* Number of graphics: first is probabilities second is incidence per year*/
       for(jk=1; jk <=m; jk++) {
-        fprintf(ficgp,"\nset out \"%s%d%d.png\" \n",subdirf2(optionfilefiname,"pe"),jk,ng);
+        fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"pe"),jk,ng);
         if (ng==2)
           fprintf(ficgp,"\nset ylabel \"Quasi-incidence per year\"\n");
         else
           fprintf(ficgp,"\nset title \"Probability\"\n");
-        fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65\nset log y\nplot  [%.f:%.f] ",ageminpar,agemaxpar);
+        fprintf(ficgp,"\nset ter png small size 320, 240\nset log y\nplot  [%.f:%.f] ",ageminpar,agemaxpar);
         i=1;
         for(k2=1; k2<=nlstate; k2++) {
           k3=i;
- Line 3853  plot [%.f:%.f] \"%s\" u ($1==%d ? ($3):1
+ Line 4339  plot [%.f:%.f] \"%s\" u ($1==%d ? ($3):1
                 fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1);
               else
                 fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
-              ij=1;
+              ij=1;/* To be checked else nbcode[0][0] wrong */
               for(j=3; j <=ncovmodel; j++) {
-                if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
+                /* if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { /\* Bug valgrind *\/ */
-                  fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
+                /*        /\*fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);*\/ */
-                  ij++;
+                /*        ij++; */
-                }
+                /* } */
-                else
+                /* else */
                   fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
               }
               fprintf(ficgp,")/(1");
- Line 3868  plot [%.f:%.f] \"%s\" u ($1==%d ? ($3):1
+ Line 4354  plot [%.f:%.f] \"%s\" u ($1==%d ? ($3):1
                 fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
                 ij=1;
                 for(j=3; j <=ncovmodel; j++){
-                  if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
+                  /* if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { */
-                    fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
+                  /*   fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]); */
-                    ij++;
+                  /*   ij++; */
-                  }
+                  /* } */
-                  else
+                  /* else */
                     fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
                 }
                 fprintf(ficgp,")");
- Line 3885  plot [%.f:%.f] \"%s\" u ($1==%d ? ($3):1
+ Line 4371  plot [%.f:%.f] \"%s\" u ($1==%d ? ($3):1
         } /* end k2 */
       } /* end jk */
     } /* end ng */
+  /* avoid: */
     fflush(ficgp);
  }  /* end gnuplot */
- Line 3938  prevforecast(char fileres[], double anpr
+ Line 4425  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 yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1;
-   int *popage;
    double agec; /* generic age */
    double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
    double *popeffectif,*popcount;
- Line 4233  void prwizard(int ncovmodel, int nlstate
+ Line 4719  void prwizard(int ncovmodel, int nlstate
    /* Wizard to print covariance matrix template */
-   char ca[32], cb[32], cc[32];
+   char ca[32], cb[32];
-   int i,j, k, l, li, lj, lk, ll, jj, npar, itimes;
+   int i,j, k, li, lj, lk, ll, jj, npar, itimes;
    int numlinepar;
    printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
- Line 4386  double gompertz(double x[])
+ Line 4872  double gompertz(double x[])
    return -2*L*num/sump;
  }
+ #ifdef GSL
+ /******************* Gompertz_f Likelihood ******************************/
+ double gompertz_f(const gsl_vector *v, void *params)
+ {
+   double A,B,LL=0.0,sump=0.,num=0.;
+   double *x= (double *) v->data;
+   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;*/
+     num=num+1;
+   }
+   /* for (i=0; i<=imx; i++)
+      if (wav[i]>0) printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d %d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]);*/
+   printf("x[0]=%lf x[1]=%lf\n",x[0],x[1]);
+   for (i=1;i<=imx ; i++)
+     {
+       if (cens[i] == 1 && wav[i]>1)
+         A=-x[0]/(x[1])*(exp(x[1]*(agecens[i]-agegomp))-exp(x[1]*(ageexmed[i]-agegomp)));
+       if (cens[i] == 0 && wav[i]>1)
+         A=-x[0]/(x[1])*(exp(x[1]*(agedc[i]-agegomp))-exp(x[1]*(ageexmed[i]-agegomp)))
+              +log(x[0]/YEARM)+x[1]*(agedc[i]-agegomp)+log(YEARM);
+       /*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */
+       if (wav[i] > 1 ) { /* ??? */
+         LL=LL+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]);*/
+       }
+     }
+  /*printf("x1=%2.9f x2=%2.9f x3=%2.9f L=%f\n",x[1],x[2],x[3],L);*/
+   printf("x[0]=%lf x[1]=%lf -2*LL*num/sump=%lf\n",x[0],x[1],-2*LL*num/sump);
+   return -2*LL*num/sump;
+ }
+ #endif
  /******************* Printing html file ***********/
  void printinghtmlmort(char fileres[], char title[], char datafile[], int firstpass, \
                    int lastpass, int stepm, int weightopt, char model[],\
- Line 4414  fprintf(fichtm,"<ul><li><h4>Life table</
+ Line 4941  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;
- Line 4427  void printinggnuplotmort(char fileres[],
+ Line 4954  void printinggnuplotmort(char fileres[],
    strcpy(optfileres,"vpl");
    fprintf(ficgp,"set out \"graphmort.png\"\n ");
    fprintf(ficgp,"set xlabel \"Age\"\n set ylabel \"Force of mortality (per year)\" \n ");
-   fprintf(ficgp, "set ter png small\n set log y\n");
+   fprintf(ficgp, "set ter png small size 320, 240\n set log y\n");
-   fprintf(ficgp, "set size 0.65,0.65\n");
+   /* fprintf(ficgp, "set size 0.65,0.65\n"); */
    fprintf(ficgp,"plot [%d:100] %lf*exp(%lf*(x-%d))",agegomp,p[1],p[2],agegomp);
  }
+ int readdata(char datafile[], int firstobs, int lastobs, int *imax)
+ {
+   /*-------- data file ----------*/
+   FILE *fic;
+   char dummy[]="                         ";
+   int i=0, j=0, n=0;
+   int linei, month, year,iout;
+   char line[MAXLINE], linetmp[MAXLINE];
+   char stra[MAXLINE], strb[MAXLINE];
+   char *stratrunc;
+   int lstra;
+   if((fic=fopen(datafile,"r"))==NULL)    {
+     printf("Problem while opening datafile: %s\n", datafile);return 1;
+     fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);return 1;
+   }
+   i=1;
+   linei=0;
+   while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs))) {
+     linei=linei+1;
+     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)||(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);
+       continue;
+     }
+     trimbb(linetmp,line); /* Trims multiple blanks in line */
+     strcpy(line, linetmp);
+     for (j=maxwav;j>=1;j--){
+       cutv(stra, strb, line, ' ');
+       if(strb[0]=='.') { /* Missing status */
+         lval=-1;
+       }else{
+         errno=0;
+         lval=strtol(strb,&endptr,10);
+       /*        if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
+         if( strb[0]=='\0' || (*endptr != '\0')){
+           printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);
+           fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog);
+           return 1;
+         }
+       }
+       s[j][i]=lval;
+       strcpy(line,stra);
+       cutv(stra, strb,line,' ');
+       if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
+       }
+       else  if(iout=sscanf(strb,"%s.",dummy) != 0){
+         month=99;
+         year=9999;
+       }else{
+         printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d.  Exiting.\n",strb, linei,i, line,j);
+         fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d.  Exiting.\n",strb, linei,i, line,j);fflush(ficlog);
+         return 1;
+       }
+       anint[j][i]= (double) year;
+       mint[j][i]= (double)month;
+       strcpy(line,stra);
+     } /* ENd Waves */
+     cutv(stra, strb,line,' ');
+     if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
+     }
+     else  if(iout=sscanf(strb,"%s.",dummy) != 0){
+       month=99;
+       year=9999;
+     }else{
+       printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);
+         fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);fflush(ficlog);
+         return 1;
+     }
+     andc[i]=(double) year;
+     moisdc[i]=(double) month;
+     strcpy(line,stra);
+     cutv(stra, strb,line,' ');
+     if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
+     }
+     else  if(iout=sscanf(strb,"%s.", dummy) != 0){
+       month=99;
+       year=9999;
+     }else{
+       printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);
+       fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);fflush(ficlog);
+         return 1;
+     }
+     if (year==9999) {
+       printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given.  Exiting.\n",strb, linei,i,line);
+       fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given. Exiting.\n",strb, linei,i,line);fflush(ficlog);
+         return 1;
+     }
+     annais[i]=(double)(year);
+     moisnais[i]=(double)(month);
+     strcpy(line,stra);
+     cutv(stra, strb,line,' ');
+     errno=0;
+     dval=strtod(strb,&endptr);
+     if( strb[0]=='\0' || (*endptr != '\0')){
+       printf("Error reading data around '%f' at line number %d, \"%s\" for individual %d\nShould be a weight.  Exiting.\n",dval, i,line,linei);
+       fprintf(ficlog,"Error reading data around '%f' at line number %d, \"%s\" for individual %d\nShould be a weight.  Exiting.\n",dval, i,line,linei);
+       fflush(ficlog);
+       return 1;
+     }
+     weight[i]=dval;
+     strcpy(line,stra);
+     for (j=ncovcol;j>=1;j--){
+       cutv(stra, strb,line,' ');
+       if(strb[0]=='.') { /* Missing status */
+         lval=-1;
+       }else{
+         errno=0;
+         lval=strtol(strb,&endptr,10);
+         if( strb[0]=='\0' || (*endptr != '\0')){
+           printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative).  Exiting.\n",lval, linei,i, line);
+           fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative).  Exiting.\n",lval, linei,i, line);fflush(ficlog);
+           return 1;
+         }
+       }
+       if(lval <-1 || lval >1){
+         printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
+  Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
+  for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
+  For example, for multinomial values like 1, 2 and 3,\n \
+  build V1=0 V2=0 for the reference value (1),\n \
+         V1=1 V2=0 for (2) \n \
+  and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
+  output of IMaCh is often meaningless.\n \
+  Exiting.\n",lval,linei, i,line,j);
+         fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
+  Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
+  for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
+  For example, for multinomial values like 1, 2 and 3,\n \
+  build V1=0 V2=0 for the reference value (1),\n \
+         V1=1 V2=0 for (2) \n \
+  and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
+  output of IMaCh is often meaningless.\n \
+  Exiting.\n",lval,linei, i,line,j);fflush(ficlog);
+         return 1;
+       }
+       covar[j][i]=(double)(lval);
+       strcpy(line,stra);
+     }
+     lstra=strlen(stra);
+     if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */
+       stratrunc = &(stra[lstra-9]);
+       num[i]=atol(stratrunc);
+     }
+     else
+       num[i]=atol(stra);
+     /*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;}*/
+     i=i+1;
+   } /* End loop reading  data */
+   *imax=i-1; /* Number of individuals */
+   fclose(fic);
+   return (0);
+   /* endread: */
+     printf("Exiting readdata: ");
+     fclose(fic);
+     return (1);
+ }
+ void removespace(char *str) {
+   char *p1 = str, *p2 = str;
+   do
+     while (*p2 == ' ')
+       p2++;
+   while (*p1++ = *p2++);
+ }
+ int decodemodel ( char model[], int lastobs) /**< This routine decode the model and returns:
+    * Model  V1+V2+V3+V8+V7*V8+V5*V6+V8*age+V3*age
+    * - cptcovt total number of covariates of the model nbocc(+)+1 = 8
+    * - cptcovn or number of covariates k of the models excluding age*products =6
+    * - cptcovage number of covariates with age*products =2
+    * - cptcovs number of simple covariates
+    * - Tvar[k] is the id of the kth covariate Tvar[1]@12 {1, 2, 3, 8, 10, 11, 8, 3, 7, 8, 5, 6}, thus Tvar[5=V7*V8]=10
+    *     which is a new column after the 9 (ncovcol) variables.
+    * - if k is a product Vn*Vm covar[k][i] is filled with correct values for each individual
+    * - Tprod[l] gives the kth covariates of the product Vn*Vm l=1 to cptcovprod-cptcovage
+    *    Tprod[1]@2 {5, 6}: position of first product V7*V8 is 5, and second V5*V6 is 6.
+    * - Tvard[k]  p Tvard[1][1]@4 {7, 8, 5, 6} for V7*V8 and V5*V6 .
+  */
+ {
+   int i, j, k, ks;
+   int  j1, k1, k2;
+   char modelsav[80];
+   char stra[80], strb[80], strc[80], strd[80],stre[80];
+   /*removespace(model);*/
+   if (strlen(model) >1){ /* If there is at least 1 covariate */
+     j=0, j1=0, k1=0, k2=-1, ks=0, cptcovn=0;
+     j=nbocc(model,'+'); /**< j=Number of '+' */
+     j1=nbocc(model,'*'); /**< j1=Number of '*' */
+     cptcovs=j+1-j1; /**<  Number of simple covariates V1+V2*age+V3 +V3*V4=> V1 + V3 =2  */
+     cptcovt= j+1; /* Number of total covariates in the model V1 + V2*age+ V3 + V3*V4=> 4*/
+                   /* including age products which are counted in cptcovage.
+                   /* but the covariates which are products must be treated separately: ncovn=4- 2=2 (V1+V3). */
+     cptcovprod=j1; /**< Number of products  V1*V2 +v3*age = 2 */
+     cptcovprodnoage=0; /**< Number of covariate products without age: V3*V4 =1  */
+     strcpy(modelsav,model);
+     if (strstr(model,"AGE") !=0){
+       printf("Error. AGE must be in lower case 'age' model=%s ",model);
+       fprintf(ficlog,"Error. AGE must be in lower case model=%s ",model);fflush(ficlog);
+       return 1;
+     }
+     if (strstr(model,"v") !=0){
+       printf("Error. 'v' must be in upper case 'V' model=%s ",model);
+       fprintf(ficlog,"Error. 'v' must be in upper case model=%s ",model);fflush(ficlog);
+       return 1;
+     }
+     /*   Design
+      *  V1   V2   V3   V4  V5  V6  V7  V8  V9 Weight
+      *  <          ncovcol=8                >
+      * Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8
+      *   k=  1    2      3       4     5       6      7        8
+      *  cptcovn number of covariates (not including constant and age ) = # of + plus 1 = 7+1=8
+      *  covar[k,i], value of kth covariate if not including age for individual i:
+      *       covar[1][i]= (V2), covar[4][i]=(V3), covar[8][i]=(V8)
+      *  Tvar[k] # of the kth covariate:  Tvar[1]=2  Tvar[4]=3 Tvar[8]=8
+      *       if multiplied by age: V3*age Tvar[3=V3*age]=3 (V3) Tvar[7]=8 and
+      *  Tage[++cptcovage]=k
+      *       if products, new covar are created after ncovcol with k1
+      *  Tvar[k]=ncovcol+k1; # of the kth covariate product:  Tvar[5]=ncovcol+1=10  Tvar[6]=ncovcol+1=11
+      *  Tprod[k1]=k; Tprod[1]=5 Tprod[2]= 6; gives the position of the k1th product
+      *  Tvard[k1][1]=m Tvard[k1][2]=m; Tvard[1][1]=5 (V5) Tvard[1][2]=6 Tvard[2][1]=7 (V7) Tvard[2][2]=8
+      *  Tvar[cptcovn+k2]=Tvard[k1][1];Tvar[cptcovn+k2+1]=Tvard[k1][2];
+      *  Tvar[8+1]=5;Tvar[8+2]=6;Tvar[8+3]=7;Tvar[8+4]=8 inverted
+      *  V1   V2   V3   V4  V5  V6  V7  V8  V9  V10  V11
+      *  <          ncovcol=8                >
+      *       Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8    d1   d1   d2  d2
+      *          k=  1    2      3       4     5       6      7        8    9   10   11  12
+      *     Tvar[k]= 2    1      3       3    10      11      8        8    5    6    7   8
+      * p Tvar[1]@12={2,   1,     3,      3,   11,     10,     8,       8,   7,   8,   5,  6}
+      * p Tprod[1]@2={                         6, 5}
+      *p Tvard[1][1]@4= {7, 8, 5, 6}
+      * covar[k][i]= V2   V1      ?      V3    V5*V6?   V7*V8?  ?       V8
+      *  cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
+      *How to reorganize?
+      * Model V1 + V2 + V3 + V8 + V5*V6 + V7*V8 + V3*age + V8*age
+      * Tvars {2,   1,     3,      3,   11,     10,     8,       8,   7,   8,   5,  6}
+      *       {2,   1,     4,      8,    5,      6,     3,       7}
+      * Struct []
+      */
+     /* This loop fills the array Tvar from the string 'model'.*/
+     /* j is the number of + signs in the model V1+V2+V3 j=2 i=3 to 1 */
+     /*   modelsav=V2+V1+V4+age*V3 strb=age*V3 stra=V2+V1+V4  */
+     /*  k=4 (age*V3) Tvar[k=4]= 3 (from V3) Tage[cptcovage=1]=4 */
+     /*  k=3 V4 Tvar[k=3]= 4 (from V4) */
+     /*  k=2 V1 Tvar[k=2]= 1 (from V1) */
+     /*  k=1 Tvar[1]=2 (from V2) */
+     /*  k=5 Tvar[5] */
+     /* for (k=1; k<=cptcovn;k++) { */
+     /*  cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; */
+     /*  } */
+     /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
+     /*
+      * Treating invertedly V2+V1+V3*age+V2*V4 is as if written V2*V4 +V3*age + V1 + V2 */
+     for(k=cptcovt; k>=1;k--) /**< Number of covariates */
+         Tvar[k]=0;
+     cptcovage=0;
+     for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model */
+       cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+'
+                                      modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 */
+       if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */
+       /*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
+       /*scanf("%d",i);*/
+       if (strchr(strb,'*')) {  /**< Model includes a product V2+V1+V4+V3*age strb=V3*age */
+         cutl(strc,strd,strb,'*'); /**< strd*strc  Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */
+         if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */
+           /* covar is not filled and then is empty */
+           cptcovprod--;
+           cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */
+           Tvar[k]=atoi(stre);  /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2 */
+           cptcovage++; /* Sums the number of covariates which include age as a product */
+           Tage[cptcovage]=k;  /* Tage[1] = 4 */
+           /*printf("stre=%s ", stre);*/
+         } else if (strcmp(strd,"age")==0) { /* or age*Vn */
+           cptcovprod--;
+           cutl(stre,strb,strc,'V');
+           Tvar[k]=atoi(stre);
+           cptcovage++;
+           Tage[cptcovage]=k;
+         } else {  /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2  strb=V3*V2*/
+           /* loops on k1=1 (V3*V2) and k1=2 V4*V3 */
+           cptcovn++;
+           cptcovprodnoage++;k1++;
+           cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
+           Tvar[k]=ncovcol+k1; /* For model-covariate k tells which data-covariate to use but
+                                   because this model-covariate is a construction we invent a new column
+                                   ncovcol + k1
+                                   If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2
+                                   Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */
+           cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */
+           Tprod[k1]=k;  /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2  */
+           Tvard[k1][1] =atoi(strc); /* m 1 for V1*/
+           Tvard[k1][2] =atoi(stre); /* n 4 for V4*/
+           k2=k2+2;
+           Tvar[cptcovt+k2]=Tvard[k1][1]; /* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) */
+           Tvar[cptcovt+k2+1]=Tvard[k1][2];  /* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) */
+           for (i=1; i<=lastobs;i++){
+             /* Computes the new covariate which is a product of
+                covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */
+             covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];
+           }
+         } /* End age is not in the model */
+       } /* End if model includes a product */
+       else { /* no more sum */
+         /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
+        /*  scanf("%d",i);*/
+         cutl(strd,strc,strb,'V');
+         ks++; /**< Number of simple covariates */
+         cptcovn++;
+         Tvar[k]=atoi(strd);
+       }
+       strcpy(modelsav,stra);  /* modelsav=V2+V1+V4 stra=V2+V1+V4 */
+       /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
+         scanf("%d",i);*/
+     } /* end of loop + */
+   } /* end model */
+   /*The number n of Vn is stored in Tvar. cptcovage =number of age covariate. Tage gives the position of age. cptcovprod= number of products.
+     If model=V1+V1*age then Tvar[1]=1 Tvar[2]=1 cptcovage=1 Tage[1]=2 cptcovprod=0*/
+   /* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);
+   printf("cptcovprod=%d ", cptcovprod);
+   fprintf(ficlog,"cptcovprod=%d ", cptcovprod);
+   scanf("%d ",i);*/
+   return (0); /* with covar[new additional covariate if product] and Tage if age */
+   /*endread:*/
+     printf("Exiting decodemodel: ");
+     return (1);
+ }
+ calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn )
+ {
+   int i, m;
+   for (i=1; i<=imx; i++) {
+     for(m=2; (m<= maxwav); m++) {
+       if (((int)mint[m][i]== 99) && (s[m][i] <= nlstate)){
+         anint[m][i]=9999;
+         s[m][i]=-1;
+       }
+       if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){
+         *nberr++;
+         printf("Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i);
+         fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\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){
+         *nberr++;
+         printf("Error! Month of death of individual %ld on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]);
+         fprintf(ficlog,"Error! Month of death of individual %ld on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]);
+         s[m][i]=-1; /* We prefer to skip it (and to skip it in version 0.8a1 too */
+       }
+     }
+   }
+   for (i=1; i<=imx; i++)  {
+     agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
+     for(m=firstpass; (m<= lastpass); m++){
+       if(s[m][i] >0 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){
+         if (s[m][i] >= nlstate+1) {
+           if(agedc[i]>0)
+             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 ((int)andc[i]!=9999){
+                 nbwarn++;
+                 printf("Warning negative age at death: %ld line:%d\n",num[i],i);
+                 fprintf(ficlog,"Warning negative age at death: %ld line:%d\n",num[i],i);
+                 agev[m][i]=-1;
+               }
+             }
+         }
+         else if(s[m][i] !=9){ /* Standard case, age in fractional
+                                  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((int)mint[m][i]==99 || (int)anint[m][i]==9999)
+             agev[m][i]=1;
+           else if(agev[m][i] < *agemin){
+             *agemin=agev[m][i];
+             printf(" Min anint[%d][%d]=%.2f annais[%d]=%.2f, agemin=%.2f\n",m,i,anint[m][i], i,annais[i], *agemin);
+           }
+           else if(agev[m][i] >*agemax){
+             *agemax=agev[m][i];
+             /* printf(" Max anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.2f\n",m,i,anint[m][i], i,annais[i], *agemax);*/
+           }
+           /*agev[m][i]=anint[m][i]-annais[i];*/
+           /*     agev[m][i] = age[i]+2*m;*/
+         }
+         else { /* =9 */
+           agev[m][i]=1;
+           s[m][i]=-1;
+         }
+       }
+       else /*= 0 Unknown */
+         agev[m][i]=1;
+     }
+   }
+   for (i=1; i<=imx; i++)  {
+     for(m=firstpass; (m<=lastpass); m++){
+       if (s[m][i] > (nlstate+ndeath)) {
+         *nberr++;
+         printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);
+         fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);
+         return 1;
+       }
+     }
+   }
+   /*for (i=1; i<=imx; i++){
+   for (m=firstpass; (m<lastpass); m++){
+      printf("%ld %d %.lf %d %d\n", num[i],(covar[1][i]),agev[m][i],s[m][i],s[m+1][i]);
+ }
+ }*/
+   printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, *agemin, *agemax);
+   fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, *agemin, *agemax);
+   return (0);
+  /* endread:*/
+     printf("Exiting calandcheckages: ");
+     return (1);
+ }
  /***********************************************/
  /**************** Main Program *****************/
- Line 4443  void printinggnuplotmort(char fileres[],
+ Line 5425  void printinggnuplotmort(char fileres[],
  int main(int argc, char *argv[])
  {
+ #ifdef GSL
+   const gsl_multimin_fminimizer_type *T;
+   size_t iteri = 0, it;
+   int rval = GSL_CONTINUE;
+   int status = GSL_SUCCESS;
+   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 i,j, k, n=MAXN,iter=0,m,size=100, cptcod;
-   int linei, month, year,iout;
-   int jj, ll, li, lj, lk, imk;
+   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];
-   char dummy[]="                         ";
    /*  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 linetmp[MAXLINE];
+   char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[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;
-   int c,  h , cpt,l;
+   int jl;
-   int ju,jl, mi;
+   int i1, j1, jk, stepsize;
-   int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;
+   int *tab;
-   int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,*tab;
    int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */
    int mobilav=0,popforecast=0;
    int hstepm, nhstepm;
- Line 4486  int main(int argc, char *argv[])
+ Line 5470  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 *severity;
    double ***param; /* Matrix of parameters */
    double  *p;
    double **matcov; /* Matrix of covariance */
- Line 4498  int main(int argc, char *argv[])
+ Line 5481  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];
+   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];
-   char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80];
-   char  *strt, strtend[80];
-   char *stratrunc;
-   int lstra;
-   long total_usecs;
  /*   setlocale (LC_ALL, ""); */
  /*   bindtextdomain (PACKAGE, LOCALEDIR); */
  /*   textdomain (PACKAGE); */
- Line 4520  int main(int argc, char *argv[])
+ Line 5500  int main(int argc, char *argv[])
  /*   setlocale (LC_MESSAGES, ""); */
    /*   gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */
-   (void) gettimeofday(&start_time,&tzp);
+   rstart_time = time(NULL);
+   /*  (void) gettimeofday(&start_time,&tzp);*/
+   start_time = *localtime(&rstart_time);
    curr_time=start_time;
-   tm = *localtime(&start_time.tv_sec);
+   /*tml = *localtime(&start_time.tm_sec);*/
-   tmg = *gmtime(&start_time.tv_sec);
+   /* strcpy(strstart,asctime(&tml)); */
-   strcpy(strstart,asctime(&tm));
+   strcpy(strstart,asctime(&start_time));
  /*  printf("Localtime (at start)=%s",strstart); */
- /*  tp.tv_sec = tp.tv_sec +86400; */
+ /*  tp.tm_sec = tp.tm_sec +86400; */
- /*  tm = *localtime(&start_time.tv_sec); */
+ /*  tm = *localtime(&start_time.tm_sec); */
  /*   tmg.tm_year=tmg.tm_year +dsign*dyear; */
  /*   tmg.tm_mon=tmg.tm_mon +dsign*dmonth; */
  /*   tmg.tm_hour=tmg.tm_hour + 1; */
- /*   tp.tv_sec = mktime(&tmg); */
+ /*   tp.tm_sec = mktime(&tmg); */
  /*   strt=asctime(&tmg); */
  /*   printf("Time(after) =%s",strstart);  */
  /*  (void) time (&time_value);
- Line 4553  int main(int argc, char *argv[])
+ Line 5535  int main(int argc, char *argv[])
      i=strlen(pathr);
      if(pathr[i-1]=='\n')
        pathr[i-1]='\0';
+     i=strlen(pathr);
+     if(pathr[i-1]==' ') /* This may happen when dragging on oS/X! */
+       pathr[i-1]='\0';
     for (tok = pathr; tok != NULL; ){
        printf("Pathr |%s|\n",pathr);
        while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0');
- Line 4608  int main(int argc, char *argv[])
+ Line 5593  int main(int argc, char *argv[])
   path=%s \n\
   optionfile=%s\n\
   optionfilext=%s\n\
-  optionfilefiname=%s\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname);
+  optionfilefiname='%s'\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname);
    printf("Local time (at start):%s",strstart);
    fprintf(ficlog,"Local time (at start): %s",strstart);
    fflush(ficlog);
  /*   (void) gettimeofday(&curr_time,&tzp); */
- /*   printf("Elapsed time %d\n", asc_diff_time(curr_time.tv_sec-start_time.tv_sec,tmpout)); */
+ /*   printf("Elapsed time %d\n", asc_diff_time(curr_time.tm_sec-start_time.tm_sec,tmpout)); */
    /* */
    strcpy(fileres,"r");
- Line 4624  int main(int argc, char *argv[])
+ Line 5609  int main(int argc, char *argv[])
    /*---------arguments file --------*/
    if((ficpar=fopen(optionfile,"r"))==NULL)    {
-     printf("Problem with optionfile %s\n",optionfile);
+     printf("Problem with optionfile '%s' with errno='%s'\n",optionfile,strerror(errno));
-     fprintf(ficlog,"Problem with optionfile %s\n",optionfile);
+     fprintf(ficlog,"Problem with optionfile '%s' with errno='%s'\n",optionfile,strerror(errno));
      fflush(ficlog);
-     goto end;
+     /* goto end; */
+     exit(70);
    }
- Line 4647  int main(int argc, char *argv[])
+ Line 5633  int main(int argc, char *argv[])
      ungetc(c,ficpar);
      fgets(line, MAXLINE, ficpar);
      numlinepar++;
-     puts(line);
+     fputs(line,stdout);
      fputs(line,ficparo);
      fputs(line,ficlog);
    }
- Line 4663  int main(int argc, char *argv[])
+ Line 5649  int main(int argc, char *argv[])
      ungetc(c,ficpar);
      fgets(line, MAXLINE, ficpar);
      numlinepar++;
-     puts(line);
+     fputs(line, stdout);
+     //puts(line);
      fputs(line,ficparo);
      fputs(line,ficlog);
    }
    ungetc(c,ficpar);
-   covar=matrix(0,NCOVMAX,1,n);
+   covar=matrix(0,NCOVMAX,1,n);  /**< used in readdata */
-   cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement*/
+   cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement plus one, indepently of n in Vn*/
-   if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;
+   /* v1+v2+v3+v2*v4+v5*age makes cptcovn = 5
-   /* where is ncovprod ?*/
+      v1+v2*age+v2*v3 makes cptcovn = 3
-   ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age : v1+v2+v3+v2*v4+v5*age makes 5+2=7*/
+   */
+   if (strlen(model)>1)
+     ncovmodel=2+nbocc(model,'+')+1; /*Number of variables including intercept and age = cptcovn + intercept + age : v1+v2+v3+v2*v4+v5*age makes 5+2=7*/
+   else
+     ncovmodel=2;
    nvar=ncovmodel-1; /* Suppressing age as a basic covariate */
    nforce= (nlstate+ndeath-1)*nlstate; /* Number of forces ij from state i to j */
    npar= nforce*ncovmodel; /* Number of parameters like aij*/
- Line 4706  int main(int argc, char *argv[])
+ Line 5697  int main(int argc, char *argv[])
      matcov=matrix(1,npar,1,npar);
    }
    else{
-     /* Read guess parameters */
+     /* Read guessed parameters */
      /* Reads comments: lines beginning with '#' */
      while((c=getc(ficpar))=='#' && c!= EOF){
        ungetc(c,ficpar);
        fgets(line, MAXLINE, ficpar);
        numlinepar++;
-       puts(line);
+       fputs(line,stdout);
        fputs(line,ficparo);
        fputs(line,ficlog);
      }
- Line 4755  run imach with mle=-1 to get a correct t
+ Line 5746  run imach with mle=-1 to get a correct t
      }
      fflush(ficlog);
+     /* Reads scales values */
      p=param[1][1];
      /* Reads comments: lines beginning with '#' */
- Line 4762  run imach with mle=-1 to get a correct t
+ Line 5754  run imach with mle=-1 to get a correct t
        ungetc(c,ficpar);
        fgets(line, MAXLINE, ficpar);
        numlinepar++;
-       puts(line);
+       fputs(line,stdout);
        fputs(line,ficparo);
        fputs(line,ficlog);
      }
- Line 4771  run imach with mle=-1 to get a correct t
+ Line 5763  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);
          }
- Line 4784  run imach with mle=-1 to get a correct t
+ Line 5776  run imach with mle=-1 to get a correct t
            fprintf(ficparo," %le",delti3[i][j][k]);
            fprintf(ficlog," %le",delti3[i][j][k]);
          }
          fscanf(ficpar,"\n");
          numlinepar++;
          printf("\n");
          fprintf(ficparo,"\n");
          fprintf(ficlog,"\n");
-       }
-     }
-     fflush(ficlog);
-     delti=delti3[1][1];
-     /* free_ma3x(delti3,1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); */ /* Hasn't to to freed here otherwise delti is no more allocated */
-     /* Reads comments: lines beginning with '#' */
-     while((c=getc(ficpar))=='#' && c!= EOF){
-       ungetc(c,ficpar);
-       fgets(line, MAXLINE, ficpar);
-       numlinepar++;
-       puts(line);
-       fputs(line,ficparo);
-       fputs(line,ficlog);
-     }
-     ungetc(c,ficpar);
-     matcov=matrix(1,npar,1,npar);
-     for(i=1; i <=npar; i++)
-       for(j=1; j <=npar; j++) matcov[i][j]=0.;
-     for(i=1; i <=npar; i++){
-       fscanf(ficpar,"%s",&str);
-       if(mle==1)
-         printf("%s",str);
-       fprintf(ficlog,"%s",str);
-       fprintf(ficparo,"%s",str);
-       for(j=1; j <=i; j++){
-         fscanf(ficpar," %le",&matcov[i][j]);
-         if(mle==1){
-           printf(" %.5le",matcov[i][j]);
-         }
-         fprintf(ficlog," %.5le",matcov[i][j]);
-         fprintf(ficparo," %.5le",matcov[i][j]);
-       }
-       fscanf(ficpar,"\n");
-       numlinepar++;
-       if(mle==1)
-         printf("\n");
-       fprintf(ficlog,"\n");
-       fprintf(ficparo,"\n");
-     }
-     for(i=1; i <=npar; i++)
-       for(j=i+1;j<=npar;j++)
-         matcov[i][j]=matcov[j][i];
-     if(mle==1)
-       printf("\n");
-     fprintf(ficlog,"\n");
-     fflush(ficlog);
-     /*-------- Rewriting parameter file ----------*/
-     strcpy(rfileres,"r");    /* "Rparameterfile */
-     strcat(rfileres,optionfilefiname);    /* Parameter file first name*/
-     strcat(rfileres,".");    /* */
-     strcat(rfileres,optionfilext);    /* Other files have txt extension */
-     if((ficres =fopen(rfileres,"w"))==NULL) {
-       printf("Problem writing new parameter file: %s\n", fileres);goto end;
-       fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end;
-     }
-     fprintf(ficres,"#%s\n",version);
-   }    /* End of mle != -3 */
-   /*-------- data file ----------*/
-   if((fic=fopen(datafile,"r"))==NULL)    {
-     printf("Problem while opening datafile: %s\n", datafile);goto end;
-     fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);goto end;
-   }
-   n= lastobs;
-   severity = vector(1,maxwav);
-   outcome=imatrix(1,maxwav+1,1,n);
-   num=lvector(1,n);
-   moisnais=vector(1,n);
-   annais=vector(1,n);
-   moisdc=vector(1,n);
-   andc=vector(1,n);
-   agedc=vector(1,n);
-   cod=ivector(1,n);
-   weight=vector(1,n);
-   for(i=1;i<=n;i++) weight[i]=1.0; /* Equal weights, 1 by default */
-   mint=matrix(1,maxwav,1,n);
-   anint=matrix(1,maxwav,1,n);
-   s=imatrix(1,maxwav+1,1,n); /* s[i][j] health state for wave i and individual j */
-   tab=ivector(1,NCOVMAX);
-   ncodemax=ivector(1,8);
-   i=1;
-   linei=0;
-   while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs))) {
-     linei=linei+1;
-     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)||(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);
-       continue;
-     }
-     trimbb(linetmp,line); /* Trims multiple blanks in line */
-     for (j=0; line[j]!='\0';j++){
-       line[j]=linetmp[j];
-     }
-     for (j=maxwav;j>=1;j--){
-       cutv(stra, strb,line,' ');
-       if(strb[0]=='.') { /* Missing status */
-         lval=-1;
-       }else{
-         errno=0;
-         lval=strtol(strb,&endptr,10);
-       /*        if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
-         if( strb[0]=='\0' || (*endptr != '\0')){
-           printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);
-           fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong.  Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog);
-           goto end;
-         }
-       }
-       s[j][i]=lval;
-       strcpy(line,stra);
-       cutv(stra, strb,line,' ');
-       if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
-       }
-       else  if(iout=sscanf(strb,"%s.") != 0){
-         month=99;
-         year=9999;
-       }else{
-         printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d.  Exiting.\n",strb, linei,i, line,j);
-         fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d.  Exiting.\n",strb, linei,i, line,j);fflush(ficlog);
-         goto end;
-       }
-       anint[j][i]= (double) year;
-       mint[j][i]= (double)month;
-       strcpy(line,stra);
-     } /* ENd Waves */
-     cutv(stra, strb,line,' ');
-     if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
-     }
-     else  if(iout=sscanf(strb,"%s.",dummy) != 0){
-       month=99;
-       year=9999;
-     }else{
-       printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of death (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);
-         fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of death (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);fflush(ficlog);
-         goto end;
-     }
-     andc[i]=(double) year;
-     moisdc[i]=(double) month;
-     strcpy(line,stra);
-     cutv(stra, strb,line,' ');
-     if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
-     }
-     else  if(iout=sscanf(strb,"%s.") != 0){
-       month=99;
-       year=9999;
-     }else{
-       printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);
-       fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);fflush(ficlog);
-         goto end;
-     }
-     if (year==9999) {
-       printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given.  Exiting.\n",strb, linei,i,line);
-       fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given. Exiting.\n",strb, linei,i,line);fflush(ficlog);
-         goto end;
-     }
-     annais[i]=(double)(year);
-     moisnais[i]=(double)(month);
-     strcpy(line,stra);
-     cutv(stra, strb,line,' ');
-     errno=0;
-     dval=strtod(strb,&endptr);
-     if( strb[0]=='\0' || (*endptr != '\0')){
-       printf("Error reading data around '%f' at line number %ld, \"%s\" for individual %d\nShould be a weight.  Exiting.\n",dval, i,line,linei);
-       fprintf(ficlog,"Error reading data around '%f' at line number %ld, \"%s\" for individual %d\nShould be a weight.  Exiting.\n",dval, i,line,linei);
-       fflush(ficlog);
-       goto end;
-     }
-     weight[i]=dval;
-     strcpy(line,stra);
-     for (j=ncovcol;j>=1;j--){
-       cutv(stra, strb,line,' ');
-       if(strb[0]=='.') { /* Missing status */
-         lval=-1;
-       }else{
-         errno=0;
-         lval=strtol(strb,&endptr,10);
-         if( strb[0]=='\0' || (*endptr != '\0')){
-           printf("Error reading data around '%d' at line number %ld for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative).  Exiting.\n",lval, linei,i, line);
-           fprintf(ficlog,"Error reading data around '%d' at line number %ld for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative).  Exiting.\n",lval, linei,i, line);fflush(ficlog);
-           goto end;
-         }
-       }
-       if(lval <-1 || lval >1){
-         printf("Error reading data around '%d' at line number %ld for individual %d, '%s'\n \
-  Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
-  for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
-  For example, for multinomial values like 1, 2 and 3,\n \
-  build V1=0 V2=0 for the reference value (1),\n \
-         V1=1 V2=0 for (2) \n \
-  and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
-  output of IMaCh is often meaningless.\n \
-  Exiting.\n",lval,linei, i,line,j);
-         fprintf(ficlog,"Error reading data around '%d' at line number %ld for individual %d, '%s'\n \
-  Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
-  for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
-  For example, for multinomial values like 1, 2 and 3,\n \
-  build V1=0 V2=0 for the reference value (1),\n \
-         V1=1 V2=0 for (2) \n \
-  and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
-  output of IMaCh is often meaningless.\n \
-  Exiting.\n",lval,linei, i,line,j);fflush(ficlog);
-         goto end;
-       }
-       covar[j][i]=(double)(lval);
-       strcpy(line,stra);
-     }
-     lstra=strlen(stra);
-     if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */
-       stratrunc = &(stra[lstra-9]);
-       num[i]=atol(stratrunc);
-     }
-     else
-       num[i]=atol(stra);
-     /*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;}*/
-     i=i+1;
-   } /* End loop reading  data */
-   fclose(fic);
-   /* printf("ii=%d", ij);
-      scanf("%d",i);*/
-   imx=i-1; /* Number of individuals */
-   /* for (i=1; i<=imx; i++){
-     if ((s[1][i]==3) && (s[2][i]==2)) s[2][i]=3;
-     if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3;
-     if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;
-     }*/
-    /*  for (i=1; i<=imx; i++){
-      if (s[4][i]==9)  s[4][i]=-1;
-      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]));}*/
-   /* for (i=1; i<=imx; i++) */
-    /*if ((s[3][i]==3) ||  (s[4][i]==3)) weight[i]=0.08;
-      else weight[i]=1;*/
-   /* Calculation of the number of parameters from char model */
-   Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. Stores the number n of the covariates in Vm+Vn at 1 and m at 2 */
-   Tprod=ivector(1,15);
-   Tvaraff=ivector(1,15);
-   Tvard=imatrix(1,15,1,2);
-   Tage=ivector(1,15);
-   if (strlen(model) >1){ /* If there is at least 1 covariate */
-     j=0, j1=0, k1=1, k2=1;
-     j=nbocc(model,'+'); /* j=Number of '+' */
-     j1=nbocc(model,'*'); /* j1=Number of '*' */
-     cptcovn=j+1; /* Number of covariates V1+V2+V3 =>2+1=3 */
-     cptcovprod=j1; /*Number of products  V1*V2 =1 */
-     strcpy(modelsav,model);
-     if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){
-       printf("Error. Non available option model=%s ",model);
-       fprintf(ficlog,"Error. Non available option model=%s ",model);fflush(ficlog);
-       goto end;
-     }
-     /* This loop fills the array Tvar from the string 'model'.*/
-     /* j is the number of + signs in the model V1+V2+V3 j=2 i=3 to 1 */
-     for(i=(j+1); i>=1;i--){
-       cutv(stra,strb,modelsav,'+'); /* keeps in strb after the first '+'
-                                      modelsav=V2+V3*age+V1+V4 strb=V3*age+V1+V4
-                                      stra=V2
-                                     */
-       if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */
-       /*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
-       /*scanf("%d",i);*/
-       if (strchr(strb,'*')) {  /* Model includes a product V1+V3*age+V2 strb=V3*age*/
-         cutv(strd,strc,strb,'*'); /* strd*strc  Vm*Vn: V3*age strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */
-         if (strcmp(strc,"age")==0) { /* Vn*age */
-           cptcovprod--;
-           cutv(strb,stre,strd,'V');
-           Tvar[i]=atoi(stre);  /* V1+V3*age+V2 Tvar[2]=3 */
-           cptcovage++; /* Sums the number of covariates including age as a product */
-           Tage[cptcovage]=i;  /* Tage[1] =2 */
-           /*printf("stre=%s ", stre);*/
-         }
-         else if (strcmp(strd,"age")==0) { /* or age*Vn */
-           cptcovprod--;
-           cutv(strb,stre,strc,'V');
-           Tvar[i]=atoi(stre);
-           cptcovage++;
-           Tage[cptcovage]=i;
-         }
-         else {  /* Age is not in the model V1+V3*V2+V2  strb=V3*V2*/
-           cutv(strb,stre,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
-           Tvar[i]=ncovcol+k1;  /* find 'n' in Vn and stores in Tvar.
-                                   If already ncovcol=2 and model=V2*V1 Tvar[1]=2+1 and Tvar[2]=2+2 etc */
-           cutv(strb,strc,strd,'V'); /* strd was Vm, strc is m */
-           Tprod[k1]=i;  /* Tprod[1]  */
-           Tvard[k1][1]=atoi(strc); /* m*/
-           Tvard[k1][2]=atoi(stre); /* n */
-           Tvar[cptcovn+k2]=Tvard[k1][1];
-           Tvar[cptcovn+k2+1]=Tvard[k1][2];
-           for (k=1; k<=lastobs;k++)
-             covar[ncovcol+k1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];
-           k1++;
-           k2=k2+2;
-         }
-       }
-       else { /* no more sum */
-         /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
-        /*  scanf("%d",i);*/
-       cutv(strd,strc,strb,'V');
-       Tvar[i]=atoi(strc);
        }
-       strcpy(modelsav,stra);  /* modelsav=V2+V3*age+V1+V4 strb=V3*age+V1+V4 */
+     }
-       /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
+     fflush(ficlog);
-         scanf("%d",i);*/
-     } /* end of loop + */
-   } /* end model */
-   /*The number n of Vn is stored in Tvar. cptcovage =number of age covariate. Tage gives the position of age. cptcovprod= number of products.
-     If model=V1+V1*age then Tvar[1]=1 Tvar[2]=1 cptcovage=1 Tage[1]=2 cptcovprod=0*/
-   /* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);
-   printf("cptcovprod=%d ", cptcovprod);
-   fprintf(ficlog,"cptcovprod=%d ", cptcovprod);
-   scanf("%d ",i);*/
+     /* Reads covariance matrix */
+     delti=delti3[1][1];
-     /*  if(mle==1){*/
-   if (weightopt != 1) { /* Maximisation without weights*/
-     for(i=1;i<=n;i++) weight[i]=1.0;
-   }
-     /*-calculation of age at interview from date of interview and age at death -*/
-   agev=matrix(1,maxwav,1,imx);
-   for (i=1; i<=imx; i++) {
+     /* free_ma3x(delti3,1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); */ /* Hasn't to to freed here otherwise delti is no more allocated */
-     for(m=2; (m<= maxwav); m++) {
-       if (((int)mint[m][i]== 99) && (s[m][i] <= nlstate)){
+     /* Reads comments: lines beginning with '#' */
-         anint[m][i]=9999;
+     while((c=getc(ficpar))=='#' && c!= EOF){
-         s[m][i]=-1;
+       ungetc(c,ficpar);
-       }
+       fgets(line, MAXLINE, ficpar);
-       if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){
+       numlinepar++;
-         nberr++;
+       fputs(line,stdout);
-         printf("Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i);
+       fputs(line,ficparo);
-         fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i);
+       fputs(line,ficlog);
-         s[m][i]=-1;
-       }
-       if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){
-         nberr++;
-         printf("Error! Month of death of individual %ld on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]);
-         fprintf(ficlog,"Error! Month of death of individual %ld on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]);
-         s[m][i]=-1; /* We prefer to skip it (and to skip it in version 0.8a1 too */
-       }
      }
-   }
+     ungetc(c,ficpar);
-   for (i=1; i<=imx; i++)  {
+     matcov=matrix(1,npar,1,npar);
-     agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
+     for(i=1; i <=npar; i++)
-     for(m=firstpass; (m<= lastpass); m++){
+       for(j=1; j <=npar; j++) matcov[i][j]=0.;
-       if(s[m][i] >0 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){
-         if (s[m][i] >= nlstate+1) {
+     for(i=1; i <=npar; i++){
-           if(agedc[i]>0)
+       fscanf(ficpar,"%s",str);
-             if((int)moisdc[i]!=99 && (int)andc[i]!=9999)
+       if(mle==1)
-               agev[m][i]=agedc[i];
+         printf("%s",str);
-           /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/
+       fprintf(ficlog,"%s",str);
-             else {
+       fprintf(ficparo,"%s",str);
-               if ((int)andc[i]!=9999){
+       for(j=1; j <=i; j++){
-                 nbwarn++;
+         fscanf(ficpar," %le",&matcov[i][j]);
-                 printf("Warning negative age at death: %ld line:%d\n",num[i],i);
+         if(mle==1){
-                 fprintf(ficlog,"Warning negative age at death: %ld line:%d\n",num[i],i);
+           printf(" %.5le",matcov[i][j]);
-                 agev[m][i]=-1;
-               }
-             }
-         }
-         else if(s[m][i] !=9){ /* Standard case, age in fractional
-                                  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((int)mint[m][i]==99 || (int)anint[m][i]==9999)
-             agev[m][i]=1;
-           else if(agev[m][i] <agemin){
-             agemin=agev[m][i];
-             printf(" Min anint[%d][%d]=%.2f annais[%d]=%.2f, agemin=%.2f\n",m,i,anint[m][i], i,annais[i], agemin);
-           }
-           else if(agev[m][i] >agemax){
-             agemax=agev[m][i];
-             /* printf(" anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.0f\n",m,i,anint[m][i], i,annais[i], agemax);*/
-           }
-           /*agev[m][i]=anint[m][i]-annais[i];*/
-           /*     agev[m][i] = age[i]+2*m;*/
-         }
-         else { /* =9 */
-           agev[m][i]=1;
-           s[m][i]=-1;
          }
+         fprintf(ficlog," %.5le",matcov[i][j]);
+         fprintf(ficparo," %.5le",matcov[i][j]);
        }
-       else /*= 0 Unknown */
+       fscanf(ficpar,"\n");
-         agev[m][i]=1;
+       numlinepar++;
+       if(mle==1)
+         printf("\n");
+       fprintf(ficlog,"\n");
+       fprintf(ficparo,"\n");
      }
+     for(i=1; i <=npar; i++)
+       for(j=i+1;j<=npar;j++)
+         matcov[i][j]=matcov[j][i];
-   }
+     if(mle==1)
-   for (i=1; i<=imx; i++)  {
+       printf("\n");
-     for(m=firstpass; (m<=lastpass); m++){
+     fprintf(ficlog,"\n");
-       if (s[m][i] > (nlstate+ndeath)) {
-         nberr++;
+     fflush(ficlog);
-         printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);
-         fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);
+     /*-------- Rewriting parameter file ----------*/
-         goto end;
+     strcpy(rfileres,"r");    /* "Rparameterfile */
-       }
+     strcat(rfileres,optionfilefiname);    /* Parameter file first name*/
+     strcat(rfileres,".");    /* */
+     strcat(rfileres,optionfilext);    /* Other files have txt extension */
+     if((ficres =fopen(rfileres,"w"))==NULL) {
+       printf("Problem writing new parameter file: %s\n", fileres);goto end;
+       fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end;
      }
-   }
+     fprintf(ficres,"#%s\n",version);
+   }    /* End of mle != -3 */
-   /*for (i=1; i<=imx; i++){
-   for (m=firstpass; (m<lastpass); m++){
-      printf("%ld %d %.lf %d %d\n", num[i],(covar[1][i]),agev[m][i],s[m][i],s[m+1][i]);
- }
- }*/
+   n= lastobs;
+   num=lvector(1,n);
+   moisnais=vector(1,n);
+   annais=vector(1,n);
+   moisdc=vector(1,n);
+   andc=vector(1,n);
+   agedc=vector(1,n);
+   cod=ivector(1,n);
+   weight=vector(1,n);
+   for(i=1;i<=n;i++) weight[i]=1.0; /* Equal weights, 1 by default */
+   mint=matrix(1,maxwav,1,n);
+   anint=matrix(1,maxwav,1,n);
+   s=imatrix(1,maxwav+1,1,n); /* s[i][j] health state for wave i and individual j */
+   tab=ivector(1,NCOVMAX);
+   ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */
+   /* Reads data from file datafile */
+   if (readdata(datafile, firstobs, lastobs, &imx)==1)
+     goto end;
+   /* Calculation of the number of parameters from char model */
+     /*    modelsav=V2+V1+V4+age*V3 strb=age*V3 stra=V2+V1+V4
+         k=4 (age*V3) Tvar[k=4]= 3 (from V3) Tag[cptcovage=1]=4
+         k=3 V4 Tvar[k=3]= 4 (from V4)
+         k=2 V1 Tvar[k=2]= 1 (from V1)
+         k=1 Tvar[1]=2 (from V2)
+     */
+   Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. */
+   /*  V2+V1+V4+age*V3 is a model with 4 covariates (3 plus signs).
+       For each model-covariate stores the data-covariate id. Tvar[1]=2, Tvar[2]=1, Tvar[3]=4,
+       Tvar[4=age*V3] is 3 and 'age' is recorded in Tage.
+   */
+   /* For model-covariate k tells which data-covariate to use but
+     because this model-covariate is a construction we invent a new column
+     ncovcol + k1
+     If already ncovcol=4 and model=V2+V1+V1*V4+age*V3
+     Tvar[3=V1*V4]=4+1 etc */
+   Tprod=ivector(1,NCOVMAX); /* Gives the position of a product */
+   /* Tprod[k1=1]=3(=V1*V4) for V2+V1+V1*V4+age*V3
+      if  V2+V1+V1*V4+age*V3+V3*V2   TProd[k1=2]=5 (V3*V2)
+   */
+   Tvaraff=ivector(1,NCOVMAX); /* Unclear */
+   Tvard=imatrix(1,NCOVMAX,1,2); /* n=Tvard[k1][1]  and m=Tvard[k1][2] gives the couple n,m of the k1 th product Vn*Vm
+                             * For V3*V2 (in V2+V1+V1*V4+age*V3+V3*V2), V3*V2 position is 2nd.
+                             * Tvard[k1=2][1]=3 (V3) Tvard[k1=2][2]=2(V2) */
+   Tage=ivector(1,NCOVMAX); /* Gives the covariate id of covariates associated with age: V2 + V1 + age*V4 + V3*age
+covariates (3 plus signs)
+                          Tage[1=V3*age]= 4; Tage[2=age*V4] = 3
+                       */
+   if(decodemodel(model, lastobs) == 1)
+     goto end;
+   if((double)(lastobs-imx)/(double)imx > 1.10){
+     nbwarn++;
+     printf("Warning: The value of parameter lastobs=%d is big compared to the \n  effective number of cases imx=%d, please adjust, \n  otherwise you are allocating more memory than necessary.\n",lastobs, imx);
+     fprintf(ficlog,"Warning: The value of parameter lastobs=%d is big compared to the \n  effective number of cases imx=%d, please adjust, \n  otherwise you are allocating more memory than necessary.\n",lastobs, imx);
+   }
+     /*  if(mle==1){*/
+   if (weightopt != 1) { /* Maximisation without weights. We can have weights different from 1 but want no weight*/
+     for(i=1;i<=imx;i++) weight[i]=1.0; /* changed to imx */
+   }
+     /*-calculation of age at interview from date of interview and age at death -*/
+   agev=matrix(1,maxwav,1,imx);
+   if(calandcheckages(imx, maxwav, &agemin, &agemax, &nberr, &nbwarn) == 1)
+     goto end;
-   printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);
-   fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);
    agegomp=(int)agemin;
-   free_vector(severity,1,maxwav);
-   free_imatrix(outcome,1,maxwav+1,1,n);
    free_vector(moisnais,1,n);
    free_vector(annais,1,n);
    /* free_matrix(mint,1,maxwav,1,n);
       free_matrix(anint,1,maxwav,1,n);*/
    free_vector(moisdc,1,n);
    free_vector(andc,1,n);
+   /* */
    wav=ivector(1,imx);
    dh=imatrix(1,lastpass-firstpass+1,1,imx);
    bh=imatrix(1,lastpass-firstpass+1,1,imx);
- Line 5250  run imach with mle=-1 to get a correct t
+ Line 5936  run imach with mle=-1 to get a correct t
    /* Concatenates waves */
    concatwav(wav, dh, bh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);
+   /* */
    /* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */
    nbcode=imatrix(0,NCOVMAX,0,NCOVMAX);
    ncodemax[1]=1;
-   if (cptcovn > 0) tricode(Tvar,nbcode,imx);
+   Ndum =ivector(-1,NCOVMAX);
+   if (ncovmodel > 2)
-   codtab=imatrix(1,100,1,10); /* Cross tabulation to get the order of
+     tricode(Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */
-                                  the estimations*/
+   codtab=imatrix(1,100,1,10); /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */
+   /*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtab[100][10]);*/
    h=0;
+   /*if (cptcovn > 0) */
    m=pow(2,cptcoveff);
    for(k=1;k<=cptcoveff; k++){ /* scans any effective covariate */
-     for(i=1; i <=(m/pow(2,k));i++){ /* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 */
+     for(i=1; i <=pow(2,cptcoveff-k);i++){ /* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 */
-       for(j=1; j <= ncodemax[k]; j++){ /* For each modality of this covariate */
+       for(j=1; j <= ncodemax[k]; j++){ /* For each modality of this covariate ncodemax=2*/
-         for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){  /* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 */
+         for(cpt=1; cpt <=pow(2,k-1); cpt++){  /* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 */
            h++;
            if (h>m)
-             h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j;
+             h=1;
+           /**< codtab(h,k)  k   = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1
+            *     h     1     2     3     4
+            *______________________________
+            *     1 i=1 1 i=1 1 i=1 1 i=1 1
+            *     2     2     1     1     1
+            *     3 i=2 1     2     1     1
+            *     4     2     2     1     1
+            *     5 i=3 1 i=2 1     2     1
+            *     6     2     1     2     1
+            *     7 i=4 1     2     2     1
+            *     8     2     2     2     1
+            *     9 i=5 1 i=3 1 i=2 1     1
+            *    10     2     1     1     1
+            *    11 i=6 1     2     1     1
+            *    12     2     2     1     1
+            *    13 i=7 1 i=4 1     2     1
+            *    14     2     1     2     1
+            *    15 i=8 1     2     2     1
+            *    16     2     2     2     1
+            */
+           codtab[h][k]=j;
+           /*codtab[h][Tvar[k]]=j;*/
            printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][Tvar[k]]);
          }
        }
- Line 5283  run imach with mle=-1 to get a correct t
+ Line 5999  run imach with mle=-1 to get a correct t
       printf("\n");
       }
       scanf("%d",i);*/
+  free_ivector(Ndum,-1,NCOVMAX);
    /*------------ gnuplot -------------*/
    strcpy(optionfilegnuplot,optionfilefiname);
- Line 5296  run imach with mle=-1 to get a correct t
+ Line 6016  run imach with mle=-1 to get a correct t
    else{
      fprintf(ficgp,"\n# %s\n", version);
      fprintf(ficgp,"# %s\n", optionfilegnuplot);
-     fprintf(ficgp,"set missing 'NaNq'\n");
+     //fprintf(ficgp,"set missing 'NaNq'\n");
+     fprintf(ficgp,"set datafile missing 'NaNq'\n");
    }
    /*  fclose(ficgp);*/
    /*--------- index.htm --------*/
- Line 5366  Interval (in months) between two waves:
+ Line 6087  Interval (in months) between two waves:
    globpr=0; /* To get the number ipmx of contributions and the sum of weights*/
    if (mle==-3){
      ximort=matrix(1,NDIM,1,NDIM);
+ /*     ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */
      cens=ivector(1,n);
      ageexmed=vector(1,n);
      agecens=vector(1,n);
- Line 5404  Interval (in months) between two waves:
+ Line 6126  Interval (in months) between two waves:
          ximort[i][j]=(i == j ? 1.0 : 0.0);
      }
-     p[1]=0.0268; p[NDIM]=0.083;
+     /*p[1]=0.0268; p[NDIM]=0.083;*/
      /*printf("%lf %lf", p[1], p[2]);*/
+ #ifdef GSL
+     printf("GSL optimization\n");  fprintf(ficlog,"Powell\n");
+ #else
      printf("Powell\n");  fprintf(ficlog,"Powell\n");
+ #endif
      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);
      }
+ #ifdef GSL
+     fprintf(ficrespow,"# GSL optimization\n# iter -2*LL");
+ #else
      fprintf(ficrespow,"# Powell\n# iter -2*LL");
+ #endif
      /*  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");
+ #ifdef GSL
+     /* gsl starts here */
+     T = gsl_multimin_fminimizer_nmsimplex;
+     gsl_multimin_fminimizer *sfm = NULL;
+     gsl_vector *ss, *x;
+     gsl_multimin_function minex_func;
+     /* Initial vertex size vector */
+     ss = gsl_vector_alloc (NDIM);
+     if (ss == NULL){
+       GSL_ERROR_VAL ("failed to allocate space for ss", GSL_ENOMEM, 0);
+     }
+     /* Set all step sizes to 1 */
+     gsl_vector_set_all (ss, 0.001);
+     /* Starting point */
+     x = gsl_vector_alloc (NDIM);
+     if (x == NULL){
+       gsl_vector_free(ss);
+       GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0);
+     }
+     /* Initialize method and iterate */
+     /*     p[1]=0.0268; p[NDIM]=0.083; */
+ /*     gsl_vector_set(x, 0, 0.0268); */
+ /*     gsl_vector_set(x, 1, 0.083); */
+     gsl_vector_set(x, 0, p[1]);
+     gsl_vector_set(x, 1, p[2]);
+     minex_func.f = &gompertz_f;
+     minex_func.n = NDIM;
+     minex_func.params = (void *)&p; /* ??? */
+     sfm = gsl_multimin_fminimizer_alloc (T, NDIM);
+     gsl_multimin_fminimizer_set (sfm, &minex_func, x, ss);
+     printf("Iterations beginning .....\n\n");
+     printf("Iter. #    Intercept       Slope     -Log Likelihood     Simplex size\n");
+     iteri=0;
+     while (rval == GSL_CONTINUE){
+       iteri++;
+       status = gsl_multimin_fminimizer_iterate(sfm);
+       if (status) printf("error: %s\n", gsl_strerror (status));
+       fflush(0);
+       if (status)
+         break;
+       rval = gsl_multimin_test_size (gsl_multimin_fminimizer_size (sfm), 1e-6);
+       ssval = gsl_multimin_fminimizer_size (sfm);
+       if (rval == GSL_SUCCESS)
+         printf ("converged to a local maximum at\n");
+       printf("%5d ", iteri);
+       for (it = 0; it < NDIM; it++){
+         printf ("%10.5f ", gsl_vector_get (sfm->x, it));
+       }
+       printf("f() = %-10.5f ssize = %.7f\n", sfm->fval, ssval);
+     }
+     printf("\n\n Please note: Program should be run many times with varying starting points to detemine global maximum\n\n");
-     powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz);
+     gsl_vector_free(x); /* initial values */
+     gsl_vector_free(ss); /* inital step size */
+     for (it=0; it<NDIM; it++){
+       p[it+1]=gsl_vector_get(sfm->x,it);
+       fprintf(ficrespow," %.12lf", p[it]);
+     }
+     gsl_multimin_fminimizer_free (sfm); /* p *(sfm.x.data) et p *(sfm.x.data+1)  */
+ #endif
+ #ifdef POWELL
+      powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz);
+ #endif
      fclose(ficrespow);
      hesscov(matcov, p, NDIM, delti, 1e-4, gompertz);
- Line 5483  Interval (in months) between two waves:
+ Line 6290  Interval (in months) between two waves:
      free_vector(lsurv,1,AGESUP);
      free_vector(lpop,1,AGESUP);
      free_vector(tpop,1,AGESUP);
+ #ifdef GSL
+     free_ivector(cens,1,n);
+     free_vector(agecens,1,n);
+     free_ivector(dcwave,1,n);
+     free_matrix(ximort,1,NDIM,1,NDIM);
+ #endif
    } /* Endof if mle==-3 */
    else{ /* For mle >=1 */
- Line 5643  Interval (in months) between two waves:
+ Line 6456  Interval (in months) between two waves:
      while((c=getc(ficpar))=='#' && c!= EOF){
        ungetc(c,ficpar);
        fgets(line, MAXLINE, ficpar);
-       puts(line);
+       fputs(line,stdout);
        fputs(line,ficparo);
      }
      ungetc(c,ficpar);
- Line 5663  Interval (in months) between two waves:
+ Line 6476  Interval (in months) between two waves:
      while((c=getc(ficpar))=='#' && c!= EOF){
        ungetc(c,ficpar);
        fgets(line, MAXLINE, ficpar);
-       puts(line);
+       fputs(line,stdout);
        fputs(line,ficparo);
      }
      ungetc(c,ficpar);
- Line 5677  Interval (in months) between two waves:
+ Line 6490  Interval (in months) between two waves:
      while((c=getc(ficpar))=='#' && c!= EOF){
        ungetc(c,ficpar);
        fgets(line, MAXLINE, ficpar);
-       puts(line);
+       fputs(line,stdout);
        fputs(line,ficparo);
      }
      ungetc(c,ficpar);
- Line 5693  Interval (in months) between two waves:
+ Line 6506  Interval (in months) between two waves:
      while((c=getc(ficpar))=='#' && c!= EOF){
        ungetc(c,ficpar);
        fgets(line, MAXLINE, ficpar);
-       puts(line);
+       fputs(line,stdout);
        fputs(line,ficparo);
      }
      ungetc(c,ficpar);
- Line 5707  Interval (in months) between two waves:
+ Line 6520  Interval (in months) between two waves:
-     /*  freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint);*/
+      /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint); */
-     /*,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
+     /* ,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); */
      replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */
      printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);
- Line 5733  Interval (in months) between two waves:
+ Line 6546  Interval (in months) between two waves:
      /*--------------- Prevalence limit  (period or stable prevalence) --------------*/
+ #include "prevlim.h"  /* Use ficrespl, ficlog */
-     strcpy(filerespl,"pl");
-     strcat(filerespl,fileres);
-     if((ficrespl=fopen(filerespl,"w"))==NULL) {
-       printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);goto end;
-       fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);goto end;
-     }
-     printf("Computing period (stable) prevalence: result on file '%s' \n", filerespl);
-     fprintf(ficlog,"Computing period (stable) prevalence: result on file '%s' \n", filerespl);
-     pstamp(ficrespl);
-     fprintf(ficrespl,"# Period (stable) prevalence \n");
-     fprintf(ficrespl,"#Age ");
-     for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
-     fprintf(ficrespl,"\n");
-     prlim=matrix(1,nlstate,1,nlstate);
-     agebase=ageminpar;
-     agelim=agemaxpar;
-     ftolpl=1.e-10;
-     i1=cptcoveff;
-     if (cptcovn < 1){i1=1;}
-     for(cptcov=1,k=0;cptcov<=i1;cptcov++){
-       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
-         k=k+1;
-         /* to clean */
-         printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,codtab[cptcod][cptcov],nbcode);
-         fprintf(ficrespl,"\n#******");
-         printf("\n#******");
-         fprintf(ficlog,"\n#******");
-         for(j=1;j<=cptcoveff;j++) {
-           fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
-           printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
-           fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
-         }
-         fprintf(ficrespl,"******\n");
-         printf("******\n");
-         fprintf(ficlog,"******\n");
-         for (age=agebase; age<=agelim; age++){
-           prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
-           fprintf(ficrespl,"%.0f ",age );
-           for(j=1;j<=cptcoveff;j++)
-             fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
-           for(i=1; i<=nlstate;i++)
-             fprintf(ficrespl," %.5f", prlim[i][i]);
-           fprintf(ficrespl,"\n");
-         }
-       }
-     }
      fclose(ficrespl);
-     /*------------- h Pij x at various ages ------------*/
+ #ifdef FREEEXIT2
+ #include "freeexit2.h"
-     strcpy(filerespij,"pij");  strcat(filerespij,fileres);
+ #endif
-     if((ficrespij=fopen(filerespij,"w"))==NULL) {
-       printf("Problem with Pij resultfile: %s\n", filerespij);goto end;
-       fprintf(ficlog,"Problem with Pij resultfile: %s\n", filerespij);goto end;
-     }
-     printf("Computing pij: result on file '%s' \n", filerespij);
-     fprintf(ficlog,"Computing pij: result on file '%s' \n", filerespij);
-     stepsize=(int) (stepm+YEARM-1)/YEARM;
-     /*if (stepm<=24) stepsize=2;*/
-     agelim=AGESUP;
-     hstepm=stepsize*YEARM; /* Every year of age */
-     hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */
-     /* hstepm=1;   aff par mois*/
-     pstamp(ficrespij);
-     fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");
-     for(cptcov=1,k=0;cptcov<=i1;cptcov++){
-       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
-         k=k+1;
-         fprintf(ficrespij,"\n#****** ");
-         for(j=1;j<=cptcoveff;j++)
-           fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
-         fprintf(ficrespij,"******\n");
-         for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
-           nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
-           nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
-           /*      nhstepm=nhstepm*YEARM; aff par mois*/
-           p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+     /*------------- h Pij x at various ages ------------*/
-           oldm=oldms;savm=savms;
+ #include "hpijx.h"
-           hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);
+     fclose(ficrespij);
-           fprintf(ficrespij,"# Cov Agex agex+h hpijx with i,j=");
-           for(i=1; i<=nlstate;i++)
-             for(j=1; j<=nlstate+ndeath;j++)
-               fprintf(ficrespij," %1d-%1d",i,j);
-           fprintf(ficrespij,"\n");
-           for (h=0; h<=nhstepm; h++){
-             fprintf(ficrespij,"%d %3.f %3.f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm );
-             for(i=1; i<=nlstate;i++)
-               for(j=1; j<=nlstate+ndeath;j++)
-                 fprintf(ficrespij," %.5f", p3mat[i][j][h]);
-             fprintf(ficrespij,"\n");
-           }
-           free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
-           fprintf(ficrespij,"\n");
-         }
-       }
-     }
+   /*-------------- Variance of one-step probabilities---*/
+     k=1;
      varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart);
-     fclose(ficrespij);
      probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
      for(i=1;i<=AGESUP;i++)
- Line 5892  Interval (in months) between two waves:
+ Line 6609  Interval (in months) between two waves:
      }
      printf("Computing Health Expectancies: result on file '%s' \n", filerese);
      fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese);
-     for(cptcov=1,k=0;cptcov<=i1;cptcov++){
+     /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
-       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
+       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
-         k=k+1;
+     for (k=1; k <= (int) pow(2,cptcoveff); k++){
          fprintf(ficreseij,"\n#****** ");
          for(j=1;j<=cptcoveff;j++) {
            fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
- Line 5906  Interval (in months) between two waves:
+ Line 6624  Interval (in months) between two waves:
          evsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, strstart);
          free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
-       }
+       /*}*/
      }
      fclose(ficreseij);
- Line 5951  Interval (in months) between two waves:
+ Line 6669  Interval (in months) between two waves:
      printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
      fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
-     for(cptcov=1,k=0;cptcov<=i1;cptcov++){
+     /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
-       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
+       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
-         k=k+1;
-         fprintf(ficrest,"\n#****** ");
+     for (k=1; k <= (int) pow(2,cptcoveff); k++){
+         fprintf(ficrest,"\n#****** ");
          for(j=1;j<=cptcoveff;j++)
            fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
          fprintf(ficrest,"******\n");
- Line 5976  Interval (in months) between two waves:
+ Line 6695  Interval (in months) between two waves:
          eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
          oldm=oldms;savm=savms;
          cvevsij(eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);
+         /*
+          */
+         /* goto endfree; */
          vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
          pstamp(ficrest);
          for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
-           oldm=oldms;savm=savms;
+           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);   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 ");
+           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);
            else
- Line 6025  Interval (in months) between two waves:
+ Line 6751  Interval (in months) between two waves:
          free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
          free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
          free_vector(epj,1,nlstate+1);
-       }
+       /*}*/
      }
      free_vector(weight,1,n);
-     free_imatrix(Tvard,1,15,1,2);
+     free_imatrix(Tvard,1,NCOVMAX,1,2);
      free_imatrix(s,1,maxwav+1,1,n);
      free_matrix(anint,1,maxwav,1,n);
      free_matrix(mint,1,maxwav,1,n);
- Line 6050  Interval (in months) between two waves:
+ Line 6776  Interval (in months) between two waves:
      }
      printf("Computing Variance-covariance of period (stable) prevalence: file '%s' \n", fileresvpl);
-     for(cptcov=1,k=0;cptcov<=i1;cptcov++){
+     /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
-       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
+       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
-         k=k+1;
-         fprintf(ficresvpl,"\n#****** ");
+     for (k=1; k <= (int) pow(2,cptcoveff); k++){
+         fprintf(ficresvpl,"\n#****** ");
          for(j=1;j<=cptcoveff;j++)
            fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
          fprintf(ficresvpl,"******\n");
- Line 6062  Interval (in months) between two waves:
+ Line 6789  Interval (in months) between two waves:
          oldm=oldms;savm=savms;
          varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k,strstart);
          free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
-       }
+       /*}*/
      }
      fclose(ficresvpl);
- Line 6070  Interval (in months) between two waves:
+ Line 6797  Interval (in months) between two waves:
      /*---------- End : free ----------------*/
      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);
+     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);
      free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
- Line 6085  Interval (in months) between two waves:
+ Line 6811  Interval (in months) between two waves:
      free_matrix(agev,1,maxwav,1,imx);
      free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
-     free_ivector(ncodemax,1,8);
+     free_ivector(ncodemax,1,NCOVMAX);
-     free_ivector(Tvar,1,15);
+     free_ivector(Tvar,1,NCOVMAX);
-     free_ivector(Tprod,1,15);
+     free_ivector(Tprod,1,NCOVMAX);
-     free_ivector(Tvaraff,1,15);
+     free_ivector(Tvaraff,1,NCOVMAX);
-     free_ivector(Tage,1,15);
+     free_ivector(Tage,1,NCOVMAX);
      free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX);
      free_imatrix(codtab,1,100,1,10);
- Line 6106  Interval (in months) between two waves:
+ Line 6832  Interval (in months) between two waves:
    }
    printf("See log file on %s\n",filelog);
    /*  gettimeofday(&end_time, (struct timezone*)0);*/  /* after time */
-   (void) gettimeofday(&end_time,&tzp);
+   /*(void) gettimeofday(&end_time,&tzp);*/
-   tm = *localtime(&end_time.tv_sec);
+   rend_time = time(NULL);
-   tmg = *gmtime(&end_time.tv_sec);
+   end_time = *localtime(&rend_time);
-   strcpy(strtend,asctime(&tm));
+   /* tml = *localtime(&end_time.tm_sec); */
+   strcpy(strtend,asctime(&end_time));
    printf("Local time at start %s\nLocal time at end   %s",strstart, strtend);
    fprintf(ficlog,"Local time at start %s\nLocal time at end   %s\n",strstart, strtend);
-   printf("Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout));
+   printf("Total time used %s\n", asc_diff_time(rend_time -rstart_time,tmpout));
-   printf("Total time was %d Sec.\n", end_time.tv_sec -start_time.tv_sec);
+   printf("Total time was %.0lf Sec.\n", difftime(rend_time,rstart_time));
-   fprintf(ficlog,"Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout));
+   fprintf(ficlog,"Total time used %s\n", asc_diff_time(rend_time -rstart_time,tmpout));
-   fprintf(ficlog,"Total time was %d Sec.\n", end_time.tv_sec -start_time.tv_sec);
+   fprintf(ficlog,"Total time was %.0lf Sec.\n", difftime(rend_time,rstart_time));
    /*  printf("Total time was %d uSec.\n", total_usecs);*/
  /*   if(fileappend(fichtm,optionfilehtm)){ */
    fprintf(fichtm,"<br>Local time at start %s<br>Local time at end   %s<br>\n</body></html>",strstart, strtend);
- Line 6135  Interval (in months) between two waves:
+ Line 6862  Interval (in months) between two waves:
      printf("Current directory %s!\n",pathcd);
    /*strcat(plotcmd,CHARSEPARATOR);*/
    sprintf(plotcmd,"gnuplot");
- #ifndef UNIX
+ #ifdef _WIN32
    sprintf(plotcmd,"\"%sgnuplot.exe\"",pathimach);
  #endif
    if(!stat(plotcmd,&info)){
-     printf("Error gnuplot program not found: %s\n",plotcmd);fflush(stdout);
+     printf("Error or gnuplot program not found: '%s'\n",plotcmd);fflush(stdout);
      if(!stat(getenv("GNUPLOTBIN"),&info)){
-       printf("Error gnuplot program not found: %s Environment GNUPLOTBIN not set.\n",plotcmd);fflush(stdout);
+       printf("Error or gnuplot program not found: '%s' Environment GNUPLOTBIN not set.\n",plotcmd);fflush(stdout);
      }else
        strcpy(pplotcmd,plotcmd);
- #ifdef UNIX
+ #ifdef __unix
      strcpy(plotcmd,GNUPLOTPROGRAM);
      if(!stat(plotcmd,&info)){
-       printf("Error gnuplot program not found: %s\n",plotcmd);fflush(stdout);
+       printf("Error gnuplot program not found: '%s'\n",plotcmd);fflush(stdout);
      }else
        strcpy(pplotcmd,plotcmd);
  #endif
- Line 6155  Interval (in months) between two waves:
+ Line 6882  Interval (in months) between two waves:
      strcpy(pplotcmd,plotcmd);
    sprintf(plotcmd,"%s %s",pplotcmd, optionfilegnuplot);
-   printf("Starting graphs with: %s\n",plotcmd);fflush(stdout);
+   printf("Starting graphs with: '%s'\n",plotcmd);fflush(stdout);
    if((outcmd=system(plotcmd)) != 0){
-     printf("\n Problem with gnuplot\n");
+     printf("gnuplot command might not be in your path: '%s', err=%d\n", plotcmd, outcmd);
+     printf("\n Trying if gnuplot resides on the same directory that IMaCh\n");
+     sprintf(plotcmd,"%sgnuplot %s", pathimach, optionfilegnuplot);
+     if((outcmd=system(plotcmd)) != 0)
+       printf("\n Still a problem with gnuplot command %s, err=%d\n", plotcmd, outcmd);
    }
-   printf(" Wait...");
+   printf(" Successful, please wait...");
    while (z[0] != 'q') {
      /* chdir(path); */
-     printf("\nType e to edit output files, g to graph again and q for exiting: ");
+     printf("\nType e to edit results with your browser, g to graph again and q for exit: ");
      scanf("%s",z);
  /*     if (z[0] == 'c') system("./imach"); */
      if (z[0] == 'e') {
-       printf("Starting browser with: %s",optionfilehtm);fflush(stdout);
+ #ifdef __APPLE__
-       system(optionfilehtm);
+       sprintf(pplotcmd, "open %s", optionfilehtm);
+ #elif __linux
+       sprintf(pplotcmd, "xdg-open %s", optionfilehtm);
+ #else
+       sprintf(pplotcmd, "%s", optionfilehtm);
+ #endif
+       printf("Starting browser with: %s",pplotcmd);fflush(stdout);
+       system(pplotcmd);
      }
      else if (z[0] == 'g') system(plotcmd);
      else if (z[0] == 'q') exit(0);
- Line 6179  Interval (in months) between two waves:
+ Line 6917  Interval (in months) between two waves:
      scanf("%s",z);
    }
  }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>

Removed from v.1.135
changed lines
	Added in v.1.166