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.350 and 1.359

-version 1.350, 2023/04/24 11:38:06
+version 1.359, 2024/04/24 21:21:17
  Line 1
  /* $Id$
    $State$
    $Log$
+   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 ***
+   Revision 1.352  2023/04/29 10:46:21  brouard
+   *** empty log message ***
+   Revision 1.351  2023/04/29 10:43:47  brouard
+   Summary: 099r45
    Revision 1.350  2023/04/24 11:38:06  brouard
    *** empty log message ***
- Line 1269  Important routines
+ Line 1296  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 1365  double gnuplotversion=GNUPLOTVERSION;
+ Line 1394  double gnuplotversion=GNUPLOTVERSION;
  /* $State$ */
  #include "version.h"
  char version[]=__IMACH_VERSION__;
- char copyright[]="January 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 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 fullversion[]="$Revision$ $Date$";
  char strstart[80];
  char optionfilext[10], optionfilefiname[FILENAMELENGTH];
- Line 1413  int *wav; /* Number of waves for this in
+ Line 1442  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 1603  int **nbcode, *Tvar; /**< model=V2 => Tv
+ Line 1633  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 1804  char *trimbb(char *out, char *in)
+ Line 1834  char *trimbb(char *out, char *in)
    return s;
  }
+ char *trimbtab(char *out, char *in)
+ { /* Trim  blanks or tabs in line but keeps first blanks if line starts with blanks */
+   char *s;
+   s=out;
+   while (*in != '\0'){
+     while( (*in == ' ' || *in == '\t')){ /* && *(in+1) != '\0'){*/
+       in++;
+     }
+     *out++ = *in++;
+   }
+   *out='\0';
+   return s;
+ }
  /* char *substrchaine(char *out, char *in, char *chain) */
  /* { */
  /*   /\* Substract chain 'chain' from 'in', return and output 'out' *\/ */
- Line 2586  void linmin(double p[], double xi[], int
+ Line 2630  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 );
+   /* 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 );
+ #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]);
+      /* r8vec_print ( n, x, "  X:" ); */
+    }
+    printf("\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=%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 */
+ #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 2637  void powell(double p[], double **xi, int
+ Line 4187  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 2708  void powell(double p[], double **xi, int
+ Line 4259  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 2718  void powell(double p[], double **xi, int
+ Line 4270  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 2803  void powell(double p[], double **xi, int
+ Line 4357  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 2828  void powell(double p[], double **xi, int
+ Line 4388  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 2853  void powell(double p[], double **xi, int
+ Line 4413  void powell(double p[], double **xi, int
        if (t < 0.0) { /* Then we use it for new direction */
  #else
        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);
- Line 2932  void powell(double p[], double **xi, int
+ Line 4492  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 3140  void powell(double p[], double **xi, int
+ Line 4702  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 3175  void powell(double p[], double **xi, int
+ Line 4739  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 3442  double **pmij(double **ps, double *cov,
+ Line 5006  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 3657  double ***hpxij(double ***po, int nhstep
+ Line 5221  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 3838  double ***hbxij(double ***po, int nhstep
+ Line 5402  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 4361  double func( double *x)
+ Line 5925  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 4530  double funcone( double *x)
+ Line 6094  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 4545  double funcone( double *x)
+ Line 6109  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 4556  double funcone( double *x)
+ Line 6121  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 4568  double funcone( double *x)
+ Line 6137  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 4905  void likelione(FILE *ficres,double p[],
+ Line 6475  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 4990  void mlikeli(FILE *ficres,double p[], in
+ Line 6561  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=1;
+   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); */
+   fmin = praxis(ftol,macheps, h0, npar, prin, p, func);
+ printf("End Praxis\n");
  #endif  /* FLATSUP */
  #ifdef LINMINORIGINAL
- Line 5252  double hessij( double x[], double **hess
+ Line 6839  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 6120  void prevalence(double ***probs, double
+ Line 7708  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 6359  void  concatwav(int wav[], int **dh, int
+ Line 7947  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 6396  void  concatwav(int wav[], int **dh, int
+ Line 7984  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 7589  void varprob(char optionfilefiname[], do
+ Line 9177  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 8007  void printinghtml(char fileresu[], char
+ Line 9595  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 8144  divided by h: <sub>h</sub>P<sub>ij</sub>
+ Line 9732  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 8183  divided by h: <sub>h</sub>P<sub>ij</sub>
+ Line 9771  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 8337  true period expectancies (those weighted
+ Line 9925  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 ng=0;
- Line 8409  void printinggnuplot(char fileresu[], ch
+ Line 9997  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 9095  set ter svg size 640, 480\nunset log y\n
+ Line 10684  set ter svg size 640, 480\nunset log y\n
              fprintf(ficgp," u %d:(",ioffset);
              kl=0;
              strcpy(gplotcondition,"(");
-             for (k=1; k<=cptcoveff; k++){    /* For each covariate writing the chain of conditions */
+             /* for (k=1; k<=cptcoveff; k++){    /\* For each covariate writing the chain of conditions *\/ */
                /* lv= decodtabm(k1,k,cptcoveff); /\* Should be the covariate value corresponding to combination k1 and covariate k *\/ */
-               lv=codtabm(k1,TnsdVar[Tvaraff[k]]);
+             for (k=1; k<=cptcovs; k++){    /* For each covariate k get corresponding value lv for combination k1 */
+               /* lv=codtabm(k1,TnsdVar[Tvaraff[k]]); */
+               lv=Tvresult[nres][k];
+               vlv=TinvDoQresult[nres][Tvresult[nres][k]];
                /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
                /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
                /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
                /* vlv= nbcode[Tvaraff[k]][lv]; /\* Value of the modality of Tvaraff[k] *\/ */
-               vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])];
+               /* vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])]; */
                kl++;
-               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,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]); */
+               sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,lv, kl+1, vlv );
                kl++;
-               if(k <cptcoveff && cptcoveff>1)
+               if(k <cptcovs && cptcovs>1)
                  sprintf(gplotcondition+strlen(gplotcondition)," && ");
              }
              strcpy(gplotcondition+strlen(gplotcondition),")");
- Line 9190  set ter svg size 640, 480\nunset log y\n
+ Line 10783  set ter svg size 640, 480\nunset log y\n
            }else{
              fprintf(ficgp,",\\\n '' ");
            }
-           if(cptcoveff ==0){ /* No covariate */
+           /* if(cptcoveff ==0){ /\* No covariate *\/ */
+           if(cptcovs ==0){ /* No covariate */
              ioffset=2; /* Age is in 2 */
              /*# yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/
              /*#   1       2   3   4   5  6    7  8   9   10  11  12  13  14  15  16  17  18 */
- Line 9302  set ter svg size 640, 480\nunset log y\n
+ Line 10896  set ter svg size 640, 480\nunset log y\n
      fprintf(ficgp,"#Number of graphics: first is logit, 2nd is probabilities, third is incidences per year\n");
      fprintf(ficgp,"#model=1+age+%s \n",model);
      fprintf(ficgp,"# Type of graphic ng=%d\n",ng);
-     fprintf(ficgp,"#   k1=1 to 2^%d=%d\n",cptcoveff,m);/* to be checked */
+     /* fprintf(ficgp,"#   k1=1 to 2^%d=%d\n",cptcoveff,m);/\* to be checked *\/ */
+     fprintf(ficgp,"#   k1=1 to 2^%d=%d\n",cptcovs,m);/* to be checked */
      /* for(k1=1; k1 <=m; k1++)  /\* For each combination of covariate *\/ */
      for(nres=1; nres <= nresult; nres++){ /* For each resultline */
       /* k1=nres; */
- Line 9864  void prevforecast(char fileres[], double
+ Line 11459  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 9920  void prevforecast(char fileres[], double
+ Line 11515  void prevforecast(char fileres[], double
    /* date2dmy(dateintmean,&jintmean,&mintmean,&aintmean); */
    /* date2dmy(dateprojd,&jprojd, &mprojd, &anprojd); */
    /* date2dmy(dateprojf,&jprojf, &mprojf, &anprojf); */
-   i1=pow(2,cptcoveff);
+   /* i1=pow(2,cptcoveff); */
-   if (cptcovn < 1){i1=1;}
+   /* if (cptcovn < 1){i1=1;} */
    fprintf(ficresf,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jintmean,mintmean,aintmean,dateintmean,dateprev1,dateprev2);
    fprintf(ficresf,"#****** Routine prevforecast **\n");
  /*            if (h==(int)(YEARM*yearp)){ */
-   for(nres=1; nres <= nresult; nres++) /* For each resultline */
+   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-     for(k=1; k<=i1;k++){ /* We want to find the combination k corresponding to the values of the dummies given in this resut line (to be cleaned one day) */
+     k=TKresult[nres];
-     if(i1 != 1 && TKresult[nres]!= k)
+     if(TKresult[nres]==0) k=1; /* To be checked for noresult */
-       continue;
+     /*  for(k=1; k<=i1;k++){ /\* We want to find the combination k corresponding to the values of the dummies given in this resut line (to be cleaned one day) *\/ */
-     if(invalidvarcomb[k]){
+     /* if(i1 != 1 && TKresult[nres]!= k) */
-       printf("\nCombination (%d) projection ignored because no cases \n",k);
+     /*   continue; */
-       continue;
+     /* if(invalidvarcomb[k]){ */
-     }
+     /*   printf("\nCombination (%d) projection ignored because no cases \n",k);  */
+     /*   continue; */
+     /* } */
      fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#");
-     for(j=1;j<=cptcoveff;j++) {
+     for(j=1;j<=cptcovs;j++){
-       /* fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); */
+       /* for(j=1;j<=cptcoveff;j++) { */
-       fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
+     /*   /\* fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); *\/ */
-     }
+     /*   fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
-     for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+     /* } */
-       fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+     /* for (k4=1; k4<= nsq; k4++){ /\* For each selected (single) quantitative value *\/ */
+     /*   fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); */
+     /* } */
+       fprintf(ficresf," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
      }
      fprintf(ficresf," yearproj age");
      for(j=1; j<=nlstate+ndeath;j++){
        for(i=1; i<=nlstate;i++)
- Line 9968  void prevforecast(char fileres[], double
+ Line 11569  void prevforecast(char fileres[], double
            }
          }
          fprintf(ficresf,"\n");
-         for(j=1;j<=cptcoveff;j++)
+         /* for(j=1;j<=cptcoveff;j++)  */
+         for(j=1;j<=cptcovs;j++)
+           fprintf(ficresf,"%d %lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
            /* fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); /\* Tvaraff not correct *\/ */
-           fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); /* TnsdVar[Tvaraff]  correct */
+           /* fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); /\* TnsdVar[Tvaraff]  correct *\/ */
          fprintf(ficresf,"%.f %.f ",anprojd+yearp,agec+h*hstepm/YEARM*stepm);
          for(j=1; j<=nlstate+ndeath;j++) {
- Line 10007  void prevforecast(char fileres[], double
+ Line 11610  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 10062  void prevforecast(char fileres[], double
+ Line 11665  void prevforecast(char fileres[], double
    /* if(jintmean==0) jintmean=1; */
    /* if(mintmean==0) jintmean=1; */
-   i1=pow(2,cptcoveff);
+   /* i1=pow(2,cptcoveff); */
-   if (cptcovn < 1){i1=1;}
+   /* if (cptcovn < 1){i1=1;} */
    fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jintmean,mintmean,aintmean,dateintmean,dateprev1,dateprev2);
    printf("# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jintmean,mintmean,aintmean,dateintmean,dateprev1,dateprev2);
    fprintf(ficresfb,"#****** Routine prevbackforecast **\n");
-   for(nres=1; nres <= nresult; nres++) /* For each resultline */
+   for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-   for(k=1; k<=i1;k++){
+     k=TKresult[nres];
-     if(i1 != 1 && TKresult[nres]!= k)
+     if(TKresult[nres]==0) k=1; /* To be checked for noresult */
-       continue;
+   /* for(k=1; k<=i1;k++){ */
-     if(invalidvarcomb[k]){
+   /*   if(i1 != 1 && TKresult[nres]!= k) */
-       printf("\nCombination (%d) projection ignored because no cases \n",k);
+   /*     continue; */
-       continue;
+   /*   if(invalidvarcomb[k]){ */
-     }
+   /*     printf("\nCombination (%d) projection ignored because no cases \n",k);  */
+   /*     continue; */
+   /*   } */
      fprintf(ficresfb,"\n#****** hbijx=probability over h years, hb.jx is weighted by observed prev \n#");
-     for(j=1;j<=cptcoveff;j++) {
+     for(j=1;j<=cptcovs;j++){
-       fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
+     /* for(j=1;j<=cptcoveff;j++) { */
-     }
+     /*   fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
-     for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
+     /* } */
-       fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
+       fprintf(ficresfb," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
      }
+    /*  fprintf(ficrespij,"******\n"); */
+    /* for (k4=1; k4<= nsq; k4++){ /\* For each selected (single) quantitative value *\/ */
+    /*    fprintf(ficresfb," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); */
+    /*  } */
      fprintf(ficresfb," yearbproj age");
      for(j=1; j<=nlstate+ndeath;j++){
        for(i=1; i<=nlstate;i++)
- Line 10115  void prevforecast(char fileres[], double
+ Line 11724  void prevforecast(char fileres[], double
            }
          }
          fprintf(ficresfb,"\n");
-         for(j=1;j<=cptcoveff;j++)
+         /* for(j=1;j<=cptcoveff;j++) */
-           fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
+         for(j=1;j<=cptcovs;j++)
+           fprintf(ficresfb,"%d %lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
+           /* fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
          fprintf(ficresfb,"%.f %.f ",anbackd+yearp,agec-h*hstepm/YEARM*stepm);
          for(i=1; i<=nlstate+ndeath;i++) {
            ppij=0.;ppi=0.;
- Line 10685  void printinggnuplotmort(char fileresu[]
+ Line 12296  void printinggnuplotmort(char fileresu[]
    char dirfileres[132],optfileres[132];
-   int ng;
+   /*int ng;*/
    /*#ifdef windows */
- Line 10709  int readdata(char datafile[], int firsto
+ Line 12320  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 11080  int decoderesult( char resultline[], int
+ Line 12691  int decoderesult( char resultline[], int
    if (strlen(resultsav) >1){
      j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' in this resultline */
    }
-   if(j == 0){ /* Resultline but no = */
+   if(j == 0 && cptcovs== 0){ /* Resultline but no =  and no covariate in the model */
      TKresult[nres]=0; /* Combination for the nresult and the model */
      return (0);
    }
    if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */
-     printf("ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, %s.\n",j, cptcovs, model);
+     fprintf(ficlog,"ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, 1+age+%s.\n",j, cptcovs, model);fflush(ficlog);
-     fprintf(ficlog,"ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, %s.\n",j, cptcovs, model);
+     printf("ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, 1+age+%s.\n",j, cptcovs, model);fflush(stdout);
-     /* return 1;*/
+     if(j==0)
+       return 1;
    }
    for(k=1; k<=j;k++){ /* Loop on any covariate of the RESULT LINE */
      if(nbocc(resultsav,'=') >1){
- Line 11325  int decodemodel( char model[], int lasto
+ Line 12937  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 11375  int decodemodel( char model[], int lasto
+ Line 12987  int decodemodel( char model[], int lasto
      if (strlen(modelsav) >1){ /* 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 +V7*V2 +age*V6*V3 +age*V7*V3 +age*V6*V4 +age*V7*V4 */
        j=nbocc(modelsav,'+'); /**< j=Number of '+' */
        j1=nbocc(modelsav,'*'); /**< j1=Number of '*' */
-       cptcovs=j+1-j1; /**<  Number of simple covariates V1 +V1*age +V3 +V3*V4 +age*age => V1 + V3 =4+1-3=2  */
+       cptcovs=0; /**<  Number of simple covariates V1 +V1*age +V3 +V3*V4 +age*age => V1 + V3 =4+1-3=2  Wrong */
        cptcovt= j+1; /* Number of total covariates in the model, not including
                       * cst, age and age*age
                       * V1+V1*age+ V3 + V3*V4+age*age=> 3+1=4*/
- Line 11440  int decodemodel( char model[], int lasto
+ Line 13052  int decodemodel( char model[], int lasto
          Tvar[k]=0; Tprod[k]=0; Tposprod[k]=0;
        }
        cptcovage=0;
+       /* First loop in order to calculate */
+       /* for age*VN*Vm
+        * Provides, Typevar[k], Tage[cptcovage], existcomb[n][m], FixedV[ncovcolt+k12]
+        * Tprod[k1]=k  Tposprod[k]=k1;    Tvard[k1][1] =m;
+       */
+       /* Needs  FixedV[Tvardk[k][1]] */
+       /* For others:
+        * Sets   Typevar[k];
+        * Tvar[k]=ncovcol+nqv+ntv+nqtv+k11;
+        *        Tposprod[k]=k11;
+        *        Tprod[k11]=k;
+        *        Tvardk[k][1] =m;
+        * Needs FixedV[Tvardk[k][1]] == 0
+       */
        for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model line */
          cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' cutl from left to right
                                           modelsav==V2+V1+V5*age+V4+V3*age strb=V3*age stra=V2+V1V5*age+V4 */    /* <model> "V5+V4+V3+V4*V3+V5*age+V1*age+V1" strb="V5" stra="V4+V3+V4*V3+V5*age+V1*age+V1" */
- Line 11465  int decodemodel( char model[], int lasto
+ Line 13093  int decodemodel( char model[], int lasto
                  strcat(strb,stre);
                  strcpy(strd,strb); /* in order for strd to not be "age"  for next test (will be Vn*Vm */
                }
-               printf("DEBUG FIXED k=%d, Tage[k]=%d, Tvar[Tage[k]=%d,FixedV[Tvar[Tage[k]]]=%d\n",k,Tage[k],Tvar[Tage[k]],FixedV[Tvar[Tage[k]]]);
+               /* printf("DEBUG FIXED k=%d, Tage[k]=%d, Tvar[Tage[k]=%d,FixedV[Tvar[Tage[k]]]=%d\n",k,Tage[k],Tvar[Tage[k]],FixedV[Tvar[Tage[k]]]); */
-               FixedV[Tvar[Tage[k]]]=0; /* HERY not sure */
+               /* FixedV[Tvar[Tage[k]]]=0; /\* HERY not sure if V7*V4*age Fixed might not exist  yet*\/ */
              }else{  /* strc=Vn*Vm (and strd=age) and should be strb=Vn*Vm but want to keep original strb double product  */
                strcpy(stre,strb); /* save full b in stre */
                strcpy(strb,strc); /* save short c in new short b for next block strb=Vn*Vm*/
- Line 11499  int decodemodel( char model[], int lasto
+ Line 13127  int decodemodel( char model[], int lasto
                Tvardk[k][1] =m; /* m 1 for V1*/
                Tvard[k1][2] =n; /* n 4 for V4*/
                Tvardk[k][2] =n; /* n 4 for V4*/
- /*            Tvar[Tage[cptcovage]]=k1;*/ /* Tvar[6=age*V3*V2]=9 (new fixed covariate) */
+ /*            Tvar[Tage[cptcovage]]=k1;*/ /* Tvar[6=age*V3*V2]=9 (new fixed covariate) */ /* We don't know about Fixed yet HERE */
                if( FixedV[Tvardk[k][1]] == 0 && FixedV[Tvardk[k][2]] == 0){ /* If the product is a fixed covariate then we feed the new column with Vn*Vm */
                  for (i=1; i<=lastobs;i++){/* For fixed product */
                    /* Computes the new covariate which is a product of
- Line 11649  int decodemodel( char model[], int lasto
+ Line 13277  int decodemodel( char model[], int lasto
                                  /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
                                    scanf("%d",i);*/
        } /* end of loop + on total covariates */
      } /* end if strlen(modelsave == 0) age*age might exist */
    } /* end if strlen(model == 0) */
    cptcovs=cptcovt - cptcovdageprod - cptcovprod;/**<  Number of simple covariates V1 +V1*age +V3 +V3*V4 +age*age + age*v4*V3=> V1 + V3 =4+1-3=2  */
- Line 11687  Fixed[k] 0=fixed (product or simple), 1
+ Line 13317  Fixed[k] 0=fixed (product or simple), 1
  Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model);
    for(k=-1;k<=NCOVMAX; k++){ Fixed[k]=0; Dummy[k]=0;}
    for(k=1;k<=NCOVMAX; k++){TvarFind[k]=0; TvarVind[k]=0;}
+   /* Second loop for calculating  Fixed[k], Dummy[k]*/
    for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0,ncovva=0,ncovvta=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0, ncovvt=0;k<=cptcovt; k++){ /* or cptocvt loop on k from model */
      if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */
        Fixed[k]= 0;
- Line 12730  int hPijx(double *p, int bage, int fage)
+ Line 14365  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 12936  int main(int argc, char *argv[])
+ Line 14571  int main(int argc, char *argv[])
    double fret;
    double dum=0.; /* Dummy variable */
-   double ***p3mat;
+   /* double*** p3mat;*/
    /* double ***mobaverage; */
    double wald;
-   char line[MAXLINE];
+   char line[MAXLINE], linetmp[MAXLINE];
    char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE];
    char  modeltemp[MAXLINE];
- Line 12949  int main(int argc, char *argv[])
+ Line 14584  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 12984  int main(int argc, char *argv[])
+ Line 14619  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 13042  int main(int argc, char *argv[])
+ Line 14677  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 13273  int main(int argc, char *argv[])
+ Line 14908  int main(int argc, char *argv[])
      }else
        break;
    }
-   if((num_filled=sscanf(line,"model=1+age%[^.\n]", model)) !=EOF){
+   if((num_filled=sscanf(line,"model=%[^.\n]", model)) !=EOF){ /* Every character after model but dot and  return */
+     if (num_filled != 1){
+       printf("ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line);
+       fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line);
+       model[0]='\0';
+       goto end;
+     }else{
+       trimbtab(linetmp,line); /* Trims multiple blanks in line */
+       strcpy(line, linetmp);
+     }
+   }
+   if((num_filled=sscanf(line,"model=1+age%[^.\n]", model)) !=EOF){ /* Every character after 1+age but dot and  return */
      if (num_filled != 1){
        printf("ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line);
        fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line);
- Line 13646  Please run with mle=-1 to get a correct
+ Line 15292  Please run with mle=-1 to get a correct
    Tvard=imatrix(1,NCOVMAX,1,2); /* n=Tvard[k1][1]  and m=Tvard[k1][2] gives the couple n,m of the k1 th product Vn*Vm
                              * For V3*V2 (in V2+V1+V1*V4+age*V3+V3*V2), V3*V2 position is 2nd.
                              * Tvard[k1=2][1]=3 (V3) Tvard[k1=2][2]=2(V2) */
-   Tvardk=imatrix(-1,NCOVMAX,1,2);
+   Tvardk=imatrix(0,NCOVMAX,1,2);
    Tage=ivector(1,NCOVMAX); /* Gives the covariate id of covariates associated with age: V2 + V1 + age*V4 + V3*age
 covariates (3 plus signs)
                           Tage[1=V3*age]= 4; Tage[2=age*V4] = 3
- Line 13933  This file: <a href=\"%s\">%s</a></br>Tit
+ Line 15579  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 14024  Interval (in months) between two waves:
+ Line 15670  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 14127  Interval (in months) between two waves:
+ Line 15773  Interval (in months) between two waves:
      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 14581  Please run with mle=-1 to get a correct
+ Line 16227  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 14901  Please run with mle=-1 to get a correct
+ Line 16547  Please run with mle=-1 to get a correct
      pstamp(ficreseij);
-     i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */
+     /* i1=pow(2,cptcoveff); /\* Number of combination of dummy covariates *\/ */
-     if (cptcovn < 1){i1=1;}
+     /* if (cptcovn < 1){i1=1;} */
-     for(nres=1; nres <= nresult; nres++) /* For each resultline */
+     for(nres=1; nres <= nresult; nres++){ /* For each resultline */
-     for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
+     /* for(k=1; k<=i1;k++){ /\* For any combination of dummy covariates, fixed and varying *\/ */
-       if(i1 != 1 && TKresult[nres]!= k)
+       /* if(i1 != 1 && TKresult[nres]!= k) */
-         continue;
+       /*        continue; */
        fprintf(ficreseij,"\n#****** ");
        printf("\n#****** ");
-       for(j=1;j<=cptcoveff;j++) {
+       for(j=1;j<=cptcovs;j++){
-         fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
+       /* for(j=1;j<=cptcoveff;j++) { */
-         printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);
+         /* fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
+         fprintf(ficreseij," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
+         printf(" V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]);
+         /* printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */
        }
        for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
          printf(" V%d=%lg ",TvarsQ[j], TinvDoQresult[nres][TvarsQ[j]]); /* TvarsQ[j] gives the name of the jth quantitative (fixed or time v) */
- Line 14981  Please run with mle=-1 to get a correct
+ Line 16630  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 15156  Please run with mle=-1 to get a correct
+ Line 16805  Please run with mle=-1 to get a correct
      free_vector(weight,firstobs,lastobs);
-     free_imatrix(Tvardk,-1,NCOVMAX,1,2);
+     free_imatrix(Tvardk,0,NCOVMAX,1,2);
      free_imatrix(Tvard,1,NCOVMAX,1,2);
      free_imatrix(s,1,maxwav+1,firstobs,lastobs);
      free_matrix(anint,1,maxwav,firstobs,lastobs);
- Line 15178  Please run with mle=-1 to get a correct
+ Line 16827  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 15239  Please run with mle=-1 to get a correct
+ Line 16889  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>

Removed from v.1.350
changed lines
	Added in v.1.359