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.353 and 1.365

-version 1.353, 2023/05/08 18:48:22
+version 1.365, 2024/06/28 13:53:38
  Line 1
  /* $Id$
    $State$
    $Log$
+   Revision 1.365  2024/06/28 13:53:38  brouard
+   * imach.c (Module): fixing some bugs in gnuplot and quantitative variables, but not completely solved
+   Revision 1.364  2024/06/28 12:27:05  brouard
+   * imach.c (Module): fixing some bugs in gnuplot and quantitative variables, but not completely solved
+   Revision 1.363  2024/06/28 09:31:55  brouard
+   Summary: Adding log lines too
+   Revision 1.362  2024/06/28 08:00:31  brouard
+   Summary: 0.99s6
+   * imach.c (Module): s6 errors with age*age (harmless).
+   Revision 1.361  2024/05/12 20:29:32  brouard
+   Summary: Version 0.99s5
+   * src/imach.c Version 0.99s5 In fact, the covariance of total life
+   expectancy e.. with a partial life expectancy e.j is high,
+   therefore the complete matrix of variance covariance has to be
+   included in the formula of the standard error of the proportion of
+   total life expectancy spent in a specific state:
+   var(X/Y)=mu_x^2/mu_y^2*(sigma_x^2/mu_x^2 -2
+   sigma_xy/mu_x/mu_y+sigma^2/mu_y^2).  Also an error with mle=-3
+   made the program core dump. It is fixed in this version.
+   Revision 1.360  2024/04/30 10:59:22  brouard
+   Summary: Version 0.99s4 and estimation of std of e.j/e..
+   Revision 1.359  2024/04/24 21:21:17  brouard
+   Summary: First IMaCh version using Brent Praxis software based on Buckhardt and Gegenfürtner C codes
+   Revision 1.6  2024/04/24 21:10:29  brouard
+   Summary: First IMaCh version using Brent Praxis software based on Buckhardt and Gegenfürtner C codes
+   Revision 1.5  2023/10/09 09:10:01  brouard
+   Summary: trying to reconsider
+   Revision 1.4  2023/06/22 12:50:51  brouard
+   Summary: stil on going
+   Revision 1.3  2023/06/22 11:28:07  brouard
+   *** empty log message ***
+   Revision 1.2  2023/06/22 11:22:40  brouard
+   Summary: with svd but not working yet
    Revision 1.353  2023/05/08 18:48:22  brouard
    *** empty log message ***
- Line 1207  Important routines
+ Line 1254  Important routines
  - Tricode which tests the modality of dummy variables (in order to warn with wrong or empty modalities)
    and returns the number of efficient covariates cptcoveff and modalities nbcode[Tvar[k]][1]= 0 and nbcode[Tvar[k]][2]= 1 usually.
  - printinghtml which outputs results like life expectancy in and from a state for a combination of modalities of dummy variables
-   o There are 2**cptcoveff combinations of (0,1) for cptcoveff variables. Outputting only combinations with people, éliminating 1 1 if
+   o There are 2**cptcoveff combinations of (0,1) for cptcoveff variables. Outputting only combinations with people, eliminating 1 1 if
      race White (0 0), Black vs White (1 0), Hispanic (0 1) and 1 1 being meaningless.
- Line 1278  Important routines
+ Line 1325  Important routines
  /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */
  /* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */
  /* #define FLATSUP  *//* Suppresses directions where likelihood is flat */
+ /* #define POWELLORIGINCONJUGATE  /\* Don't use conjugate but biggest decrease if valuable *\/ */
+ /* #define NOTMINFIT */
  #include <math.h>
  #include <stdio.h>
- Line 1374  double gnuplotversion=GNUPLOTVERSION;
+ Line 1423  double gnuplotversion=GNUPLOTVERSION;
  /* $State$ */
  #include "version.h"
  char version[]=__IMACH_VERSION__;
- char copyright[]="April 2023,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2022";
+ char copyright[]="April 2024,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2024";
  char fullversion[]="$Revision$ $Date$";
  char strstart[80];
  char optionfilext[10], optionfilefiname[FILENAMELENGTH];
- Line 1422  int *wav; /* Number of waves for this in
+ Line 1471  int *wav; /* Number of waves for this in
  int maxwav=0; /* Maxim number of waves */
  int jmin=0, jmax=0; /* min, max spacing between 2 waves */
  int ijmin=0, ijmax=0; /* Individuals having jmin and jmax */
- int gipmx=0, gsw=0; /* Global variables on the number of contributions
+ int gipmx = 0;
+ double gsw = 0; /* Global variables on the number of contributions
                     to the likelihood and the sum of weights (done by funcone)*/
  int mle=1, weightopt=0;
  int **mw; /* mw[mi][i] is number of the mi wave for this individual */
- Line 1495  char *endptr;
+ Line 1545  char *endptr;
  long lval;
  double dval;
+ /* This for praxis gegen */
+   /* int prin=1; */
+   double h0=0.25;
+   double macheps;
+   double ffmin;
  #define NR_END 1
  #define FREE_ARG char*
  #define FTOL 1.0e-10
- Line 1612  int **nbcode, *Tvar; /**< model=V2 => Tv
+ Line 1668  int **nbcode, *Tvar; /**< model=V2 => Tv
  /* Tage[cptcovage]=k            5               8         10 */ /* Position in the model of ith cov*age */
  /* model="V2+V3+V4+V6+V7+V6*V2+V7*V2+V6*V3+V7*V3+V6*V4+V7*V4+age*V2+age*V3+age*V4+age*V6+age*V7+age*V6*V2+age*V6*V3+age*V7*V3+age*V6*V4+age*V7*V4\r"*/
  /*  p Tvard[1][1]@21 = {6, 2, 7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0}*/
- /*  p Tvard[2][1]@21 = {7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0 <repeats 11 times>}
+ /*  p Tvard[2][1]@21 = {7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0 <repeats 11 times>} */
  /* p Tvardk[1][1]@24 = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 2, 7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0, 0}*/
  /* p Tvardk[1][1]@22 = {0, 0, 0, 0, 0, 0, 0, 0, 6, 2, 7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0, 0} */
  /* Tvard[1][1]@4={4,3,1,2}    V4*V3 V1*V2               */ /* Position in model of the ith prod without age */
- Line 2609  void linmin(double p[], double xi[], int
+ Line 2665  void linmin(double p[], double xi[], int
    free_vector(pcom,1,n);
  }
+ /**** praxis gegen ****/
+ /* This has been tested by Visual C from Microsoft and works */
+ /* meaning tha valgrind could be wrong */
+ /*********************************************************************/
+ /*      f u n c t i o n     p r a x i s                              */
+ /*                                                                   */
+ /* praxis is a general purpose routine for the minimization of a     */
+ /* function in several variables. the algorithm used is a modifi-    */
+ /* cation of conjugate gradient search method by powell. the changes */
+ /* are due to r.p. brent, who gives an algol-w program, which served */
+ /* as a basis for this function.                                     */
+ /*                                                                   */
+ /* references:                                                       */
+ /*     - powell, m.j.d., 1964. an efficient method for finding       */
+ /*       the minimum of a function in several variables without      */
+ /*       calculating derivatives, computer journal, 7, 155-162       */
+ /*     - brent, r.p., 1973. algorithms for minimization without      */
+ /*       derivatives, prentice hall, englewood cliffs.               */
+ /*                                                                   */
+ /*     problems, suggestions or improvements are always wellcome     */
+ /*                       karl gegenfurtner   07/08/87                */
+ /*                                           c - version             */
+ /*********************************************************************/
+ /*                                                                   */
+ /* usage: min = praxis(tol, macheps, h, n, prin, x, func)      */
+ /* macheps has been suppressed because it is replaced by DBL_EPSILON */
+ /* and if it was an argument of praxis (as it is in original brent)  */
+ /* it should be declared external */
+ /* usage: min = praxis(tol, h, n, prin, x, func)      */
+ /* was    min = praxis(fun, x, n);                                   */
+ /*                                                                   */
+ /*  fun        the function to be minimized. fun is called from      */
+ /*             praxis with x and n as arguments                      */
+ /*  x          a double array containing the initial guesses for     */
+ /*             the minimum, which will contain the solution on       */
+ /*             return                                                */
+ /*  n          an integer specifying the number of unknown           */
+ /*             parameters                                            */
+ /*  min        praxis returns the least calculated value of fun      */
+ /*                                                                   */
+ /* some additional global variables control some more aspects of     */
+ /* the inner workings of praxis. setting them is optional, they      */
+ /* are all set to some reasonable default values given below.        */
+ /*                                                                   */
+ /*   prin      controls the printed output from the routine.         */
+ /*             0 -> no output                                        */
+ /*             1 -> print only starting and final values             */
+ /*             2 -> detailed map of the minimization process         */
+ /*             3 -> print also eigenvalues and vectors of the        */
+ /*                  search directions                                */
+ /*             the default value is 1                                */
+ /*  tol        is the tolerance allowed for the precision of the     */
+ /*             solution. praxis returns if the criterion             */
+ /*             2 * ||x[k]-x[k-1]|| <= sqrt(macheps) * ||x[k]|| + tol */
+ /*             is fulfilled more than ktm times.                     */
+ /*             the default value depends on the machine precision    */
+ /*  ktm        see just above. default is 1, and a value of 4 leads  */
+ /*             to a very(!) cautious stopping criterion.             */
+ /*  h0 or step       is a steplength parameter and should be set equal     */
+ /*             to the expected distance from the solution.           */
+ /*             exceptionally small or large values of step lead to   */
+ /*             slower convergence on the first few iterations        */
+ /*             the default value for step is 1.0                     */
+ /*  scbd       is a scaling parameter. 1.0 is the default and        */
+ /*             indicates no scaling. if the scales for the different */
+ /*             parameters are very different, scbd should be set to  */
+ /*             a value of about 10.0.                                */
+ /*  illc       should be set to true (1) if the problem is known to  */
+ /*             be ill-conditioned. the default is false (0). this    */
+ /*             variable is automatically set, when praxis finds      */
+ /*             the problem to be ill-conditioned during iterations.  */
+ /*  maxfun     is the maximum number of calls to fun allowed. praxis */
+ /*             will return after maxfun calls to fun even when the   */
+ /*             minimum is not yet found. the default value of 0      */
+ /*             indicates no limit on the number of calls.            */
+ /*             this return condition is only checked every n         */
+ /*             iterations.                                           */
+ /*                                                                   */
+ /*********************************************************************/
+ #include <math.h>
+ #include <stdio.h>
+ #include <stdlib.h>
+ #include <float.h> /* for DBL_EPSILON */
+ /* #include "machine.h" */
+ /* extern void minfit(int n, double eps, double tol, double **ab, double q[]); */
+ /* extern void minfit(int n, double eps, double tol, double ab[N][N], double q[]); */
+ /* control parameters */
+ /* control parameters */
+ #define SQREPSILON 1.0e-19
+ /* #define EPSILON 1.0e-8 */ /* in main */
+ double tol = SQREPSILON,
+        scbd = 1.0,
+        step = 1.0;
+ int    ktm = 1,
+        /* prin = 2, */
+        maxfun = 0,
+        illc = 0;
+ /* some global variables */
+ static int i, j, k, k2, nl, nf, kl, kt;
+ /* static double s; */
+ double sl, dn, dmin,
+        fx, f1, lds, ldt, sf, df,
+        qf1, qd0, qd1, qa, qb, qc,
+        m2, m4, small_windows, vsmall, large,
+        vlarge, ldfac, t2;
+ /* static double d[N], y[N], z[N], */
+ /*        q0[N], q1[N], v[N][N]; */
+ static double *d, *y, *z;
+ static double  *q0, *q1, **v;
+ double *tflin; /* used in flin: return (*fun)(tflin, n); */
+ double *e; /* used in minfit, don't konw how to free memory and thus made global */
+ /* static double s, sl, dn, dmin, */
+ /*        fx, f1, lds, ldt, sf, df, */
+ /*        qf1, qd0, qd1, qa, qb, qc, */
+ /*        m2, m4, small, vsmall, large,  */
+ /*        vlarge, ldfac, t2; */
+ /* static double d[N], y[N], z[N], */
+ /*        q0[N], q1[N], v[N][N]; */
+ /* these will be set by praxis to point to it's arguments */
+ static int prin; /* added */
+ static int n;
+ static double *x;
+ static double (*fun)();
+ /* static double (*fun)(double *x, int n); */
+ /* these will be set by praxis to the global control parameters */
+ /* static double h, macheps, t; */
+ extern double macheps;
+ static double h;
+ static double t;
+ static double
+ drandom()       /* return random no between 0 and 1 */
+ {
+    return (double)(rand()%(8192*2))/(double)(8192*2);
+ }
+ static void sort()              /* d and v in descending order */
+ {
+    int k, i, j;
+    double s;
+    for (i=1; i<=n-1; i++) {
+        k = i; s = d[i];
+        for (j=i+1; j<=n; j++) {
+            if (d[j] > s) {
+               k = j;
+               s = d[j];
+            }
+        }
+        if (k > i) {
+           d[k] = d[i];
+           d[i] = s;
+           for (j=1; j<=n; j++) {
+               s = v[j][i];
+               v[j][i] = v[j][k];
+               v[j][k] = s;
+           }
+        }
+    }
+ }
+ double randbrent ( int *naught )
+ {
+   double ran1, ran3[127], half;
+   int ran2, q, r, i, j;
+   int init=0; /* false */
+   double rr;
+   /* REAL*8 RAN1,RAN3(127),HALF */
+   /*     INTEGER RAN2,Q,R */
+   /*     LOGICAL INIT */
+   /*     DATA INIT/.FALSE./ */
+   /*     IF (INIT) GO TO 3 */
+   if(!init){
+ /*       R = MOD(NAUGHT,8190) + 1 *//* 1804289383 rand () */
+     r = *naught % 8190 + 1;/* printf(" naught r %d %d",*naught,r); */
+     ran2=127;
+     for(i=ran2; i>0; i--){
+ /*       RAN2 = 128 */
+ /*       DO 2 I=1,127 */
+       ran2 = ran2-1;
+ /*          RAN2 = RAN2 - 1 */
+       ran1 = -pow(2.0,55);
+ /*          RAN1 = -2.D0**55 */
+ /*          DO 1 J=1,7 */
+       for(j=1; j<=7;j++){
+ /*             R = MOD(1756*R,8191) */
+         r = (1756*r) % 8191;/* printf(" i=%d (1756*r)%8191=%d",j,r); */
+         q=r/32;
+ /*             Q = R/32 */
+ /* 1           RAN1 = (RAN1 + Q)*(1.0D0/256) */
+         ran1 =(ran1+q)*(1.0/256);
+       }
+ /* 2        RAN3(RAN2) = RAN1 */
+       ran3[ran2] = ran1; /* printf(" ran2=%d ran1=%.7g \n",ran2,ran1); */
+     }
+ /*       INIT = .TRUE. */
+     init=1;
+ /* 3     IF (RAN2.EQ.1) RAN2 = 128 */
+   }
+   if(ran2 == 0) ran2 = 126;
+   else ran2 = ran2 -1;
+   /* RAN2 = RAN2 - 1 */
+   /* RAN1 = RAN1 + RAN3(RAN2) */
+   ran1 = ran1 + ran3[ran2];/* printf("BIS ran2=%d ran1=%.7g \n",ran2,ran1);  */
+   half= 0.5;
+   /* HALF = .5D0 */
+   /* IF (RAN1.GE.0.D0) HALF = -HALF */
+   if(ran1 >= 0.) half =-half;
+   ran1 = ran1 +half;
+   ran3[ran2] = ran1;
+   rr= ran1+0.5;
+   /* RAN1 = RAN1 + HALF */
+   /*   RAN3(RAN2) = RAN1 */
+   /*   RANDOM = RAN1 + .5D0 */
+ /*   r = ( ( double ) ( *seed ) ) * 4.656612875E-10; */
+   return rr;
+ }
+ static void matprint(char *s, double **v, int m, int n)
+ /* char *s; */
+ /* double v[N][N]; */
+ {
+ #define INCX 8
+   int i;
+   int i2hi;
+   int ihi;
+   int ilo;
+   int i2lo;
+   int jlo=1;
+   int j;
+   int j2hi;
+   int jhi;
+   int j2lo;
+   ilo=1;
+   ihi=n;
+   jlo=1;
+   jhi=n;
+   printf ("\n" );
+   printf ("%s\n", s );
+   for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX )
+   {
+     j2hi = j2lo + INCX - 1;
+     if ( n < j2hi )
+     {
+       j2hi = n;
+     }
+     if ( jhi < j2hi )
+     {
+       j2hi = jhi;
+     }
+     /* fprintf ( ficlog, "\n" ); */
+     printf ("\n" );
+ /*
+   For each column J in the current range...
+   Write the header.
+ */
+     /* fprintf ( ficlog, "  Col:  "); */
+     printf ("Col:");
+     for ( j = j2lo; j <= j2hi; j++ )
+     {
+       /* fprintf ( ficlog, "  %7d     ", j - 1 ); */
+       /* printf (" %9d      ", j - 1 ); */
+       printf (" %9d      ", j );
+     }
+     /* fprintf ( ficlog, "\n" ); */
+     /* fprintf ( ficlog, "  Row\n" ); */
+     /* fprintf ( ficlog, "\n" ); */
+     printf ("\n" );
+     printf ("  Row\n" );
+     printf ("\n" );
+ /*
+   Determine the range of the rows in this strip.
+ */
+     if ( 1 < ilo ){
+       i2lo = ilo;
+     }else{
+       i2lo = 1;
+     }
+     if ( m < ihi ){
+       i2hi = m;
+     }else{
+       i2hi = ihi;
+     }
+     for ( i = i2lo; i <= i2hi; i++ ){
+ /*
+   Print out (up to) 5 entries in row I, that lie in the current strip.
+ */
+       /* fprintf ( ficlog, "%5d:", i - 1 ); */
+       /* printf ("%5d:", i - 1 ); */
+       printf ("%5d:", i );
+       for ( j = j2lo; j <= j2hi; j++ )
+       {
+         /* fprintf ( ficlog, "  %14g", a[i-1+(j-1)*m] ); */
+         /* printf ("%14.7g  ", a[i-1+(j-1)*m] ); */
+            /* printf("%14.7f  ", v[i-1][j-1]); */
+            printf("%14.7f  ", v[i][j]);
+         /* fprintf ( stdout, "  %14g", a[i-1+(j-1)*m] ); */
+       }
+       /* fprintf ( ficlog, "\n" ); */
+       printf ("\n" );
+     }
+   }
+    /* printf("%s\n", s); */
+    /* for (k=0; k<n; k++) { */
+    /*     for (i=0; i<n; i++) { */
+    /*         /\* printf("%20.10e ", v[k][i]); *\/ */
+    /*     } */
+    /*     printf("\n"); */
+    /* } */
+ #undef INCX
+ }
+ void vecprint(char *s, double *x, int n)
+ /* char *s; */
+ /* double x[N]; */
+ {
+    int i=0;
+    printf(" %s", s);
+    /* for (i=0; i<n; i++) */
+    for (i=1; i<=n; i++)
+      printf ("  %14.7g",  x[i] );
+      /* printf("  %8d: %14g\n", i, x[i]); */
+    printf ("\n" );
+ }
+ static void print()             /* print a line of traces */
+ {
+    printf("\n");
+    /* printf("... chi square reduced to ... %20.10e\n", fx); */
+    /* printf("... after %u function calls ...\n", nf); */
+    /* printf("... including %u linear searches ...\n", nl); */
+    printf("%10d    %10d%14.7g",nl, nf, fx);
+    vecprint("... current values of x ...", x, n);
+ }
+ /* static void print2(int n, double *x, int prin, double fx, int nf, int nl) */ /* print a line of traces */
+ static void print2() /* print a line of traces */
+ {
+   int i; double fmin=0.;
+    /* printf("\n"); */
+    /* printf("... chi square reduced to ... %20.10e\n", fx); */
+    /* printf("... after %u function calls ...\n", nf); */
+    /* printf("... including %u linear searches ...\n", nl); */
+    /* printf("%10d    %10d%14.7g",nl, nf, fx); */
+   /* printf ( "\n" ); */
+   printf ( "  Linear searches      %d", nl );
+   fprintf (ficlog, "  Linear searches      %d", nl );
+   /* printf ( "  Linear searches      %d\n", nl ); */
+   /* printf ( "  Function evaluations %d\n", nf ); */
+   /* printf ( "  Function value FX = %g\n", fx ); */
+   printf ( "  Function evaluations %d", nf );
+   printf ( "  Function value FX = %.12lf\n", fx );
+   fprintf (ficlog, "  Function evaluations %d", nf );
+   fprintf (ficlog, "  Function value FX = %.12lf\n", fx );
+ #ifdef DEBUGPRAX
+    printf("n=%d prin=%d\n",n,prin);
+ #endif
+    /* if(fx <= fmin) printf(" UNDEFINED "); else  printf("%14.7g",log(fx-fmin)); */
+    if ( n <= 4 || 2 < prin )
+    {
+      /* for(i=1;i<=n;i++)printf("%14.7g",x[i-1]); */
+      for(i=1;i<=n;i++){
+        printf(" %14.7g",x[i]);
+        fprintf(ficlog," %14.7g",x[i]);
+      }
+      /* r8vec_print ( n, x, "  X:" ); */
+    }
+    printf("\n");
+    fprintf(ficlog,"\n");
+  }
+ /* #ifdef MSDOS */
+ /* static double tflin[N]; */
+ /* #endif */
+ static double flin(double l, int j)
+ /* double l; */
+ {
+    int i;
+    /* #ifndef MSDOS */
+    /*    double tflin[N]; */
+    /* #endif    */
+    /* double *tflin; */ /* Be careful to put tflin on a vector n */
+    /* j is used from 0 to n-1 and can be -1 for parabolic search */
+    /* if (j != -1) {            /\* linear search *\/ */
+    if (j > 0) {         /* linear search */
+      /* for (i=0; i<n; i++){ */
+      for (i=1; i<=n; i++){
+           tflin[i] = x[i] + l *v[i][j];
+ #ifdef DEBUGPRAX
+           /* printf("     flin i=%14d t=%14.7f x=%14.7f l=%14.7f v[%d,%d]=%14.7f nf=%14d\n",i+1, tflin[i],x[i],l,i,j,v[i][j],nf); */
+           printf("     flin i=%14d t=%14.7f x=%14.7f l=%14.7f v[%d,%d]=%14.7f nf=%14d\n",i, tflin[i],x[i],l,i,j,v[i][j],nf);
+ #endif
+      }
+    }
+    else {                       /* search along parabolic space curve */
+       qa = l*(l-qd1)/(qd0*(qd0+qd1));
+       qb = (l+qd0)*(qd1-l)/(qd0*qd1);
+       qc = l*(l+qd0)/(qd1*(qd0+qd1));
+ #ifdef DEBUGPRAX
+       printf("     search along a parabolic space curve. j=%14d nf=%14d l=%14.7f qd0=%14.7f qd1=%14.7f\n",j,nf,l,qd0,qd1);
+ #endif
+       /* for (i=0; i<n; i++){ */
+       for (i=1; i<=n; i++){
+           tflin[i] = qa*q0[i]+qb*x[i]+qc*q1[i];
+ #ifdef DEBUGPRAX
+           /* printf("      parabole i=%14d t(i)=%14.7f q0=%14.7f x=%14.7f q1=%14.7f\n",i+1,tflin[i],q0[i],x[i],q1[i]); */
+           printf("      parabole i=%14d t(i)=%14.7e q0=%14.7e x=%14.7e q1=%14.7e\n",i,tflin[i],q0[i],x[i],q1[i]);
+ #endif
+       }
+    }
+    nf++;
+ #ifdef NR_SHIFT
+       return (*fun)((tflin-1), n);
+ #else
+      /* return (*fun)(tflin, n);*/
+       return (*fun)(tflin);
+ #endif
+ }
+ void minny(int j, int nits, double *d2, double *x1, double f1, int fk)
+ /* double *d2, *x1, f1; */
+ {
+ /* here j is from 0 to n-1 and can be -1 for parabolic search  */
+   /*      MINIMIZES F FROM X IN THE DIRECTION V(*,J) */
+           /*      UNLESS J<1, WHEN A QUADRATIC SEARCH IS DONE */
+           /*      IN THE PLANE DEFINED BY Q0, Q1 AND X. */
+           /*      D2 AN APPROXIMATION TO HALF F'' (OR ZERO), */
+           /*      X1 AN ESTIMATE OF DISTANCE TO MINIMUM, */
+           /*      RETURNED AS THE DISTANCE FOUND. */
+           /*       IF FK = TRUE THEN F1 IS FLIN(X1), OTHERWISE */
+           /*       X1 AND F1 ARE IGNORED ON ENTRY UNLESS FINAL */
+           /*       FX > F1. NITS CONTROLS THE NUMBER OF TIMES */
+           /*       AN ATTEMPT IS MADE TO HALVE THE INTERVAL. */
+           /* SIDE EFFECTS: USES AND ALTERS X, FX, NF, NL. */
+           /*       IF J < 1 USES VARIABLES Q... . */
+           /*       USES H, N, T, M2, M4, LDT, DMIN, MACHEPS; */
+    int k, i, dz;
+    double x2, xm, f0, f2, fm, d1, t2, sf1, sx1;
+    double s;
+    double macheps;
+    macheps=pow(16.0,-13.0);
+    sf1 = f1; sx1 = *x1;
+    k = 0; xm = 0.0; fm = f0 = fx; dz = *d2 < macheps;
+    /* h=1.0;*/ /* To be revised */
+ #ifdef DEBUGPRAX
+    /* printf("min macheps=%14g h=%14g step=%14g t=%14g fx=%14g\n",macheps,h, step,t, fx);  */
+    /* Where is fx coming from */
+    printf("   min macheps=%14g h=%14g  t=%14g fx=%.9lf dirj=%d\n",macheps, h, t, fx, j);
+    matprint("  min vectors:",v,n,n);
+ #endif
+    /* find step size */
+    s = 0.;
+    /* for (i=0; i<n; i++) s += x[i]*x[i]; */
+    for (i=1; i<=n; i++) s += x[i]*x[i];
+    s = sqrt(s);
+    if (dz)
+       t2 = m4*sqrt(fabs(fx)/dmin + s*ldt) + m2*ldt;
+    else
+       t2 = m4*sqrt(fabs(fx)/(*d2) + s*ldt) + m2*ldt;
+    s = s*m4 + t;
+    if (dz && t2 > s) t2 = s;
+    if (t2 < small_windows) t2 = small_windows;
+    if (t2 > 0.01*h) t2 = 0.01 * h;
+    if (fk && f1 <= fm) {
+       xm = *x1;
+       fm = f1;
+    }
+ #ifdef DEBUGPRAX
+    printf("   additional flin X1=%14.7f t2=%14.7f *f1=%14.7f fm=%14.7f fk=%d\n",*x1,t2,f1,fm,fk);
+ #endif
+    if (!fk || fabs(*x1) < t2) {
+      *x1 = (*x1 >= 0 ? t2 : -t2);
+       /* *x1 = (*x1 > 0 ? t2 : -t2); */ /* kind of error */
+ #ifdef DEBUGPRAX
+      printf("    additional flin X1=%16.10e dirj=%d fk=%d\n",*x1, j, fk);
+ #endif
+       f1 = flin(*x1, j);
+ #ifdef DEBUGPRAX
+     printf("    after flin f1=%18.12e dirj=%d fk=%d\n",f1, j,fk);
+ #endif
+    }
+    if (f1 <= fm) {
+       xm = *x1;
+       fm = f1;
+    }
+ L0: /*L0 loop or next */
+ /*
+   Evaluate FLIN at another point and estimate the second derivative.
+ */
+    if (dz) {
+       x2 = (f0 < f1 ? -(*x1) : 2*(*x1));
+ #ifdef DEBUGPRAX
+       printf("     additional second flin x2=%14.8e x1=%14.8e f0=%14.8e f1=%18.12e dirj=%d\n",x2,*x1,f0,f1,j);
+ #endif
+       f2 = flin(x2, j);
+ #ifdef DEBUGPRAX
+       printf("     additional second flin x2=%16.10e x1=%16.10e f1=%18.12e f0=%18.10e f2=%18.10e fm=%18.10e\n",x2, *x1, f1,f0,f2,fm);
+ #endif
+       if (f2 <= fm) {
+          xm = x2;
+          fm = f2;
+       }
+       /* d2 is the curvature or double difference f1 doesn't seem to be accurately computed */
+       *d2 = (x2*(f1-f0) - (*x1)*(f2-f0))/((*x1)*x2*((*x1)-x2));
+ #ifdef DEBUGPRAX
+       double d11,d12;
+       d11=(f1-f0)/(*x1);d12=(f2-f0)/x2;
+       printf(" d11=%18.12e d12=%18.12e d11-d12=%18.12e x1-x2=%18.12e (d11-d12)/(x2-(*x1))=%18.12e\n", d11 ,d12, d11-d12, x2-(*x1), (d11-d12)/(x2-(*x1)));
+       printf(" original computing f1=%18.12e *d2=%16.10e f0=%18.12e f1-f0=%16.10e f2-f0=%16.10e\n",f1,*d2,f0,f1-f0, f2-f0);
+       double ff1=7.783920622852e+04;
+       double f1mf0=9.0344736236e-05;
+       *d2 = (f1mf0)/ (*x1)/((*x1)-x2) - (f2-f0)/x2/((*x1)-x2);
+       /* *d2 = (ff1-f0)/ (*x1)/((*x1)-x2) - (f2-f0)/x2/((*x1)-x2); */
+       printf(" simpliff computing *d2=%16.10e f1mf0=%18.12e,f1=f0+f1mf0=%18.12e\n",*d2,f1mf0,f0+f1mf0);
+       *d2 = ((f1-f0)/ (*x1) - (f2-f0)/x2)/((*x1)-x2);
+       printf(" overlifi computing *d2=%16.10e\n",*d2);
+ #endif
+       *d2 = ((f1-f0)/ (*x1) - (f2-f0)/x2)/((*x1)-x2);
+    }
+ #ifdef DEBUGPRAX
+       printf("    additional second flin xm=%14.8e fm=%14.8e *d2=%14.8e\n",xm, fm,*d2);
+ #endif
+    /*
+      Estimate the first derivative at 0.
+    */
+    d1 = (f1-f0)/(*x1) - *x1**d2; dz = 1;
+    /*
+       Predict the minimum.
+     */
+    if (*d2 <= small_windows) {
+      x2 = (d1 < 0 ? h : -h);
+    }
+    else {
+       x2 = - 0.5*d1/(*d2);
+    }
+ #ifdef DEBUGPRAX
+     printf("   AT d1=%14.8e d2=%14.8e small=%14.8e dz=%d x1=%14.8e x2=%14.8e\n",d1,*d2,small_windows,dz,*x1,x2);
+ #endif
+     if (fabs(x2) > h)
+       x2 = (x2 > 0 ? h : -h);
+ L1:  /* L1 or try loop */
+ #ifdef DEBUGPRAX
+     printf("   AT predicted minimum flin x2=%14.8e x1=%14.8e K=%14d NITS=%14d dirj=%d\n",x2,*x1,k,nits,j);
+ #endif
+    f2 = flin(x2, j); /* x[i]+x2*v[i][j] */
+ #ifdef DEBUGPRAX
+    printf("   after flin f0=%14.8e f1=%14.8e f2=%14.8e fm=%14.8e\n",f0,f1,f2, fm);
+ #endif
+    if ((k < nits) && (f2 > f0)) {
+ #ifdef DEBUGPRAX
+      printf("  NO SUCCESS SO TRY AGAIN;\n");
+ #endif
+      k++;
+      if ((f0 < f1) && (*x1*x2 > 0.0))
+        goto L0; /* or next */
+      x2 *= 0.5;
+      goto L1;
+    }
+    nl++;
+ #ifdef DEBUGPRAX
+    printf(" bebeBE end of min x1=%14.8e x2=%14.8e f1=%14.8e f2=%14.8e f0=%14.8e fm=%14.8e d2=%14.8e\n",*x1, x2, f1, f2, f0, fm, *d2);
+ #endif
+    if (f2 > fm) x2 = xm; else fm = f2;
+    if (fabs(x2*(x2-*x1)) > small_windows) {
+       *d2 = (x2*(f1-f0) - *x1*(fm-f0))/(*x1*x2*(*x1-x2));
+    }
+    else {
+       if (k > 0) *d2 = 0;
+    }
+ #ifdef DEBUGPRAX
+    printf(" bebe end of min x1 might be very wrong x1=%14.8e fx=%14.8e d2=%14.8e\n",*x1, fx, *d2);
+ #endif
+    if (*d2 <= small_windows) *d2 = small_windows;
+    *x1 = x2; fx = fm;
+    if (sf1 < fx) {
+       fx = sf1;
+       *x1 = sx1;
+    }
+   /*
+     Update X for linear search.
+   */
+ #ifdef DEBUGPRAX
+    printf("  end of min x1=%14.8e fx=%14.8e d2=%14.8e\n",*x1, fx, *d2);
+ #endif
+    /* if (j != -1) */
+    /*    for (i=0; i<n; i++) */
+    /*        x[i] += (*x1)*v[i][j]; */
+    if (j > 0)
+       for (i=1; i<=n; i++)
+           x[i] += (*x1)*v[i][j];
+ }
+ void quad()     /* look for a minimum along the curve q0, q1, q2        */
+ {
+    int i;
+    double l, s;
+    s = fx; fx = qf1; qf1 = s; qd1 = 0.0;
+    /* for (i=0; i<n; i++) { */
+    for (i=1; i<=n; i++) {
+        s = x[i]; l = q1[i]; x[i] = l; q1[i] = s;
+        qd1 = qd1 + (s-l)*(s-l);
+    }
+    s = 0.0; qd1 = sqrt(qd1); l = qd1;
+ #ifdef DEBUGPRAX
+   printf("  QUAD after sqrt qd1=%14.8e \n",qd1);
+ #endif
+    if (qd0>0.0 && qd1>0.0 &&nl>=3*n*n) {
+ #ifdef DEBUGPRAX
+      printf(" QUAD before min value=%14.8e \n",qf1);
+ #endif
+       /* min(-1, 2, &s, &l, qf1, 1); */
+       minny(0, 2, &s, &l, qf1, 1);
+       qa = l*(l-qd1)/(qd0*(qd0+qd1));
+       qb = (l+qd0)*(qd1-l)/(qd0*qd1);
+       qc = l*(l+qd0)/(qd1*(qd0+qd1));
+    }
+    else {
+       fx = qf1; qa = qb = 0.0; qc = 1.0;
+    }
+ #ifdef DEBUGPRAX
+   printf("after eventual min qd0=%14.8e qd1=%14.8e nl=%d\n",qd0, qd1,nl);
+ #endif
+    qd0 = qd1;
+    /* for (i=0; i<n; i++) { */
+    for (i=1; i<=n; i++) {
+        s = q0[i]; q0[i] = x[i];
+        x[i] = qa*s + qb*x[i] + qc*q1[i];
+    }
+ #ifdef DEBUGQUAD
+    vecprint ( " X after QUAD:" , x, n );
+ #endif
+ }
+ /* void minfit(int n, double eps, double tol, double ab[N][N], double q[]) */
+ void minfit(int n, double eps, double tol, double **ab, double q[])
+ /* int n; */
+ /* double eps, tol, ab[N][N], q[N]; */
+ {
+    int l, kt, l2, i, j, k;
+    double c, f, g, h, s, x, y, z;
+    /* double eps; */
+ /* #ifndef MSDOS */
+ /*    double e[N];              /\* plenty of stack on a vax *\/ */
+ /* #endif */
+    /* double *e; */
+    /* e=vector(0,n-1); /\* should be freed somewhere but gotos *\/ */
+    /* householder's reduction to bidiagonal form */
+    if(n==1){
+      /* q[1-1]=ab[1-1][1-1]; */
+      /* ab[1-1][1-1]=1.0; */
+      q[1]=ab[1][1];
+      ab[1][1]=1.0;
+      return; /* added from hardt */
+    }
+    /* eps=macheps; */ /* added */
+    x = g = 0.0;
+ #ifdef DEBUGPRAX
+    matprint (" HOUSE holder:", ab, n, n);
+ #endif
+    /* for (i=0; i<n; i++) {  /\* FOR I := 1 UNTIL N DO *\/ */
+    for (i=1; i<=n; i++) {  /* FOR I := 1 UNTIL N DO */
+      e[i] = g; s = 0.0; l = i+1;
+      /* for (j=i; j<n; j++)  /\* FOR J := I UNTIL N DO S := S*AB(J,I)**2; *\/ /\* not correct *\/ */
+      for (j=i; j<=n; j++)  /* FOR J := I UNTIL N DO S := S*AB(J,I)**2; */ /* not correct */
+        s += ab[j][i] * ab[j][i];
+ #ifdef DEBUGPRAXFIN
+      printf("i=%d s=%d %.7g tol=%.7g",i,s,tol);
+ #endif
+      if (s < tol) {
+        g = 0.0;
+      }
+      else {
+        /* f = ab[i][i]; */
+        f = ab[i][i];
+        if (f < 0.0)
+          g = sqrt(s);
+        else
+          g = -sqrt(s);
+        /* h = f*g - s; ab[i][i] = f - g; */
+        h = f*g - s; ab[i][i] = f - g;
+        /* for (j=l; j<n; j++) { */ /* FOR J := L UNTIL N DO */ /* wrong */
+        for (j=l; j<=n; j++) {
+          f = 0.0;
+          /* for (k=i; k<n; k++) /\* FOR K := I UNTIL N DO *\/ /\* wrong *\/ */
+          for (k=i; k<=n; k++) /* FOR K := I UNTIL N DO */
+            /* f += ab[k][i] * ab[k][j]; */
+            f += ab[k][i] * ab[k][j];
+          f /= h;
+          for (k=i; k<=n; k++) /* FOR K := I UNTIL N DO */
+            /* for (k=i; k<n; k++)/\* FOR K := I UNTIL N DO *\/ /\* wrong *\/ */
+            ab[k][j] += f * ab[k][i];
+          /* ab[k][j] += f * ab[k][i]; */
+ #ifdef DEBUGPRAX
+          printf("Holder J=%d F=%.7g",j,f);
+ #endif
+        }
+      } /* end s */
+      /* q[i] = g; s = 0.0; */
+      q[i] = g; s = 0.0;
+ #ifdef DEBUGPRAX
+      printf(" I Q=%d %.7g",i,q[i]);
+ #endif
+      /* if (i < n) */
+      /* if (i <= n)  /\* I is always lower or equal to n wasn't in golub reinsch*\/ */
+      /* for (j=l; j<n; j++) */
+      for (j=l; j<=n; j++)
+        s += ab[i][j] * ab[i][j];
+      /* s += ab[i][j] * ab[i][j]; */
+      if (s < tol) {
+        g = 0.0;
+      }
+      else {
+        if(i<n)
+          /* f = ab[i][i+1]; */ /* Brent golub overflow */
+          f = ab[i][i+1];
+        if (f < 0.0)
+          g = sqrt(s);
+        else
+          g = - sqrt(s);
+        h = f*g - s;
+        /* h = f*g - s; ab[i][i+1] = f - g; */ /* Overflow for i=n Error in Golub too but not Burkardt*/
+        /* for (j=l; j<n; j++) */
+        /*     e[j] = ab[i][j]/h; */
+        if(i<n){
+          ab[i][i+1] = f - g;
+          for (j=l; j<=n; j++)
+            e[j] = ab[i][j]/h;
+          /* for (j=l; j<n; j++) { */
+          for (j=l; j<=n; j++) {
+            s = 0.0;
+            /* for (k=l; k<n; k++) s += ab[j][k]*ab[i][k]; */
+            for (k=l; k<=n; k++) s += ab[j][k]*ab[i][k];
+            /* for (k=l; k<n; k++) ab[j][k] += s * e[k]; */
+            for (k=l; k<=n; k++) ab[j][k] += s * e[k];
+          } /* END J */
+        } /* END i <n */
+      } /* end s */
+        /* y = fabs(q[i]) + fabs(e[i]); */
+      y = fabs(q[i]) + fabs(e[i]);
+      if (y > x) x = y;
+ #ifdef DEBUGPRAX
+      printf(" I Y=%d %.7g",i,y);
+ #endif
+ #ifdef DEBUGPRAX
+      printf(" i=%d e(i) %.7g",i,e[i]);
+ #endif
+    } /* end i */
+    /*
+      Accumulation of right hand transformations */
+    /* for (i=n-1; i >= 0; i--) { */ /* FOR I := N STEP -1 UNTIL 1 DO */
+    /* We should avoid the overflow in Golub */
+    /* ab[n-1][n-1] = 1.0; */
+    /* g = e[n-1]; */
+    ab[n][n] = 1.0;
+    g = e[n];
+    l = n;
+    /* for (i=n; i >= 1; i--) { */
+    for (i=n-1; i >= 1; i--) { /* n-1 loops, different from brent and golub*/
+      if (g != 0.0) {
+        /* h = ab[i-1][i]*g; */
+        h = ab[i][i+1]*g;
+        for (j=l; j<=n; j++) ab[j][i] = ab[i][j] / h;
+        for (j=l; j<=n; j++) {
+          /* h = ab[i][i+1]*g; */
+          /* for (j=l; j<n; j++) ab[j][i] = ab[i][j] / h; */
+          /* for (j=l; j<n; j++) { */
+          s = 0.0;
+          /* for (k=l; k<n; k++) s += ab[i][k] * ab[k][j]; */
+          /* for (k=l; k<n; k++) ab[k][j] += s * ab[k][i]; */
+          for (k=l; k<=n; k++) s += ab[i][k] * ab[k][j];
+          for (k=l; k<=n; k++) ab[k][j] += s * ab[k][i];
+        }/* END J */
+      }/* END G */
+      /* for (j=l; j<n; j++) */
+      /*     ab[i][j] = ab[j][i] = 0.0; */
+      /* ab[i][i] = 1.0; g = e[i]; l = i; */
+      for (j=l; j<=n; j++)
+        ab[i][j] = ab[j][i] = 0.0;
+      ab[i][i] = 1.0; g = e[i]; l = i;
+    }/* END I */
+ #ifdef DEBUGPRAX
+    matprint (" HOUSE accumulation:",ab,n, n );
+ #endif
+    /* diagonalization to bidiagonal form */
+    eps *= x;
+    /* for (k=n-1; k>= 0; k--) { */
+    for (k=n; k>= 1; k--) {
+      kt = 0;
+ TestFsplitting:
+ #ifdef DEBUGPRAX
+      printf(" TestFsplitting: k=%d kt=%d\n",k,kt);
+      /* for(i=1;i<=n;i++)printf(" e(%d)=%.14f",i,e[i]);printf("\n"); */
+ #endif
+      kt = kt+1;
+ /* TestFsplitting: */
+      /* if (++kt > 30) { */
+      if (kt > 30) {
+        e[k] = 0.0;
+        fprintf(stderr, "\n+++ MINFIT - Fatal error\n");
+        fprintf ( stderr, "  The QR algorithm failed to converge.\n" );
+      }
+      /* for (l2=k; l2>=0; l2--) { */
+      for (l2=k; l2>=1; l2--) {
+        l = l2;
+ #ifdef DEBUGPRAX
+        printf(" l e(l)< eps %d %.7g %.7g ",l,e[l], eps);
+ #endif
+        /* if (fabs(e[l]) <= eps) */
+        if (fabs(e[l]) <= eps)
+          goto TestFconvergence;
+        /* if (fabs(q[l-1]) <= eps)*/ /* missing if ( 1 < l ){ *//* printf(" q(l-1)< eps %d %.7g %.7g ",l-1,q[l-2], eps); */
+        if (fabs(q[l-1]) <= eps)
+          break; /* goto Cancellation; */
+      }
+    Cancellation:
+ #ifdef DEBUGPRAX
+      printf(" Cancellation:\n");
+ #endif
+      c = 0.0; s = 1.0;
+      for (i=l; i<=k; i++) {
+        f = s * e[i]; e[i] *= c;
+        /* f = s * e[i]; e[i] *= c; */
+        if (fabs(f) <= eps)
+          goto TestFconvergence;
+        /* g = q[i]; */
+        g = q[i];
+        if (fabs(f) < fabs(g)) {
+          double fg = f/g;
+          h = fabs(g)*sqrt(1.0+fg*fg);
+        }
+        else {
+          double gf = g/f;
+          h = (f!=0.0 ? fabs(f)*sqrt(1.0+gf*gf) : 0.0);
+        }
+        /*    COMMENT: THE ABOVE REPLACES Q(I):=H:=LONGSQRT(G*G+F*F) */
+        /* WHICH MAY GIVE INCORRECT RESULTS IF THE */
+        /* SQUARES UNDERFLOW OR IF F = G = 0; */
+        /* q[i] = h; */
+        q[i] = h;
+        if (h == 0.0) { h = 1.0; g = 1.0; }
+        c = g/h; s = -f/h;
+      }
+ TestFconvergence:
+  #ifdef DEBUGPRAX
+      printf(" TestFconvergence: l=%d k=%d\n",l,k);
+ #endif
+      /* z = q[k]; */
+      z = q[k];
+      if (l == k)
+        goto Convergence;
+      /* shift from bottom 2x2 minor */
+      /* x = q[l]; y = q[k-l]; g = e[k-1]; h = e[k]; */ /* Error */
+      x = q[l]; y = q[k-1]; g = e[k-1]; h = e[k];
+      f = ((y-z)*(y+z) + (g-h)*(g+h)) / (2.0*h*y);
+      g = sqrt(f*f+1.0);
+      if (f <= 0.0)
+        f = ((x-z)*(x+z) + h*(y/(f-g)-h))/x;
+      else
+        f = ((x-z)*(x+z) + h*(y/(f+g)-h))/x;
+      /* next qr transformation */
+      s = c = 1.0;
+      for (i=l+1; i<=k; i++) {
+ #ifdef DEBUGPRAXQR
+        printf(" Before Mid TestFconvergence: l+1=%d i=%d k=%d h=%.6e e(i)=%14.8f e(i-1)=%14.8f\n",l+1,i,k, h, e[i],e[i-1]);
+ #endif
+        /* g = e[i]; y = q[i]; h = s*g; g *= c; */
+        g = e[i]; y = q[i]; h = s*g; g *= c;
+        if (fabs(f) < fabs(h)) {
+          double fh = f/h;
+          z = fabs(h) * sqrt(1.0 + fh*fh);
+        }
+        else {
+          double hf = h/f;
+          z = (f!=0.0 ? fabs(f)*sqrt(1.0+hf*hf) : 0.0);
+        }
+        /* e[i-1] = z; */
+        e[i-1] = z;
+ #ifdef DEBUGPRAXQR
+        printf(" Mid TestFconvergence: l+1=%d i=%d k=%d h=%.6e e(i)=%14.8f e(i-1)=%14.8f\n",l+1,i,k, h, e[i],e[i-1]);
+ #endif
+        if (z == 0.0)
+          f = z = 1.0;
+        c = f/z; s = h/z;
+        f = x*c + g*s; g = - x*s + g*c; h = y*s;
+        y *= c;
+        /* for (j=0; j<n; j++) { */
+        /*     x = ab[j][i-1]; z = ab[j][i]; */
+        /*     ab[j][i-1] = x*c + z*s; */
+        /*     ab[j][i] = - x*s + z*c; */
+        /* } */
+        for (j=1; j<=n; j++) {
+          x = ab[j][i-1]; z = ab[j][i];
+          ab[j][i-1] = x*c + z*s;
+          ab[j][i] = - x*s + z*c;
+        }
+        if (fabs(f) < fabs(h)) {
+          double fh = f/h;
+          z = fabs(h) * sqrt(1.0 + fh*fh);
+        }
+        else {
+          double hf = h/f;
+          z = (f!=0.0 ? fabs(f)*sqrt(1.0+hf*hf) : 0.0);
+        }
+ #ifdef DEBUGPRAXQR
+        printf(" qr transformation z f h=%.7g %.7g %.7g i=%d k=%d\n",z,f,h, i, k);
+ #endif
+        q[i-1] = z;
+        if (z == 0.0)
+          z = f = 1.0;
+        c = f/z; s = h/z;
+        f = c*g + s*y;  /* f can be very small */
+        x = - s*g + c*y;
+      }
+      /* e[l] = 0.0; e[k] = f; q[k] = x; */
+      e[l] = 0.0; e[k] = f; q[k] = x;
+ #ifdef DEBUGPRAXQR
+      printf(" aftermid loop l=%d k=%d e(l)=%7g e(k)=%.7g q(k)=%.7g x=%.7g\n",l,k,e[l],e[k],q[k],x);
+ #endif
+      goto TestFsplitting;
+    Convergence:
+ #ifdef DEBUGPRAX
+      printf(" Convergence:\n");
+ #endif
+      if (z < 0.0) {
+        /* q[k] = - z; */
+        /* for (j=0; j<n; j++) ab[j][k] = - ab[j][k]; */
+        q[k] = - z;
+        for (j=1; j<=n; j++) ab[j][k] = - ab[j][k];
+      }/* END Z */
+    }/* END K */
+ } /* END MINFIT */
+ double praxis(double tol, double macheps, double h0, int _n, int _prin, double *_x, double (*_fun)(double *_x))
+ /* double praxis(double tol, double macheps, double h0, int _n, int _prin, double *_x, double (*_fun)(double *_x, int _n)) */
+ /* double praxis(double (*_fun)(), double _x[], int _n) */
+ /* double (*_fun)(); */
+ /* double _x[N]; */
+ /* double (*_fun)(); */
+ /* double _x[N]; */
+ {
+    /* init global extern variables and parameters */
+    /* double *d, *y, *z, */
+    /*   *q0, *q1, **v; */
+    /* double *tflin; /\* used in flin: return (*fun)(tflin, n); *\/ */
+    /* double *e; /\* used in minfit, don't konw how to free memory and thus made global *\/ */
+   int seed; /* added */
+   int biter=0;
+   double r;
+   double randbrent( int (*));
+   double s, sf;
+    h = h0; /* step; */
+    t = tol;
+    scbd = 1.0;
+    illc = 0;
+    ktm = 1;
+    macheps = DBL_EPSILON;
+    /* prin=4; */
+ #ifdef DEBUGPRAX
+    printf("Praxis macheps=%14g h=%14g step=%14g tol=%14g\n",macheps,h, h0,tol);
+ #endif
+    n = _n;
+    x = _x;
+    prin = _prin;
+    fun = _fun;
+    d=vector(1, n);
+    y=vector(1, n);
+    z=vector(1, n);
+    q0=vector(1, n);
+    q1=vector(1, n);
+    e=vector(1, n);
+    tflin=vector(1, n);
+    v=matrix(1, n, 1, n);
+    for(i=1;i<=n;i++){d[i]=y[i]=z[i]=q0[0]=e[i]=tflin[i]=0.;}
+    small_windows = (macheps) * (macheps); vsmall = small_windows*small_windows;
+    large = 1.0/small_windows; vlarge = 1.0/vsmall;
+    m2 = sqrt(macheps); m4 = sqrt(m2);
+    seed = 123456789; /* added */
+    ldfac = (illc ? 0.1 : 0.01);
+    for(i=1;i<=n;i++) z[i]=0.; /* Was missing in Gegenfurtner as well as Brent's algol or fortran  */
+    nl = kt = 0; nf = 1;
+ #ifdef NR_SHIFT
+    fx = (*fun)((x-1), n);
+ #else
+    fx = (*fun)(x);
+ #endif
+    qf1 = fx;
+    t2 = small_windows + fabs(t); t = t2; dmin = small_windows;
+ #ifdef DEBUGPRAX
+    printf("praxis2 macheps=%14g h=%14g step=%14g small=%14g t=%14g\n",macheps,h, h0,small_windows, t);
+ #endif
+    if (h < 100.0*t) h = 100.0*t;
+ #ifdef DEBUGPRAX
+    printf("praxis3 macheps=%14g h=%14g step=%14g small=%14g t=%14g\n",macheps,h, h0,small_windows, t);
+ #endif
+    ldt = h;
+    /* for (i=0; i<n; i++) for (j=0; j<n; j++) */
+    for (i=1; i<=n; i++) for (j=1; j<=n; j++)
+        v[i][j] = (i == j ? 1.0 : 0.0);
+    d[1] = 0.0; qd0 = 0.0;
+    /* for (i=0; i<n; i++) q1[i] = x[i]; */
+    for (i=1; i<=n; i++) q1[i] = x[i];
+    if (prin > 1) {
+       printf("\n------------- enter function praxis -----------\n");
+       printf("... current parameter settings ...\n");
+       printf("... scaling ... %20.10e\n", scbd);
+       printf("...   tol   ... %20.10e\n", t);
+       printf("... maxstep ... %20.10e\n", h);
+       printf("...   illc  ... %20u\n", illc);
+       printf("...   ktm   ... %20u\n", ktm);
+       printf("... maxfun  ... %20u\n", maxfun);
+    }
+    if (prin) print2();
+ mloop:
+     biter++;  /* Added to count the loops */
+    /* sf = d[0]; */
+    /* s = d[0] = 0.0; */
+     printf("\n Big iteration %d \n",biter);
+     fprintf(ficlog,"\n Big iteration %d \n",biter);
+     sf = d[1];
+    s = d[1] = 0.0;
+    /* minimize along first direction V(*,1) */
+ #ifdef DEBUGPRAX
+    printf("  Minimize along the first direction V(*,1). illc=%d\n",illc);
+    /* fprintf(ficlog,"  Minimize along the first direction V(*,1).\n"); */
+ #endif
+ #ifdef DEBUGPRAX2
+    printf("praxis4 macheps=%14g h=%14g step=%14g small=%14g t=%14g\n",macheps,h, h0,small_windows, t);
+ #endif
+    /* min(0, 2, &d[0], &s, fx, 0); /\* mac heps not global *\/ */
+    minny(1, 2, &d[1], &s, fx, 0); /* mac heps not global it seems that fx doesn't correspond to f(s=*x1) */
+ #ifdef DEBUGPRAX
+    printf("praxis5 macheps=%14g h=%14g looks at sign of s=%14g fx=%14g\n",macheps,h, s,fx);
+ #endif
+    if (s <= 0.0)
+       /* for (i=0; i < n; i++) */
+       for (i=1; i <= n; i++)
+           v[i][1] = -v[i][1];
+    /* if ((sf <= (0.9 * d[0])) || ((0.9 * sf) >= d[0])) */
+    if ((sf <= (0.9 * d[1])) || ((0.9 * sf) >= d[1]))
+       /* for (i=1; i<n; i++) */
+       for (i=2; i<=n; i++)
+           d[i] = 0.0;
+    /* for (k=1; k<n; k++) { */
+    for (k=2; k<=n; k++) {
+     /*
+       The inner loop starts here.
+     */
+ #ifdef DEBUGPRAX
+       printf("      The inner loop  here from k=%d to n=%d.\n",k,n);
+       /* fprintf(ficlog,"      The inner loop  here from k=%d to n=%d.\n",k,n); */
+ #endif
+        /* for (i=0; i<n; i++) */
+        for (i=1; i<=n; i++)
+            y[i] = x[i];
+        sf = fx;
+ #ifdef DEBUGPRAX
+        printf(" illc=%d and kt=%d and ktm=%d\n", illc, kt, ktm);
+ #endif
+        illc = illc || (kt > 0);
+ next:
+        kl = k;
+        df = 0.0;
+        if (illc) {        /* random step to get off resolution valley */
+ #ifdef DEBUGPRAX
+           printf("  A random step follows, to avoid resolution valleys.\n");
+           matprint("  before rand, vectors:",v,n,n);
+ #endif
+           for (i=1; i<=n; i++) {
+ #ifdef NOBRENTRAND
+             r = drandom();
+ #else
+             seed=i;
+             /* seed=i+1; */
+ #ifdef DEBUGRAND
+             printf(" Random seed=%d, brent i=%d",seed,i); /* YYYY i=5 j=1 vji= -0.0001170073 */
+ #endif
+             r = randbrent ( &seed );
+ #endif
+ #ifdef DEBUGRAND
+             printf(" Random r=%.7g \n",r);
+ #endif
+             z[i] = (0.1 * ldt + t2 * pow(10.0,(double)kt)) * (r - 0.5);
+             /* z[i] = (0.1 * ldt + t2 * pow(10.0,(double)kt)) * (drandom() - 0.5); */
+             s = z[i];
+               for (j=1; j <= n; j++)
+                   x[j] += s * v[j][i];
+           }
+ #ifdef DEBUGRAND
+           matprint("  after rand, vectors:",v,n,n);
+ #endif
+ #ifdef NR_SHIFT
+           fx = (*fun)((x-1), n);
+ #else
+           fx = (*fun)(x, n);
+ #endif
+           /* fx = (*func) ( (x-1) ); *//* This for func which is computed from x[1] and not from x[0] xm1=(x-1)*/
+           nf++;
+        }
+        /* minimize along non-conjugate directions */
+ #ifdef DEBUGPRAX
+         printf(" Minimize along the 'non-conjugate' directions (dots printed) V(*,%d),...,V(*,%d).\n",k,n);
+         /* fprintf(ficlog," Minimize along the 'non-conjugate' directions  (dots printed) V(*,%d),...,V(*,%d).\n",k,n); */
+ #endif
+         /* for (k2=k; k2<n; k2++) {  /\* Be careful here k2 <=n ? *\/ */
+         for (k2=k; k2<=n; k2++) {  /* Be careful here k2 <=n ? */
+            sl = fx;
+            s = 0.0;
+ #ifdef DEBUGPRAX
+            printf(" Minimize along the 'NON-CONJUGATE' true direction k2=%14d fx=%14.7f\n",k2, fx);
+    matprint("  before min vectors:",v,n,n);
+ #endif
+            /* min(k2, 2, &d[k2], &s, fx, 0); */
+    /*     jsearch=k2-1; */
+    /* min(jsearch, 2, &d[jsearch], &s, fx, 0); */
+    minny(k2, 2, &d[k2], &s, fx, 0);
+ #ifdef DEBUGPRAX
+            printf(" . D(%d)=%14.7f d[k2]=%14.7f z[k2]=%14.7f illc=%14d fx=%14.7f\n",k2,d[k2],d[k2],z[k2],illc,fx);
+ #endif
+           if (illc) {
+               /* double szk = s + z[k2]; */
+               /* s = d[k2] * szk*szk; */
+               double szk = s + z[k2];
+               s = d[k2] * szk*szk;
+            }
+            else
+               s = sl - fx;
+            /* if (df < s) { */
+            if (df <= s) {
+               df = s;
+               kl = k2;
+ #ifdef DEBUGPRAX
+             printf(" df=%.7g and choose kl=%d \n",df,kl); /* UUUU */
+ #endif
+            }
+         } /* end loop k2 */
+         /*
+           If there was not much improvement on the first try, set
+           ILLC = true and start the inner loop again.
+         */
+ #ifdef DEBUGPRAX
+         printf(" If there was not much improvement on the first try, set ILLC = true and start the inner loop again. illc=%d\n",illc);
+         /* fprintf(ficlog,"  If there was not much improvement on the first try, set ILLC = true and start the inner loop again.\n"); */
+ #endif
+         if (!illc && (df < fabs(100.0 * (macheps) * fx))) {
+ #ifdef DEBUGPRAX
+           printf("\n NO SUCCESS because DF is small, starts inner loop with same K(=%d), fabs(  100.0 * machep(=%.10e) * fx(=%.9e) )=%.9e > df(=%.9e) break illc=%d\n", k, macheps, fx, fabs ( 100.0 * macheps * fx ), df, illc);
+ #endif
+           illc = 1;
+           goto next;
+         }
+ #ifdef DEBUGPRAX
+         printf("\n SUCCESS, BREAKS inner loop K(=%d) because DF is big, fabs(  100.0 * machep(=%.10e) * fx(=%.9e) )=%.9e <= df(=%.9e) break illc=%d\n", k, macheps, fx, fabs ( 100.0 * macheps * fx ), df, illc);
+ #endif
+        /* if ((k == 1) && (prin > 1)){ /\* be careful k=2 *\/ */
+        if ((k == 2) && (prin > 1)){ /* be careful k=2 */
+ #ifdef DEBUGPRAX
+         printf("  NEW D The second difference array d:\n" );
+         /* fprintf(ficlog, " NEW D The second difference array d:\n" ); */
+ #endif
+          vecprint(" NEW D The second difference array d:",d,n);
+        }
+        /* minimize along conjugate directions */
+        /*
+          Minimize along the "conjugate" directions V(*,1),...,V(*,K-1).
+        */
+ #ifdef DEBUGPRAX
+       printf("Minimize along the 'conjugate' directions V(*,1),...,V(*,K-1=%d).\n",k-1);
+       /* fprintf(ficlog,"Minimize along the 'conjugate' directions V(*,1),...,V(*,K-1=%d).\n",k-1); */
+ #endif
+       /* for (k2=0; k2<=k-1; k2++) { */
+       for (k2=1; k2<=k-1; k2++) {
+            s = 0.0;
+            /* min(k2-1, 2, &d[k2-1], &s, fx, 0); */
+            minny(k2, 2, &d[k2], &s, fx, 0);
+        }
+        f1 = fx;
+        fx = sf;
+        lds = 0.0;
+        /* for (i=0; i<n; i++) { */
+        for (i=1; i<=n; i++) {
+            sl = x[i];
+            x[i] = y[i];
+            y[i] = sl - y[i];
+            sl = y[i];
+            lds = lds + sl*sl;
+        }
+        lds = sqrt(lds);
+ #ifdef DEBUGPRAX
+        printf("Minimization done 'conjugate', shifted all points, computed lds=%.8f\n",lds);
+ #endif
+       /*
+         Discard direction V(*,kl).
+         If no random step was taken, V(*,KL) is the "non-conjugate"
+         direction along which the greatest improvement was made.
+       */
+        if (lds > small_windows) {
+ #ifdef DEBUGPRAX
+        printf("lds big enough to throw direction  V(*,kl=%d). If no random step was taken, V(*,KL) is the 'non-conjugate' direction along which the greatest improvement was made.\n",kl);
+          matprint("  before shift new conjugate vectors:",v,n,n);
+ #endif
+          for (i=kl-1; i>=k; i--) {
+            /* for (j=0; j < n; j++) */
+            for (j=1; j <= n; j++)
+              /* v[j][i+1] = v[j][i]; */ /* This is v[j][i+1]=v[j][i] i=kl-1 to k */
+              v[j][i+1] = v[j][i]; /* This is v[j][i+1]=v[j][i] i=kl-1 to k */
+            /* v[j][i+1] = v[j][i]; */
+            /* d[i+1] = d[i];*/  /* last  is d[k+1]= d[k] */
+            d[i+1] = d[i];  /* last  is d[k]= d[k-1] */
+          }
+ #ifdef DEBUGPRAX
+          matprint("  after shift new conjugate vectors:",v,n,n);
+ #endif   /* d[k] = 0.0; */
+          d[k] = 0.0;
+          for (i=1; i <= n; i++)
+            v[i][k] = y[i] / lds;
+          /* v[i][k] = y[i] / lds; */
+ #ifdef DEBUGPRAX
+          printf("Minimize along the new 'conjugate' direction V(*,k=%d), which is the normalized vector:  (new x) - (old x). d2=%14.7g lds=%.10f\n",k,d[k],lds);
+          /* fprintf(ficlog,"Minimize along the new 'conjugate' direction V(*,k=%d), which is the normalized vector:  (new x) - (old x).\n",k); */
+     matprint("  before min new conjugate vectors:",v,n,n);
+ #endif
+          /* min(k-1, 4, &d[k-1], &lds, f1, 1); */
+          minny(k, 4, &d[k], &lds, f1, 1);
+ #ifdef DEBUGPRAX
+          printf(" after min d(k)=%d %.7g lds=%14f\n",k,d[k],lds);
+    matprint("  after min vectors:",v,n,n);
+ #endif
+          if (lds <= 0.0) {
+            lds = -lds;
+ #ifdef DEBUGPRAX
+           printf(" lds changed sign lds=%.14f k=%d\n",lds,k);
+ #endif
+            /* for (i=0; i<n; i++) */
+            /*   v[i][k] = -v[i][k]; */
+            for (i=1; i<=n; i++)
+              v[i][k] = -v[i][k];
+          }
+        }
+        ldt = ldfac * ldt;
+        if (ldt < lds)
+           ldt = lds;
+        if (prin > 0){
+ #ifdef DEBUGPRAX
+         printf(" k=%d",k);
+         /* fprintf(ficlog," k=%d",k); */
+ #endif
+         print2();/* n, x, prin, fx, nf, nl ); */
+        }
+        t2 = 0.0;
+        /* for (i=0; i<n; i++) */
+        for (i=1; i<=n; i++)
+            t2 += x[i]*x[i];
+        t2 = m2 * sqrt(t2) + t;
+        /*
+         See whether the length of the step taken since starting the
+         inner loop exceeds half the tolerance.
+       */
+ #ifdef DEBUGPRAX
+        printf("See if step length exceeds half the tolerance.\n"); /* ZZZZZ */
+       /* fprintf(ficlog,"See if step length exceeds half the tolerance.\n"); */
+ #endif
+        if (ldt > (0.5 * t2))
+           kt = 0;
+        else
+           kt++;
+ #ifdef DEBUGPRAX
+        printf("if kt=%d >? ktm=%d gotoL2 loop\n",kt,ktm);
+ #endif
+        if (kt > ktm){
+          if ( 0 < prin ){
+            /* printf("\nr8vec_print\n X:\n"); */
+            /* fprintf(ficlog,"\nr8vec_print\n X:\n"); */
+            vecprint ("END  X:", x, n );
+          }
+            goto fret;
+        }
+ #ifdef DEBUGPRAX
+    matprint("  end of L2 loop vectors:",v,n,n);
+ #endif
+    }
+    /* printf("The inner loop ends here.\n"); */
+    /* fprintf(ficlog,"The inner loop ends here.\n"); */
+    /*
+      The inner loop ends here.
+      Try quadratic extrapolation in case we are in a curved valley.
+    */
+ #ifdef DEBUGPRAX
+    printf("Try QUAD ratic extrapolation in case we are in a curved valley.\n");
+ #endif
+    /*  try quadratic extrapolation in case    */
+    /*  we are stuck in a curved valley        */
+    quad();
+    dn = 0.0;
+    /* for (i=0; i<n; i++) { */
+    for (i=1; i<=n; i++) {
+        d[i] = 1.0 / sqrt(d[i]);
+        if (dn < d[i])
+           dn = d[i];
+    }
+    if (prin > 2)
+      matprint("  NEW DIRECTIONS vectors:",v,n,n);
+    /* for (j=0; j<n; j++) { */
+    for (j=1; j<=n; j++) {
+        s = d[j] / dn;
+        /* for (i=0; i < n; i++) */
+        for (i=1; i <= n; i++)
+            v[i][j] *= s;
+    }
+    if (scbd > 1.0) {       /* scale axis to reduce condition number */
+ #ifdef DEBUGPRAX
+      printf("Scale the axes to try to reduce the condition number.\n");
+ #endif
+      /* fprintf(ficlog,"Scale the axes to try to reduce the condition number.\n"); */
+       s = vlarge;
+       /* for (i=0; i<n; i++) { */
+       for (i=1; i<=n; i++) {
+           sl = 0.0;
+           /* for (j=0; j < n; j++) */
+           for (j=1; j <= n; j++)
+               sl += v[i][j]*v[i][j];
+           z[i] = sqrt(sl);
+           if (z[i] < m4)
+              z[i] = m4;
+           if (s > z[i])
+              s = z[i];
+       }
+       /* for (i=0; i<n; i++) { */
+       for (i=1; i<=n; i++) {
+           sl = s / z[i];
+           z[i] = 1.0 / sl;
+           if (z[i] > scbd) {
+              sl = 1.0 / scbd;
+              z[i] = scbd;
+           }
+       }
+    }
+    for (i=1; i<=n; i++)
+        /* for (j=0; j<=i-1; j++) { */
+        /* for (j=1; j<=i; j++) { */
+        for (j=1; j<=i-1; j++) {
+            s = v[i][j];
+            v[i][j] = v[j][i];
+            v[j][i] = s;
+        }
+ #ifdef DEBUGPRAX
+     printf(" Calculate a new set of orthogonal directions before repeating  the main loop.\n  Transpose V for MINFIT:...\n");
+ #endif
+       /*
+       MINFIT finds the singular value decomposition of V.
+       This gives the principal values and principal directions of the
+       approximating quadratic form without squaring the condition number.
+     */
+  #ifdef DEBUGPRAX
+     printf(" MINFIT finds the singular value decomposition of V. \n This gives the principal values and principal directions of the\n  approximating quadratic form without squaring the condition number...\n");
+ #endif
+    minfit(n, macheps, vsmall, v, d);
+     /* for(i=0; i<n;i++)printf(" %14.7g",d[i]); */
+     /* v is overwritten with R. */
+     /*
+       Unscale the axes.
+     */
+    if (scbd > 1.0) {
+ #ifdef DEBUGPRAX
+       printf(" Unscale the axes.\n");
+ #endif
+       /* for (i=0; i<n; i++) { */
+       for (i=1; i<=n; i++) {
+           s = z[i];
+           /* for (j=0; j<n; j++) */
+           for (j=1; j<=n; j++)
+               v[i][j] *= s;
+       }
+       /* for (i=0; i<n; i++) { */
+       for (i=1; i<=n; i++) {
+           s = 0.0;
+           /* for (j=0; j<n; j++) */
+           for (j=1; j<=n; j++)
+               s += v[j][i]*v[j][i];
+           s = sqrt(s);
+           d[i] *= s;
+           s = 1.0 / s;
+           /* for (j=0; j<n; j++) */
+           for (j=1; j<=n; j++)
+               v[j][i] *= s;
+       }
+    }
+    /* for (i=0; i<n; i++) { */
+    double dni; /* added for compatibility with buckhardt but not brent */
+    for (i=1; i<=n; i++) {
+      dni=dn*d[i]; /* added for compatibility with buckhardt but not brent */
+        if ((dn * d[i]) > large)
+           d[i] = vsmall;
+        else if ((dn * d[i]) < small_windows)
+           d[i] = vlarge;
+        else
+         d[i] = 1.0 / dni / dni; /* added for compatibility with buckhardt but not brent */
+           /* d[i] = pow(dn * d[i],-2.0); */
+    }
+ #ifdef DEBUGPRAX
+    vecprint ("\n Before sort Eigenvalues of a:",d,n );
+ #endif
+    sort();               /* the new eigenvalues and eigenvectors */
+ #ifdef DEBUGPRAX
+    vecprint( " After sort the eigenvalues ....\n", d, n);
+    matprint( " After sort the eigenvectors....\n", v, n,n);
+ #endif
+ #ifdef DEBUGPRAX
+     printf("  Determine the smallest eigenvalue.\n");
+ #endif
+    /* dmin = d[n-1]; */
+    dmin = d[n];
+    if (dmin < small_windows)
+       dmin = small_windows;
+     /*
+      The ratio of the smallest to largest eigenvalue determines whether
+      the system is ill conditioned.
+    */
+    /* illc = (m2 * d[0]) > dmin; */
+    illc = (m2 * d[1]) > dmin;
+ #ifdef DEBUGPRAX
+     printf("  The ratio of the smallest to largest eigenvalue determines whether\n  the system is ill conditioned=%d . dmin=%.10lf < m2=%.10lf * d[1]=%.10lf \n",illc, dmin,m2, d[1]);
+ #endif
+    if ((prin > 2) && (scbd > 1.0))
+       vecprint("\n The scale factors:",z,n);
+    if (prin > 2)
+       vecprint("  Principal values (EIGEN VALUES OF A) of the quadratic form:",d,n);
+    if (prin > 2)
+      matprint("  The principal axes (EIGEN VECTORS OF A:",v,n, n);
+    if ((maxfun > 0) && (nf > maxfun)) {
+       if (prin)
+          printf("\n... maximum number of function calls reached ...\n");
+       goto fret;
+    }
+ #ifdef DEBUGPRAX
+    printf("Goto main loop\n");
+ #endif
+    goto mloop;   /* back to main loop */
+ fret:
+    if (prin > 0) {
+          vecprint("\n  X:", x, n);
+          /* printf("\n... ChiSq reduced to %20.10e ...\n", fx); */
+          /* printf("... after %20u function calls.\n", nf); */
+    }
+    free_vector(d, 1, n);
+    free_vector(y, 1, n);
+    free_vector(z, 1, n);
+    free_vector(q0, 1, n);
+    free_vector(q1, 1, n);
+    free_matrix(v, 1, n, 1, n);
+    /*   double *d, *y, *z, */
+    /* *q0, *q1, **v; */
+    free_vector(tflin, 1, n);
+    /* double *tflin; /\* used in flin: return (*fun)(tflin, n); *\/ */
+    free_vector(e, 1, n);
+    /* double *e; /\* used in minfit, don't konw how to free memory and thus made global *\/ */
+    return(fx);
+ }
+ /* end praxis gegen */
  /*************** powell ************************/
  /*
- Line 2660  void powell(double p[], double **xi, int
+ Line 4229  void powell(double p[], double **xi, int
      curr_time = *localtime(&rcurr_time);
      /* printf("\nPowell iter=%d -2*LL=%.12f gain=%.12f=%.3g %ld sec. %ld sec.",*iter,*fret, fp-*fret,fp-*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout); */
      /* fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f gain=%.12f=%.3g %ld sec. %ld sec.",*iter,*fret, fp-*fret,fp-*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog); */
-     Bigter=(*iter - *iter % ncovmodel)/ncovmodel +1; /* Big iteration, i.e on ncovmodel cycle */
+     /* Bigter=(*iter - *iter % ncovmodel)/ncovmodel +1; /\* Big iteration, i.e on ncovmodel cycle *\/ */
+     Bigter=(*iter - (*iter-1) % n)/n +1; /* Big iteration, i.e on ncovmodel cycle */
      printf("\nPowell iter=%d Big Iter=%d -2*LL=%.12f gain=%.3lg %ld sec. %ld sec.",*iter,Bigter,*fret,fp-*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout);
      fprintf(ficlog,"\nPowell iter=%d Big Iter=%d -2*LL=%.12f gain=%.3lg %ld sec. %ld sec.",*iter,Bigter,*fret,fp-*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog);
      fprintf(ficrespow,"%d %d %.12f %d",*iter,Bigter, *fret,curr_time.tm_sec-start_time.tm_sec);
- Line 2675  void powell(double p[], double **xi, int
+ Line 4245  void powell(double p[], double **xi, int
          printf("  + age*age  ");
          fprintf(ficlog,"  + age*age  ");
      }
-     for(j=1;j <=ncovmodel-2;j++){
+     for(j=1;j <=ncovmodel-2-nagesqr;j++){
        if(Typevar[j]==0) {
          printf("  +      V%d  ",Tvar[j]);
          fprintf(ficlog,"  +      V%d  ",Tvar[j]);
- Line 2731  void powell(double p[], double **xi, int
+ Line 4301  void powell(double p[], double **xi, int
          fprintf(ficlog,"   - if your program needs %d BIG iterations  (%d iterations) to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",nBigterf, niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
        }
      }
-     for (i=1;i<=n;i++) { /* For each direction i */
+     for (i=1;i<=n;i++) { /* For each direction i, maximisation after loading directions */
-       for (j=1;j<=n;j++) xit[j]=xi[j][i]; /* Directions stored from previous iteration with previous scales */
+       for (j=1;j<=n;j++) xit[j]=xi[j][i]; /* Directions stored from previous iteration with previous scales. xi is not changed but one dim xit  */
-       fptt=(*fret);
+       fptt=(*fret); /* Computes likelihood for parameters xit */
  #ifdef DEBUG
        printf("fret=%lf, %lf, %lf \n", *fret, *fret, *fret);
        fprintf(ficlog, "fret=%lf, %lf, %lf \n", *fret, *fret, *fret);
- Line 2741  void powell(double p[], double **xi, int
+ Line 4312  void powell(double p[], double **xi, int
        printf("%d",i);fflush(stdout); /* print direction (parameter) i */
        fprintf(ficlog,"%d",i);fflush(ficlog);
  #ifdef LINMINORIGINAL
-       linmin(p,xit,n,fret,func); /* Point p[n]. xit[n] has been loaded for direction i as input.*/
+       linmin(p,xit,n,fret,func); /* New point i minimizing in direction xit, i has coordinates p[j].*/
+       /* xit[j] gives the n coordinates of direction i as input.*/
+       /* *fret gives the maximum value on direction xit */
  #else
        linmin(p,xit,n,fret,func,&flat); /* Point p[n]. xit[n] has been loaded for direction i as input.*/
-                         flatdir[i]=flat; /* Function is vanishing in that direction i */
+       flatdir[i]=flat; /* Function is vanishing in that direction i */
  #endif
-                         /* Outputs are fret(new point p) p is updated and xit rescaled */
+       /* Outputs are fret(new point p) p is updated and xit rescaled */
        if (fabs(fptt-(*fret)) > del) { /* We are keeping the max gain on each of the n directions */
-                                 /* because that direction will be replaced unless the gain del is small */
+         /* because that direction will be replaced unless the gain del is small */
-                                 /* in comparison with the 'probable' gain, mu^2, with the last average direction. */
+         /* in comparison with the 'probable' gain, mu^2, with the last average direction. */
-                                 /* Unless the n directions are conjugate some gain in the determinant may be obtained */
+         /* Unless the n directions are conjugate some gain in the determinant may be obtained */
-                                 /* with the new direction. */
+         /* with the new direction. */
-                                 del=fabs(fptt-(*fret));
+         del=fabs(fptt-(*fret));
-                                 ibig=i;
+         ibig=i;
        }
  #ifdef DEBUG
        printf("%d %.12e",i,(*fret));
        fprintf(ficlog,"%d %.12e",i,(*fret));
        for (j=1;j<=n;j++) {
-                                 xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5);
+         xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5);
-                                 printf(" x(%d)=%.12e",j,xit[j]);
+         printf(" x(%d)=%.12e",j,xit[j]);
-                                 fprintf(ficlog," x(%d)=%.12e",j,xit[j]);
+         fprintf(ficlog," x(%d)=%.12e",j,xit[j]);
        }
        for(j=1;j<=n;j++) {
-                                 printf(" p(%d)=%.12e",j,p[j]);
+         printf(" p(%d)=%.12e",j,p[j]);
-                                 fprintf(ficlog," p(%d)=%.12e",j,p[j]);
+         fprintf(ficlog," p(%d)=%.12e",j,p[j]);
        }
        printf("\n");
        fprintf(ficlog,"\n");
  #endif
      } /* end loop on each direction i */
-     /* Convergence test will use last linmin estimation (fret) and compare former iteration (fp) */
+     /* Convergence test will use last linmin estimation (fret) and compare to former iteration (fp) */
      /* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit  */
      /* New value of last point Pn is not computed, P(n-1) */
      for(j=1;j<=n;j++) {
- Line 2826  void powell(double p[], double **xi, int
+ Line 4399  void powell(double p[], double **xi, int
        return;
      } /* enough precision */
      if (*iter == ITMAX*n) nrerror("powell exceeding maximum iterations.");
-     for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */
+     for (j=1;j<=n;j++) { /* Computes the extrapolated point and value f3, P_0 + 2 (P_n-P_0)=2Pn-P0 and xit is direction Pn-P0 */
        ptt[j]=2.0*p[j]-pt[j];
-       xit[j]=p[j]-pt[j];
+       xit[j]=p[j]-pt[j]; /* Coordinate j of last direction xi_n=P_n-P_0 */
-       pt[j]=p[j];
+ #ifdef DEBUG
-     }
+       printf("\n %d xit=%12.7g p=%12.7g pt=%12.7g ",j,xit[j],p[j],pt[j]);
+ #endif
+       pt[j]=p[j]; /* New P0 is Pn */
+     }
+ #ifdef DEBUG
+     printf("\n");
+ #endif
      fptt=(*func)(ptt); /* f_3 */
- #ifdef NODIRECTIONCHANGEDUNTILNITER  /* No change in drections until some iterations are done */
+ #ifdef NODIRECTIONCHANGEDUNTILNITER  /* No change in directions until some iterations are done */
                  if (*iter <=4) {
  #else
  #endif
- Line 2851  void powell(double p[], double **xi, int
+ Line 4430  void powell(double p[], double **xi, int
        /* t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); */
        /*  Even if f3 <f1, directest can be negative and t >0 */
        /* mu² and del² are equal when f3=f1 */
-                         /* f3 < f1 : mu² < del <= lambda^2 both test are equivalent */
+       /* f3 < f1 : mu² < del <= lambda^2 both test are equivalent */
-                         /* f3 < f1 : mu² < lambda^2 < del then directtest is negative and powell t is positive */
+       /* f3 < f1 : mu² < lambda^2 < del then directtest is negative and powell t is positive */
-                         /* f3 > f1 : lambda² < mu^2 < del then t is negative and directest >0  */
+       /* f3 > f1 : lambda² < mu^2 < del then t is negative and directest >0  */
-                         /* f3 > f1 : lambda² < del < mu^2 then t is positive and directest >0  */
+       /* f3 > f1 : lambda² < del < mu^2 then t is positive and directest >0  */
  #ifdef NRCORIGINAL
        t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)- del*SQR(fp-fptt); /* Original Numerical Recipes in C*/
  #else
- Line 2874  void powell(double p[], double **xi, int
+ Line 4453  void powell(double p[], double **xi, int
  #endif
  #ifdef POWELLORIGINAL
        if (t < 0.0) { /* Then we use it for new direction */
- #else
+ #else  /* Not POWELLOriginal but Brouard's */
        if (directest*t < 0.0) { /* Contradiction between both tests */
-                                 printf("directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del);
+         printf("directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del);
          printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);
          fprintf(ficlog,"directest= %.12lf (if directest<0 or t<0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del);
          fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);
        }
-       if (directest < 0.0) { /* Then we use it for new direction */
+       if (directest < 0.0) { /* Then we use (P0, Pn) for new direction Xi_n or Xi_iBig */
  #endif
  #ifdef DEBUGLINMIN
          printf("Before linmin in direction P%d-P0\n",n);
- Line 2915  void powell(double p[], double **xi, int
+ Line 4494  void powell(double p[], double **xi, int
            xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */
            xi[j][n]=xit[j];      /* and this nth direction by the by the average p_0 p_n */
          }
+ /* #else */
+ /*      for (i=1;i<=n-1;i++) {  */
+ /*        for (j=1;j<=n;j++) {  */
+ /*          xi[j][i]=xi[j][i+1]; /\* Standard method of conjugate directions, not Powell who changes the nth direction by p0 pn . *\/ */
+ /*        } */
+ /*      } */
+ /*      for (j=1;j<=n;j++) {  */
+ /*        xi[j][n]=xit[j];      /\* and this nth direction by the by the average p_0 p_n *\/ */
+ /*      } */
+ /*      /\* for (j=1;j<=n-1;j++) {  *\/ */
+ /*      /\*   xi[j][1]=xi[j][j+1]; /\\* Standard method of conjugate directions *\\/ *\/ */
+ /*      /\*   xi[j][n]=xit[j];      /\\* and this nth direction by the by the average p_0 p_n *\\/ *\/ */
+ /*      /\* } *\/ */
+ /* #endif */
  #ifdef LINMINORIGINAL
  #else
          for (j=1, flatd=0;j<=n;j++) {
- Line 2939  void powell(double p[], double **xi, int
+ Line 4533  void powell(double p[], double **xi, int
            free_vector(pt,1,n);
            return;
  #endif
-         }
+         }  /* endif(flatd >0) */
- #endif
+ #endif /* LINMINORIGINAL */
          printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
          fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
- Line 2955  void powell(double p[], double **xi, int
+ Line 4549  void powell(double p[], double **xi, int
          fprintf(ficlog,"\n");
  #endif
        } /* end of t or directest negative */
+       printf(" Directest is positive, P_n-P_0 does not increase the conjugacy. n=%d\n",n);
+       fprintf(ficlog," Directest is positive, P_n-P_0 does not increase the conjugacy. n=%d\n",n);
  #ifdef POWELLNOF3INFF1TEST
  #else
        } /* end if (fptt < fp)  */
- Line 3163  void powell(double p[], double **xi, int
+ Line 4759  void powell(double p[], double **xi, int
      first++;
    }
-   /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */
+   /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl,
+    * (int)age, (int)delaymax, (int)agefin, ncvloop,
+    * (int)age-(int)agefin); */
    free_vector(min,1,nlstate);
    free_vector(max,1,nlstate);
    free_vector(meandiff,1,nlstate);
- Line 3198  void powell(double p[], double **xi, int
+ Line 4796  void powell(double p[], double **xi, int
    /*  0.51326036147820708, 0.48673963852179264} */
    /* If we start from prlim again, prlim tends to a constant matrix */
-   int i, ii,j,k, k1;
+   int i, ii,j, k1;
    int first=0;
    double *min, *max, *meandiff, maxmax,sumnew=0.;
    /* double **matprod2(); */ /* test */
- Line 3465  double **pmij(double **ps, double *cov,
+ Line 5063  double **pmij(double **ps, double *cov,
    /* Computes the backward probability at age agefin, cov[2], and covariate combination 'ij'. In fact cov is already filled and x too.
     * Call to pmij(cov and x), call to cross prevalence, sums and inverses, left multiply, and returns in **ps as well as **bmij.
     */
-   int i, ii, j,k;
+   int ii, j;
-   double **out, **pmij();
+   double  **pmij();
    double sumnew=0.;
    double agefin;
    double k3=0.; /* constant of the w_x diagonal matrix (in order for B to sum to 1 even for death state) */
- Line 3680  double ***hpxij(double ***po, int nhstep
+ Line 5278  double ***hpxij(double ***po, int nhstep
       */
-   int i, j, d, h, k, k1;
+   int i, j, d, h, k1;
    double **out, cov[NCOVMAX+1];
    double **newm;
    double agexact;
-   double agebegin, ageend;
+   /*double agebegin, ageend;*/
    /* Hstepm could be zero and should return the unit matrix */
    for (i=1;i<=nlstate+ndeath;i++)
- Line 3861  double ***hbxij(double ***po, int nhstep
+ Line 5459  double ***hbxij(double ***po, int nhstep
       The addresss of po (p3mat allocated to the dimension of nhstepm) should be stored for output
    */
-   int i, j, d, h, k, k1;
+   int i, j, d, h, k1;
    double **out, cov[NCOVMAX+1], **bmij();
    double **newm, ***newmm;
    double agexact;
-   double agebegin, ageend;
+   /*double agebegin, ageend;*/
    double **oldm, **savm;
    newmm=po; /* To be saved */
- Line 4384  double func( double *x)
+ Line 5982  double func( double *x)
  double funcone( double *x)
  {
    /* Same as func but slower because of a lot of printf and if */
-   int i, ii, j, k, mi, d, kk, kv=0, kf=0;
+   int i, ii, j, k, mi, d, kv=0, kf=0;
    int ioffset=0;
    int ipos=0,iposold=0,ncovv=0;
- Line 4553  double funcone( double *x)
+ Line 6151  double funcone( double *x)
          * 3 ncovta=15    +age*V3*V2+age*V2+agev3+ageV4 +age*V6 + age*V7 + age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4
          * 3 TvarAVVA[1]@15= itva 3 2    2      3    4        6       7        6 3         7 3         6 4         7 4
          * 3 ncovta             1 2      3      4    5        6       7        8 9       10 11       12 13        14 15
-         * TvarAVVAind[1]@15= V3 is in k=2 1 1  2    3        4       5        4,2         5,2,      4,3           5 3}TvarVVAind[]
+         *?TvarAVVAind[1]@15= V3 is in k=2 1 1  2    3        4       5        4,2         5,2,      4,3           5 3}TvarVVAind[]
          * TvarAVVAind[1]@15= V3 is in k=6 6 12  13   14      15      16       18 18       19,19,     20,20        21,21}TvarVVAind[]
          * 3 ncovvta=10     +age*V6 + age*V7 + age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4
          * 3 we want to compute =cotvar[mw[mi][i]][TvarVVA[ncovva]][i] at position TvarVVAind[ncovva]
- Line 4568  double funcone( double *x)
+ Line 6166  double funcone( double *x)
          *                     6, 8, 9, 10, 11}
          * TvarFind[itv]                        0      0       0
          * FixedV[itv]                          1      1       1  0      1 0       1 0       1 0       0
+         *? FixedV[itv]                          1      1       1  0      1 0       1 0       1 0      1 0     1 0
          * Tvar[TvarFind[ncovf]]=[1]=2 [2]=3 [4]=4
          * Tvar[TvarFind[itv]]                [0]=?      ?ncovv 1 à ncovvt]
          *   Not a fixed cotvar[mw][itv][i]     6       7      6  2      7, 2,     6, 3,     7, 3,     6, 4,     7, 4}
- Line 4579  double funcone( double *x)
+ Line 6178  double funcone( double *x)
          ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
          /* if(TvarFind[itv]==0){ /\* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv *\/ */
          if(FixedV[itv]!=0){ /* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv */
+           /* printf("DEBUG ncovv=%d, Varying TvarVV[ncovv]=%d\n",ncovv, TvarVV[ncovv]); */
            cotvarv=cotvar[mw[mi][i]][TvarVV[ncovv]][i];  /* because cotvar starts now at first ncovcol+nqv+ntv+nqtv (1 to nqtv) */
+           /* printf("DEBUG Varying cov[ioffset+ipos=%d]=%g \n",ioffset+ipos,cotvarv); */
          }else{ /* fixed covariate */
            /* cotvarv=covar[Tvar[TvarFind[itv]]][i];  /\* Error: TvarFind gives the name, that is the true column of fixed covariates, but Tvar of the model *\/ */
+           /* printf("DEBUG ncovv=%d, Fixed TvarVV[ncovv]=%d\n",ncovv, TvarVV[ncovv]); */
            cotvarv=covar[itv][i];  /* Good: In V6*V3, 3 is fixed at position of the data */
+           /* printf("DEBUG Fixed cov[ioffset+ipos=%d]=%g \n",ioffset+ipos,cotvarv); */
          }
          if(ipos!=iposold){ /* Not a product or first of a product */
            cotvarvold=cotvarv;
- Line 4591  double funcone( double *x)
+ Line 6194  double funcone( double *x)
          }
          iposold=ipos;
          cov[ioffset+ipos]=cotvarv;
+         /* printf("DEBUG Product cov[ioffset+ipos=%d] \n",ioffset+ipos); */
          /* For products */
        }
        /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates single *\/ */
- Line 4928  void likelione(FILE *ficres,double p[],
+ Line 6532  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,k, jk, jkk=0, iter=0;
+   int i,j,  jkk=0, iter=0;
    double **xi;
-   double fret;
+   /*double fret;*/
-   double fretone; /* Only one call to likelihood */
+   /*double fretone;*/ /* Only one call to likelihood */
    /*  char filerespow[FILENAMELENGTH];*/
+   /*double * p1;*/ /* Shifted parameters from 0 instead of 1 */
  #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 (i=1;i<=npar;i++)  /* Starting with canonical directions j=1,n xi[i=1,n][j] */
      for (j=1;j<=npar;j++)
        xi[i][j]=(i==j ? 1.0 : 0.0);
-   printf("Powell\n");  fprintf(ficlog,"Powell\n");
+   printf("Powell-prax\n");  fprintf(ficlog,"Powell-prax\n");
    strcpy(filerespow,"POW_");
    strcat(filerespow,fileres);
    if((ficrespow=fopen(filerespow,"w"))==NULL) {
- Line 5013  void mlikeli(FILE *ficres,double p[], in
+ Line 6618  void mlikeli(FILE *ficres,double p[], in
    }
    powell(p,xi,npar,ftol,&iter,&fret,flatdir,func);
  #else  /* FLATSUP */
-   powell(p,xi,npar,ftol,&iter,&fret,func);
+ /*  powell(p,xi,npar,ftol,&iter,&fret,func);*/
+ /*   praxis ( t0, h0, n, prin, x, beale_f ); */
+  int prin=4;
+   /* double h0=0.25; */
+   /* double macheps; */
+   /* double fmin; */
+   macheps=pow(16.0,-13.0);
+ /* #include "praxis.h" */
+   /* Be careful that praxis start at x[0] and powell start at p[1] */
+    /* praxis ( ftol, h0, npar, prin, p, func ); */
+ /* p1= (p+1); */ /*  p *(p+1)@8 and p *(p1)@8 are equal p1[0]=p[1] */
+ printf("Praxis Gegenfurtner \n");
+ fprintf(ficlog, "Praxis  Gegenfurtner\n");fflush(ficlog);
+ /* praxis ( ftol, h0, npar, prin, p1, func ); */
+   /* fmin = praxis(1.e-5,macheps, h, n, prin, x, func); */
+   ffmin = praxis(ftol,macheps, h0, npar, prin, p, func);
+ printf("End Praxis\n");
  #endif  /* FLATSUP */
  #ifdef LINMINORIGINAL
- Line 5275  double hessij( double x[], double **hess
+ Line 6896  double hessij( double x[], double **hess
        kmax=kmax+10;
      }
      if(kmax >=10 || firstime ==1){
+       /* What are the thetai and thetaj? thetai/ncovmodel thetai=(thetai-thetai%ncovmodel)/ncovmodel +thetai%ncovmodel=(line,pos)  */
        printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol);
        fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol);
        printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e  res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
- Line 6143  void prevalence(double ***probs, double
+ Line 7765  void prevalence(double ***probs, double
    int i, m, jk, j1, bool, z1,j, iv;
    int mi; /* Effective wave */
    int iage;
-   double agebegin, ageend;
+   double agebegin; /*, ageend;*/
    double **prop;
    double posprop;
- Line 6382  void  concatwav(int wav[], int **dh, int
+ Line 8004  void  concatwav(int wav[], int **dh, int
              if(j==0) j=1;  /* Survives at least one month after exam */
              else if(j<0){
                nberr++;
-               printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+               printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld (around line %d) who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
                j=1; /* Temporary Dangerous patch */
                printf("   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
-               fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+               fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld (around line %d) who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
                fprintf(ficlog,"   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
              }
              k=k+1;
- Line 6419  void  concatwav(int wav[], int **dh, int
+ Line 8041  void  concatwav(int wav[], int **dh, int
            /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/
            if(j<0){
              nberr++;
-             printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+             printf("Error! Negative delay (%d) between waves %d and %d of individual %ld (around line %d) who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
-             fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+             fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld (around line %d) who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
            }
            sum=sum+j;
          }
- Line 6988  void  concatwav(int wav[], int **dh, int
+ Line 8610  void  concatwav(int wav[], int **dh, int
  /************ Variance ******************/
   void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[], int nres)
   {
-    /** Variance of health expectancies
+    /** Computes the matrix of variance covariance of health expectancies e.j= sum_i w_i e_ij where w_i depends of popbased,
+     * either cross-sectional or implied.
+     * return vareij[i][j][(int)age]=cov(e.i,e.j)=sum_h sum_k trgrad(h_p.i) V(theta) grad(k_p.k) Equation 20
      *  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);
      * double **newm;
      * int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav)
- Line 7005  void  concatwav(int wav[], int **dh, int
+ Line 8629  void  concatwav(int wav[], int **dh, int
     double ***gradg, ***trgradg; /**< for var eij */
     double **gradgp, **trgradgp; /**< for var p point j */
     double *gpp, *gmp; /**< for var p point j */
-    double **varppt; /**< for var p point j nlstate to nlstate+ndeath */
+    double **varppt; /**< for var p.3 p.death nlstate+1 to nlstate+ndeath */
     double ***p3mat;
     double age,agelim, hf;
     /* double ***mobaverage; */
- Line 7073  void  concatwav(int wav[], int **dh, int
+ Line 8697  void  concatwav(int wav[], int **dh, int
     fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");
     fprintf(fichtm,"\n<br>%s  <br>\n",digitp);
-    varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
+    varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); /* In fact, currently a double */
     pstamp(ficresvij);
     fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n#  (weighted average of eij where weights are ");
     if(popbased==1)
- Line 7142  void  concatwav(int wav[], int **dh, int
+ Line 8766  void  concatwav(int wav[], int **dh, int
               prlim[i][i]=mobaverage[(int)age][i][ij];
           }
         }
-        /**< Computes the shifted transition matrix \f$ {}{h}_p^{ij}x\f$ at horizon h.
+        /**< Computes the shifted plus (gp) transition matrix \f$ {}{h}_p^{ij}x\f$ at horizon h.
          */
         hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres);  /* Returns p3mat[i][j][h] for h=0 to nhstepm */
         /**< And for each alive state j, sums over i \f$ w^i_x {}{h}_p^{ij}x\f$, which are the probability
- Line 7151  void  concatwav(int wav[], int **dh, int
+ Line 8775  void  concatwav(int wav[], int **dh, int
         for(j=1; j<= nlstate; j++){
           for(h=0; h<=nhstepm; h++){
             for(i=1, gp[h][j]=0.;i<=nlstate;i++)
-              gp[h][j] += prlim[i][i]*p3mat[i][j][h];
+              gp[h][j] += prlim[i][i]*p3mat[i][j][h]; /* gp[h][j]= w_i h_pij */
           }
         }
         /* Next for computing shifted+ probability of death (h=1 means
            computed over hstepm matrices product = hstepm*stepm months)
            as a weighted average of prlim(i) * p(i,j) p.3=w1*p13 + w2*p23 .
         */
-        for(j=nlstate+1;j<=nlstate+ndeath;j++){
+        for(j=nlstate+1;j<=nlstate+ndeath;j++){ /* Currently only once for theta plus  p.3(age) Sum_i wi pi3*/
           for(i=1,gpp[j]=0.; i<= nlstate; i++)
             gpp[j] += prlim[i][i]*p3mat[i][j][1];
         }
- Line 7180  void  concatwav(int wav[], int **dh, int
+ Line 8804  void  concatwav(int wav[], int **dh, int
           }
         }
-        hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres);
+        hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres);  /* Still minus */
-        for(j=1; j<= nlstate; j++){  /* Sum of wi * eij = e.j */
+        for(j=1; j<= nlstate; j++){  /* gm[h][j]= Sum_i of wi * pij =  h_p.j */
           for(h=0; h<=nhstepm; h++){
             for(i=1, gm[h][j]=0.;i<=nlstate;i++)
               gm[h][j] += prlim[i][i]*p3mat[i][j][h];
- Line 7190  void  concatwav(int wav[], int **dh, int
+ Line 8814  void  concatwav(int wav[], int **dh, int
         }
         /* This for computing probability of death (h=1 means
            computed over hstepm matrices product = hstepm*stepm months)
-           as a weighted average of prlim.
+           as a weighted average of prlim. j is death. gmp[3]=sum_i w_i*p_i3=p.3 minus theta
         */
-        for(j=nlstate+1;j<=nlstate+ndeath;j++){
+        for(j=nlstate+1;j<=nlstate+ndeath;j++){  /* Currently only once theta_minus  p.3=Sum_i wi pi3*/
           for(i=1,gmp[j]=0.; i<= nlstate; i++)
             gmp[j] += prlim[i][i]*p3mat[i][j][1];
         }
         /* end shifting computations */
-        /**< Computing gradient matrix at horizon h
+        /**< Computing gradient of p.j matrix at horizon h and still for one parameter of vector theta
+         * equation 31 and 32
          */
-        for(j=1; j<= nlstate; j++) /* vareij */
+        for(j=1; j<= nlstate; j++) /* computes grad p.j(x, over each  h) where p.j is Sum_i w_i*pij(x over h)
+                                   * equation 24 */
           for(h=0; h<=nhstepm; h++){
             gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
           }
-        /**< Gradient of overall mortality p.3 (or p.j)
+        /**< Gradient of overall mortality p.3 (or p.death)
          */
-        for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu mortality from j */
+        for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* computes grad of p.3 from wi+pi3 grad p.3 (theta) */
           gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta];
         }
       } /* End theta */
-      /* We got the gradient matrix for each theta and state j */
+      /* We got the gradient matrix for each theta and each state j of gradg(h]theta][j)=grad(_hp.j(theta) */
-      trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); /* veij */
+      trgradg =ma3x(0,nhstepm,1,nlstate,1,npar);
-      for(h=0; h<=nhstepm; h++) /* veij */
+      for(h=0; h<=nhstepm; h++) /* veij */ /* computes the transposed of grad  (_hp.j(theta)*/
         for(j=1; j<=nlstate;j++)
           for(theta=1; theta <=npar; theta++)
             trgradg[h][j][theta]=gradg[h][theta][j];
-      for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */
+      for(j=nlstate+1; j<=nlstate+ndeath;j++) /* computes transposed of grad p.3 (theta)*/
         for(theta=1; theta <=npar; theta++)
           trgradgp[j][theta]=gradgp[theta][j];
       /**< as well as its transposed matrix
- Line 7232  void  concatwav(int wav[], int **dh, int
+ Line 8858  void  concatwav(int wav[], int **dh, int
           vareij[i][j][(int)age] =0.;
       /* Computing trgradg by matcov by gradg at age and summing over h
-       * and k (nhstepm) formula 15 of article
+       * and k (nhstepm) formula 32 of article
-       * Lievre-Brouard-Heathcote
+       * Lievre-Brouard-Heathcote so that for each j, computes the cov(e.j,e.k) (formula 31).
+       * for given h and k computes trgradg[h](i,j) matcov (theta) gradg(k)(i,j) into vareij[i][j] which is
+       cov(e.i,e.j) and sums on h and k
+       * including the covariances.
        */
       for(h=0;h<=nhstepm;h++){
- Line 7242  void  concatwav(int wav[], int **dh, int
+ Line 8871  void  concatwav(int wav[], int **dh, int
           matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]);
           for(i=1;i<=nlstate;i++)
             for(j=1;j<=nlstate;j++)
-              vareij[i][j][(int)age] += doldm[i][j]*hf*hf;
+              vareij[i][j][(int)age] += doldm[i][j]*hf*hf; /* This is vareij=sum_h sum_k trgrad(h_pij) V(theta) grad(k_pij)
+                                                              including the covariances of e.j */
         }
       }
-      /* pptj is p.3 or p.j = trgradgp by cov by gradgp, variance of
+      /* Mortality: pptj is p.3 or p.death = trgradgp by cov by gradgp, variance of
-       * p.j overall mortality formula 49 but computed directly because
+       * p.3=1-p..=1-sum i p.i  overall mortality computed directly because
        * we compute the grad (wix pijx) instead of grad (pijx),even if
        * wix is independent of theta.
        */
       matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov);
       matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp);
       for(j=nlstate+1;j<=nlstate+ndeath;j++)
         for(i=nlstate+1;i<=nlstate+ndeath;i++)
-          varppt[j][i]=doldmp[j][i];
+          varppt[j][i]=doldmp[j][i];  /* This is the variance of p.3 */
       /* end ppptj */
       /*  x centered again */
- Line 7278  void  concatwav(int wav[], int **dh, int
+ Line 8908  void  concatwav(int wav[], int **dh, int
       hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij, nres);
       for(j=nlstate+1;j<=nlstate+ndeath;j++){
         for(i=1,gmp[j]=0.;i<= nlstate; i++)
-          gmp[j] += prlim[i][i]*p3mat[i][j][1];
+          gmp[j] += prlim[i][i]*p3mat[i][j][1]; /* gmp[j] is p.3 */
       }
       /* end probability of death */
       fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij);
       for(j=nlstate+1; j<=(nlstate+ndeath);j++){
-        fprintf(ficresprobmorprev," %11.3e %11.3e",gmp[j], sqrt(varppt[j][j]));
+        fprintf(ficresprobmorprev," %11.3e %11.3e",gmp[j], sqrt(varppt[j][j]));/* p.3 (STD p.3) */
         for(i=1; i<=nlstate;i++){
-          fprintf(ficresprobmorprev," %11.3e %11.3e ",prlim[i][i],p3mat[i][j][1]);
+          fprintf(ficresprobmorprev," %11.3e %11.3e ",prlim[i][i],p3mat[i][j][1]); /* wi, pi3 */
         }
       }
       fprintf(ficresprobmorprev,"\n");
- Line 7612  void varprob(char optionfilefiname[], do
+ Line 9242  void varprob(char optionfilefiname[], do
     double ***varpij;
     strcpy(fileresprob,"PROB_");
-    strcat(fileresprob,fileres);
+    strcat(fileresprob,fileresu);
     if((ficresprob=fopen(fileresprob,"w"))==NULL) {
       printf("Problem with resultfile: %s\n", fileresprob);
       fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob);
- Line 8030  void printinghtml(char fileresu[], char
+ Line 9660  void printinghtml(char fileresu[], char
                    int popforecast, int mobilav, int prevfcast, int mobilavproj, int prevbcast, int estepm , \
                    double jprev1, double mprev1,double anprev1, double dateprev1, double dateprojd, double dateback1, \
                    double jprev2, double mprev2,double anprev2, double dateprev2, double dateprojf, double dateback2){
-   int jj1, k1, i1, cpt, k4, nres;
+   int jj1, k1, cpt, nres;
    /* In fact some results are already printed in fichtm which is open */
     fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \
     <li><a href='#secondorder'>Result files (second order (variance)</a>\n \
- Line 8167  divided by h: <sub>h</sub>P<sub>ij</sub>
+ Line 9797  divided by h: <sub>h</sub>P<sub>ij</sub>
  <img src=\"%s_%d-3-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);
       /* Survival functions (period) in state j */
       for(cpt=1; cpt<=nlstate;cpt++){
-        fprintf(fichtm,"<br>\n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
+        fprintf(fichtm,"<br>\n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages. Mean times spent in state (or Life Expectancy or Health Expectancy etc.) are the areas under each curve. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
         fprintf(fichtm," (data from text file  <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));
         fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
       }
       /* State specific survival functions (period) */
       for(cpt=1; cpt<=nlstate;cpt++){
         fprintf(fichtm,"<br>\n- Survival functions in state %d and in any other live state (total).\
-  And probability to be observed in various states (up to %d) being in state %d at different ages.       \
+  And probability to be observed in various states (up to %d) being in state %d at different ages.  Mean times spent in state (or Life Expectancy or Health Expectancy etc.) are the areas under each curve. \
   <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> ", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
         fprintf(fichtm," (data from text file  <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));
         fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
       }
       /* Period (forward stable) prevalence in each health state */
       for(cpt=1; cpt<=nlstate;cpt++){
-        fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability for a person being in state (1 to %d) at different ages, to be in state %d some years after. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
+        fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability for a person being in state (1 to %d) at different ages, to be alive in state %d some years after. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
         fprintf(fichtm," (data from text file  <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));
        fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">" ,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
       }
- Line 8206  divided by h: <sub>h</sub>P<sub>ij</sub>
+ Line 9836  divided by h: <sub>h</sub>P<sub>ij</sub>
        /* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */
         for(cpt=1; cpt<=nlstate;cpt++){
           fprintf(fichtm,"<br>\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), \
-  from year %.1f up to year %.1f (probably close to stable [mixed] back prevalence in state %d (randomness in cross-sectional prevalence is not taken into \
+  from year %.1f up to year %.1f (probably close to stable [mixed] back prevalence in state %d). Randomness in cross-sectional prevalence is not taken into \
-  account but can visually be appreciated). Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) \
+  account but can visually be appreciated. Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) \
  with weights corresponding to observed prevalence at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a>", dateprev1, dateprev2, mobilavproj, dateback1, dateback2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
           fprintf(fichtm," (data from text file  <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"FB_"),subdirf2(optionfilefiname,"FB_"));
           fprintf(fichtm," <img src=\"%s_%d-%d-%d.svg\">", subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
- Line 8345  prevalence (with 95%% confidence interva
+ Line 9975  prevalence (with 95%% confidence interva
         fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"V_"), cpt,k1,nres);
       }
       fprintf(fichtm,"\n<br>- Total life expectancy by age and \
- health expectancies in each live states (1 to %d). If popbased=1 the smooth (due to the model) \
+ health expectancies in each live state (1 to %d) with confidence intervals \
+ on left y-scale as well as proportions of time spent in each live state \
+ (with confidence intervals) on right y-scale 0 to 100%%.\
+  If popbased=1 the smooth (due to the model)                            \
  true period expectancies (those weighted with period prevalences are also\
   drawn in addition to the population based expectancies computed using\
   observed and cahotic prevalences:  <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a>",nlstate, subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres);
- Line 8360  true period expectancies (those weighted
+ Line 9993  true period expectancies (those weighted
  /******************* Gnuplot file **************/
  void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double bage, double fage , int prevfcast, int prevbcast, char pathc[], double p[], int offyear, int offbyear){
-   char dirfileres[132],optfileres[132];
+   char dirfileres[256],optfileres[256];
-   char gplotcondition[132], gplotlabel[132];
+   char gplotcondition[256], gplotlabel[256];
    int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,k4=0,kf=0,kvar=0,kk=0,ipos=0,iposold=0,ij=0, ijp=0, l=0;
-   int lv=0, vlv=0, kl=0;
+   /* int lv=0, vlv=0, kl=0; */
+   int lv=0, kl=0;
+   double vlv=0;
    int ng=0;
    int vpopbased;
    int ioffset; /* variable offset for columns */
- Line 8432  void printinggnuplot(char fileresu[], ch
+ Line 10067  void printinggnuplot(char fileresu[], ch
      kvar=Tvar[TvarFind[kf]]; /* variable name */
      /* k=18+Tvar[TvarFind[kf]];/\*offset because there are 18 columns in the ILK_ file but could be placed else where *\/ */
      /* k=18+kf;/\*offset because there are 18 columns in the ILK_ file *\/ */
-     k=19+kf;/*offset because there are 19 columns in the ILK_ file */
+     /* k=19+kf;/\*offset because there are 19 columns in the ILK_ file *\/ */
+     k=16+nlstate+kf;/*offset because there are 19 columns in the ILK_ file, first cov Vn on col 21 with 4 living states */
      for (i=1; i<= nlstate ; i ++) {
        fprintf(ficgp,"\nset out \"%s-p%dj-%d.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i,kvar);
        fprintf(ficgp,"unset log;\n# For each simple dummy covariate of the model \n plot  \"%s\"",subdirf(fileresilk));
- Line 8698  void printinggnuplot(char fileresu[], ch
+ Line 10334  void printinggnuplot(char fileresu[], ch
        for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
          fprintf(ficgp,"\nset label \"popbased %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",vpopbased,gplotlabel);
          if(vpopbased==0){
-           fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
+           fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nunset ytics; unset y2tics; set ytics nomirror; set y2tics 0,10,100;set y2range [0:100];\nplot [%.f:%.f] ",ageminpar,fage);
          }else
            fprintf(ficgp,"\nreplot ");
-         for (i=1; i<= nlstate+1 ; i ++) {
+         for (i=1; i<= nlstate+1 ; i ++) { /* For state i-1=0 is LE, while i-1=1 to nlstate are origin state */
            k=2*i;
-           fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1, vpopbased);
+           fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1, vpopbased); /* for fixed variables age, popbased, mobilav */
-           for (j=1; j<= nlstate+1 ; j ++) {
+           for (j=1; j<= nlstate+1 ; j ++) { /* e.. e.1 e.2 again j-1 is the state of end, wlim_i eij*/
-             if (j==i) fprintf(ficgp," %%lf (%%lf)");
+             if (j==i) fprintf(ficgp," %%lf (%%lf)"); /* We want to read e.. i=1,j=1, e.1 i=2,j=2, e.2 i=3,j=3 */
-             else fprintf(ficgp," %%*lf (%%*lf)");
+             else fprintf(ficgp," %%*lf (%%*lf)");  /* skipping that field with a star */
            }
            if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i);
-           else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1);
+           else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1); /* state=i-1=1 to nlstate  */
            fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1,vpopbased);
            for (j=1; j<= nlstate+1 ; j ++) {
              if (j==i) fprintf(ficgp," %%lf (%%lf)");
- Line 8721  void printinggnuplot(char fileresu[], ch
+ Line 10357  void printinggnuplot(char fileresu[], ch
              if (j==i) fprintf(ficgp," %%lf (%%lf)");
              else fprintf(ficgp," %%*lf (%%*lf)");
            }
-           if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
+           if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0,\\\n"); /* ,\\\n added for th percentage graphs */
            else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
          } /* state */
+         /* again for the percentag spent in state i-1=1 to i-1=nlstate */
+         for (i=2; i<= nlstate+1 ; i ++) { /* For state i-1=0 is LE, while i-1=1 to nlstate are origin state */
+           k=2*i;
+           fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d &&  ($4)<=1 && ($4)>=0 ?($4)*100. : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1, vpopbased); /* for fixed variables age, popbased, mobilav */
+           for (j=1; j<= nlstate ; j ++)
+             fprintf(ficgp," %%*lf (%%*lf)"); /* Skipping TLE and LE to read %LE only */
+           for (j=1; j<= nlstate+1 ; j ++) { /* e.. e.1 e.2 again j-1 is the state of end, wlim_i eij*/
+             if (j==i) fprintf(ficgp," %%lf (%%lf)"); /* We want to read e.. i=1,j=1, e.1 i=2,j=2, e.2 i=3,j=3 */
+             else fprintf(ficgp," %%*lf (%%*lf)");  /* skipping that field with a star */
+           }
+           if (i== 1) fprintf(ficgp,"\" t\"%%TLE\" w l lt %d axis x1y2, \\\n",i); /* Not used */
+           else fprintf(ficgp,"\" t\"%%LE in state (%d)\" w l lw 2 lt %d axis x1y2, \\\n",i-1,i+1); /* state=i-1=1 to nlstate  */
+           fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && ($4-$5*2)<=1 && ($4-$5*2)>=0? ($4-$5*2)*100. : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1,vpopbased);
+           for (j=1; j<= nlstate ; j ++)
+             fprintf(ficgp," %%*lf (%%*lf)"); /* Skipping TLE and LE to read %LE only */
+           for (j=1; j<= nlstate+1 ; j ++) {
+             if (j==i) fprintf(ficgp," %%lf (%%lf)");
+             else fprintf(ficgp," %%*lf (%%*lf)");
+           }
+           fprintf(ficgp,"\" t\"\" w l lt 0 axis x1y2,");
+           fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && ($4+$5*2)<=1 && ($4+$5*2)>=0 ? ($4+$5*2)*100. : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1,vpopbased);
+           for (j=1; j<= nlstate ; j ++)
+             fprintf(ficgp," %%*lf (%%*lf)"); /* Skipping TLE and LE to read %LE only */
+           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 lt 0 axis x1y2");
+           else fprintf(ficgp,"\" t\"\" w l lt 0 axis x1y2,\\\n");
+         } /* state for percent */
        } /* vpopbased */
        fprintf(ficgp,"\nset out;set out \"%s_%d-%d.svg\"; replot; set out; unset label;\n",subdirf2(optionfilefiname,"E_"),k1,nres); /* Buggy gnuplot */
      } /* end nres */
- Line 9130  set ter svg size 640, 480\nunset log y\n
+ Line 10796  set ter svg size 640, 480\nunset log y\n
                /* vlv= nbcode[Tvaraff[k]][lv]; /\* Value of the modality of Tvaraff[k] *\/ */
                /* vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])]; */
                kl++;
+               /* Problem with quantitative variables TinvDoQresult[nres] */
                /* sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]); */
-               sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,lv, kl+1, vlv );
+               sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%lg " ,kl,lv, kl+1, vlv );/* Solved but quantitative must be shifted */
                kl++;
                if(k <cptcovs && cptcovs>1)
                  sprintf(gplotcondition+strlen(gplotcondition)," && ");
- Line 9145  set ter svg size 640, 480\nunset log y\n
+ Line 10812  set ter svg size 640, 480\nunset log y\n
                fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0):%d t 'p.%d' with line lc variable", gplotcondition, \
                        ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,iyearc, cpt );
                fprintf(ficgp,",\\\n '' ");
-               fprintf(ficgp," u %d:(",iagec);
+               fprintf(ficgp," u %d:(",iagec); /* Below iyearc should be increades if quantitative variable in the reult line */
+               /* $7==6 && $8==2.47 ) && (($9-$10) == 1953 ) ? $12/(1.-$24) : 1/0):7 with labels center not */
+               /* but was  && $7==6 && $8==2 ) && (($7-$8) == 1953 ) ? $12/(1.-$24) : 1/0):7 with labels center not */
                fprintf(ficgp,"%s && (($%d-$%d) == %d ) ? $%d/(1.-$%d) : 1/0):%d with labels center not ", gplotcondition, \
                        iyearc, iagec, offyear,                           \
                        ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate, iyearc );
- Line 9893  void prevforecast(char fileres[], double
+ Line 11562  void prevforecast(char fileres[], double
    */
    /* double anprojd, mprojd, jprojd; */
    /* double anprojf, mprojf, jprojf; */
-   int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
+   int yearp, stepsize, hstepm, nhstepm, j, k, i, h,  nres=0;
    double agec; /* generic age */
-   double agelim, ppij, yp,yp1,yp2;
+   double agelim, ppij;
-   double *popeffectif,*popcount;
+   /*double *popcount;*/
    double ***p3mat;
    /* double ***mobaverage; */
    char fileresf[FILENAMELENGTH];
- Line 10044  void prevforecast(char fileres[], double
+ Line 11713  void prevforecast(char fileres[], double
       anback2 year of end of backprojection (same day and month as back1).
       prevacurrent and prev are prevalences.
    */
-   int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
+   int yearp, stepsize, hstepm, nhstepm, j, k,  i, h, nres=0;
    double agec; /* generic age */
-   double agelim, ppij, ppi, yp,yp1,yp2; /* ,jintmean,mintmean,aintmean;*/
+   double agelim, ppij, ppi; /* ,jintmean,mintmean,aintmean;*/
-   double *popeffectif,*popcount;
+   /*double *popcount;*/
    double ***p3mat;
    /* double ***mobaverage; */
    char fileresfb[FILENAMELENGTH];
- Line 10646  double gompertz(double x[])
+ Line 12315  double gompertz(double x[])
           A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)));
         } else if (cens[i] == 0){
          A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)))
-           +log(x[1]/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM);
+           +log(fabs(x[1])/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM);
+         /* +log(x[1]/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM); */  /* To be seen */
        } else
           printf("Gompertz cens[%d] neither 1 nor 0\n",i);
        /*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */
- Line 10730  void printinggnuplotmort(char fileresu[]
+ Line 12400  void printinggnuplotmort(char fileresu[]
    char dirfileres[132],optfileres[132];
-   int ng;
+   /*int ng;*/
    /*#ifdef windows */
- Line 10754  int readdata(char datafile[], int firsto
+ Line 12424  int readdata(char datafile[], int firsto
    /*-------- data file ----------*/
    FILE *fic;
    char dummy[]="                         ";
-   int i=0, j=0, n=0, iv=0, v;
+   int i = 0, j = 0, n = 0, iv = 0;/* , v;*/
    int lstra;
    int linei, month, year,iout;
    int noffset=0; /* This is the offset if BOM data file */
- Line 11371  int decodemodel( char model[], int lasto
+ Line 13041  int decodemodel( char model[], int lasto
          */
  /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1, Tage[1]=2 */
  {
-   int i, j, k, ks, v;
+   int i, j, k, ks;/* , v;*/
    int n,m;
    int  j1, k1, k11, k12, k2, k3, k4;
    char modelsav[300];
- Line 12799  int hPijx(double *p, int bage, int fage)
+ Line 14469  int hPijx(double *p, int bage, int fage)
    int agelim;
    int hstepm;
    int nhstepm;
-   int h, i, i1, j, k, k4, nres=0;
+   int h, i, i1, j, k, nres=0;
    double agedeb;
    double ***p3mat;
- Line 13003  int main(int argc, char *argv[])
+ Line 14673  int main(int argc, char *argv[])
    double ageminpar=AGEOVERFLOW,agemin=AGEOVERFLOW, agemaxpar=-AGEOVERFLOW, agemax=-AGEOVERFLOW;
    double ageminout=-AGEOVERFLOW,agemaxout=AGEOVERFLOW; /* Smaller Age range redefined after movingaverage */
+   double stdpercent; /* for computing the std error of percent e.i: e.i/e.. */
    double fret;
    double dum=0.; /* Dummy variable */
-   double ***p3mat;
+   /* double*** p3mat;*/
    /* double ***mobaverage; */
    double wald;
- Line 13018  int main(int argc, char *argv[])
+ Line 14689  int main(int argc, char *argv[])
    char pathr[MAXLINE], pathimach[MAXLINE];
    char *tok, *val; /* pathtot */
    /* int firstobs=1, lastobs=10; /\* nobs = lastobs-firstobs declared globally ;*\/ */
-   int c,  h , cpt, c2;
+   int c, h; /* c2; */
    int jl=0;
    int i1, j1, jk, stepsize=0;
    int count=0;
- Line 13053  int main(int argc, char *argv[])
+ Line 14724  int main(int argc, char *argv[])
    double ***delti3; /* Scale */
    double *delti; /* Scale */
    double ***eij, ***vareij;
-   double **varpl; /* Variances of prevalence limits by age */
+   //double **varpl; /* Variances of prevalence limits by age */
    double *epj, vepp;
- Line 13111  int main(int argc, char *argv[])
+ Line 14782  int main(int argc, char *argv[])
    getcwd(pathcd, size);
  #endif
    syscompilerinfo(0);
-   printf("\nIMaCh version %s, %s\n%s",version, copyright, fullversion);
+   printf("\nIMaCh prax version %s, %s\n%s",version, copyright, fullversion);
    if(argc <=1){
      printf("\nEnter the parameter file name: ");
      if(!fgets(pathr,FILENAMELENGTH,stdin)){
- Line 14013  This file: <a href=\"%s\">%s</a></br>Tit
+ Line 15684  This file: <a href=\"%s\">%s</a></br>Tit
    /* Calculates basic frequencies. Computes observed prevalence at single age
                   and for any valid combination of covariates
       and prints on file fileres'p'. */
-   freqsummary(fileres, p, pstart, agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \
+   freqsummary(fileres, p, pstart, (double)agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \
                firstpass, lastpass,  stepm,  weightopt, model);
    fprintf(fichtm,"\n");
- Line 14104  Interval (in months) between two waves:
+ Line 15775  Interval (in months) between two waves:
  #ifdef GSL
      printf("GSL optimization\n");  fprintf(ficlog,"Powell\n");
  #else
-     printf("Powell\n");  fprintf(ficlog,"Powell\n");
+     printf("Powell-mort\n");  fprintf(ficlog,"Powell-mort\n");
  #endif
      strcpy(filerespow,"POW-MORT_");
      strcat(filerespow,fileresu);
- Line 14199  Interval (in months) between two waves:
+ Line 15870  Interval (in months) between two waves:
      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);
+ #ifdef LINMINORIGINAL
- #endif
+ #else /* LINMINORIGINAL */
+   flatdir=ivector(1,npar);
+   for (j=1;j<=npar;j++) flatdir[j]=0;
+ #endif /*LINMINORIGINAL */
+     /* powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); */
+   /* double h0=0.25; */
+   macheps=pow(16.0,-13.0);
+   printf("Praxis Gegenfurtner mle=%d\n",mle);
+   fprintf(ficlog, "Praxis  Gegenfurtner mle=%d\n", mle);fflush(ficlog);
+    /* ffmin = praxis(ftol,macheps, h0, npar, prin, p, gompertz); */
+   /* For the Gompertz we use only two parameters */
+   int _npar=2;
+    ffmin = praxis(ftol,macheps, h0, _npar, 4, p, gompertz);
+   printf("End Praxis\n");
      fclose(ficrespow);
+ #ifdef LINMINORIGINAL
+ #else
+       free_ivector(flatdir,1,npar);
+ #endif  /* LINMINORIGINAL*/
+ #endif /* POWELL */
      hesscov(matcov, hess, p, NDIM, delti, 1e-4, gompertz);
      for(i=1; i <=NDIM; i++)
        for(j=i+1;j<=NDIM;j++)
-                                 matcov[i][j]=matcov[j][i];
+         matcov[i][j]=matcov[j][i];
      printf("\nCovariance matrix\n ");
      fprintf(ficlog,"\nCovariance matrix\n ");
- Line 14336  Please run with mle=-1 to get a correct
+ Line 16025  Please run with mle=-1 to get a correct
        fprintf(ficlog,"  + age*age  ");
        fprintf(fichtm, "<th>+ age*age</th>");
      }
-     for(j=1;j <=ncovmodel-2;j++){
+     for(j=1;j <=ncovmodel-2-nagesqr;j++){
        if(Typevar[j]==0) {
          printf("  +      V%d  ",Tvar[j]);
          fprintf(ficres,"  +      V%d  ",Tvar[j]);
- Line 14407  Please run with mle=-1 to get a correct
+ Line 16096  Please run with mle=-1 to get a correct
          fprintf(ficlog,"  + age*age  ");
          fprintf(fichtm, "<th>+ age*age</th>");
        }
-       for(j=1;j <=ncovmodel-2;j++){
+       for(j=1;j <=ncovmodel-2-nagesqr;j++){
          if(Typevar[j]==0) {
            printf("  +      V%d  ",Tvar[j]);
            fprintf(fichtm, "<th>+ V%d</th>",Tvar[j]);
- Line 14661  Please run with mle=-1 to get a correct
+ Line 16350  Please run with mle=-1 to get a correct
      }
      /* Results */
-     /* Value of covariate in each resultine will be compututed (if product) and sorted according to model rank */
+     /* Value of covariate in each resultine will be computed (if product) and sorted according to model rank */
      /* It is precov[] because we need the varying age in order to compute the real cov[] of the model equation */
      precov=matrix(1,MAXRESULTLINESPONE,1,NCOVMAX+1);
      endishere=0;
- Line 15064  Please run with mle=-1 to get a correct
+ Line 16753  Please run with mle=-1 to get a correct
        /* */
        if(i1 != 1 && TKresult[nres]!= k) /* TKresult[nres] is the combination of this nres resultline. All the i1 combinations are not output */
          continue;
-       printf("\n# model %s \n#****** Result for:", model);  /* HERE model is empty */
+       printf("\n# model=1+age+%s \n#****** Result for:", model);  /* HERE model is empty */
-       fprintf(ficrest,"\n# model %s \n#****** Result for:", model);
+       fprintf(ficrest,"\n# model=1+age+%s \n#****** Result for:", model);
-       fprintf(ficlog,"\n# model %s \n#****** Result for:", model);
+       fprintf(ficlog,"\n# model=1+age+%s \n#****** Result for:", model);
        /* It might not be a good idea to mix dummies and quantitative */
        /* for(j=1;j<=cptcoveff;j++){ /\* j=resultpos. Could be a loop on cptcovs: number of single dummy covariate in the result line as well as in the model *\/ */
        for(j=1;j<=cptcovs;j++){ /* j=resultpos. Could be a loop on cptcovs: number of single covariate (dummy or quantitative) in the result line as well as in the model */
- Line 15174  Please run with mle=-1 to get a correct
+ Line 16863  Please run with mle=-1 to get a correct
        for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
          oldm=oldms;savm=savms; /* ZZ Segmentation fault */
          cptcod= 0; /* To be deleted */
-         printf("varevsij vpopbased=%d \n",vpopbased);
+         printf("varevsij vpopbased=%d popbased=%d \n",vpopbased,popbased);
-         fprintf(ficlog, "varevsij vpopbased=%d \n",vpopbased);
+         fprintf(ficlog, "varevsij vpopbased=%d popbased=%d \n",vpopbased,popbased);
+         /* Call to varevsij to get cov(e.i, e.j)= vareij[i][j][(int)age]=sum_h sum_k trgrad(h_p.i) V(theta) grad(k_p.k) Equation 20 */
+         /* Depending of popbased which changes the prevalences, either cross-sectional or period */
          varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart, nres); /* 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 ");
+         fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each state\n\
+ #  (these are 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);
+           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
-           fprintf(ficrest,"the age specific forward period (stable) prevalences in each health state \n");
+           fprintf(ficrest,"the age specific forward period (stable) prevalences in each state) \n");
+         fprintf(ficrest,"# with proportions of time spent in each state with standard error (on the right of the table.\n ");
          fprintf(ficrest,"# Age popbased mobilav e.. (std) "); /* Adding covariate values? */
          for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
+         for (i=1;i<=nlstate;i++) fprintf(ficrest," %% e.%d/e.. (std) ",i);
          fprintf(ficrest,"\n");
          /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */
          printf("Computing age specific forward period (stable) prevalences in each health state \n");
- Line 15209  Please run with mle=-1 to get a correct
+ Line 16903  Please run with mle=-1 to get a correct
                /*ZZZ  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
                /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */
              }
-             epj[nlstate+1] +=epj[j];
+             epj[nlstate+1] +=epj[j]; /* epp=sum_j epj = sum_j sum_i w_i e_ij */
            }
            /* printf(" age %4.0f \n",age); */
-           for(i=1, vepp=0.;i <=nlstate;i++)
+           for(i=1, vepp=0.;i <=nlstate;i++)  /* Variance of total life expectancy e.. */
              for(j=1;j <=nlstate;j++)
-               vepp += vareij[i][j][(int)age];
+               vepp += vareij[i][j][(int)age]; /* sum_i sum_j cov(e.i, e.j) = var(e..) */
            fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
+           /* vareij[i][j] is the covariance  cov(e.i, e.j) and vareij[j][j] is the variance  of e.j  */
            for(j=1;j <=nlstate;j++){
              fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
            }
+           /* And proportion of time spent in state j */
+           /* $$ E[r(X,Y)-E(r(X,Y))]^2=[\frac{1}{\mu_y} -\frac{\mu_x}{{\mu_y}^2}]' Var(X,Y)[\frac{1}{\mu_y} -\frac{\mu_x}{{\mu_y}^2}]$$ */
+           /* \frac{\mu_x^2}{\mu_y^2} ( \frac{\sigma^2_x}{\mu_x^2}-2\frac{\sigma_{xy}}{\mu_x\mu_y} +\frac{\sigma^2_y}{\mu_y^2}) */
+           /* \frac{e_{.i}^2}{e_{..}^2} ( \frac{\Var e_{.i}}{e_{.i}^2}-2\frac{\Var e_{.i} + \sum_{j\ne i} \Cov e_{.j},e_{.i}}{e_{.i}e_{..}} +\frac{\Var e_{..}}{e_{..}^2})*/
+           /*\mu_x = epj[j], \sigma^2_x = vareij[j][j][(int)age] and \mu_y=epj[nlstate+1], \sigma^2_y=vepp \sigmaxy= */
+           /* vareij[j][j][(int)age]/epj[nlstate+1]^2 + vepp/epj[nlstate+1]^4 */
+           for(j=1;j <=nlstate;j++){
+             /* fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt( vareij[j][j][(int)age]/epj[j]/epj[j] + vepp/epj[j]/epj[j]/epj[j]/epj[j] )); */
+             /* fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt( vareij[j][j][(int)age]/epj[j]/epj[j] + vepp/epj[j]/epj[j]/epj[j]/epj[j] )); */
+             for(i=1,stdpercent=0.;i<=nlstate;i++){ /* Computing cov(e..,e.j)=cov(sum_i e.i,e.j)=sum_i cov(e.i, e.j) */
+               stdpercent += vareij[i][j][(int)age];
+             }
+             stdpercent= epj[j]*epj[j]/epj[nlstate+1]/epj[nlstate+1]* (vareij[j][j][(int)age]/epj[j]/epj[j]-2.*stdpercent/epj[j]/epj[nlstate+1]+ vepp/epj[nlstate+1]/epj[nlstate+1]);
+             /* stdpercent= epj[j]*epj[j]/epj[nlstate+1]/epj[nlstate+1]*(vareij[j][j][(int)age]/epj[j]/epj[j] + vepp/epj[nlstate+1]/epj[nlstate+1]); */ /* Without covariance */
+             /* fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt( vareij[j][j][(int)age]/epj[nlstate+1]/epj[nlstate+1] + epj[j]*epj[j]*vepp/epj[nlstate+1]/epj[nlstate+1]/epj[nlstate+1]/epj[nlstate+1] )); */
+             fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt(stdpercent));
+           }
            fprintf(ficrest,"\n");
          }
        } /* End vpopbased */
- Line 15261  Please run with mle=-1 to get a correct
+ Line 16974  Please run with mle=-1 to get a correct
      free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);
    }  /* mle==-3 arrives here for freeing */
    /* endfree:*/
+   if(mle!=-3) free_matrix(precov, 1,MAXRESULTLINESPONE,1,NCOVMAX+1); /* Could be elsewhere ?*/
    free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
    free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
    free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);
- Line 15322  Please run with mle=-1 to get a correct
+ Line 17036  Please run with mle=-1 to get a correct
    free_ivector(TmodelInvind,1,NCOVMAX);
    free_ivector(TmodelInvQind,1,NCOVMAX);
-   free_matrix(precov, 1,MAXRESULTLINESPONE,1,NCOVMAX+1); /* Could be elsewhere ?*/
+   /* free_matrix(precov, 1,MAXRESULTLINESPONE,1,NCOVMAX+1); /\* Could be elsewhere ?*\/ */
    free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX);
    /* free_imatrix(codtab,1,100,1,10); */

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