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 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $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    Revision 1.350  2023/04/24 11:38:06  brouard
   *** empty log message ***    *** 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 POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */
 /* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */  /* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */
 /* #define FLATSUP  *//* Suppresses directions where likelihood is flat */  /* #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 <math.h>
 #include <stdio.h>  #include <stdio.h>
Line 1365  double gnuplotversion=GNUPLOTVERSION; Line 1394  double gnuplotversion=GNUPLOTVERSION;
 /* $State$ */  /* $State$ */
 #include "version.h"  #include "version.h"
 char version[]=__IMACH_VERSION__;  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 fullversion[]="$Revision$ $Date$"; 
 char strstart[80];  char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];  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 maxwav=0; /* Maxim number of waves */
 int jmin=0, jmax=0; /* min, max spacing between 2 waves */  int jmin=0, jmax=0; /* min, max spacing between 2 waves */
 int ijmin=0, ijmax=0; /* Individuals having jmin and jmax */   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)*/                     to the likelihood and the sum of weights (done by funcone)*/
 int mle=1, weightopt=0;  int mle=1, weightopt=0;
 int **mw; /* mw[mi][i] is number of the mi wave for this individual */  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 */  /* 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"*/  /* 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[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]@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} */  /* 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 */  /* 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;    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) */  /* char *substrchaine(char *out, char *in, char *chain) */
 /* { */  /* { */
 /*   /\* Substract chain 'chain' from 'in', return and output 'out' *\/ */  /*   /\* 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);     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 ************************/  /*************** powell ************************/
 /*  /*
Line 2637  void powell(double p[], double **xi, int Line 4187  void powell(double p[], double **xi, int
     curr_time = *localtime(&rcurr_time);      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); */      /* 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); */      /* 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);      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(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);      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);          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  #ifdef DEBUG
       printf("fret=%lf, %lf, %lf \n", *fret, *fret, *fret);        printf("fret=%lf, %lf, %lf \n", *fret, *fret, *fret);
       fprintf(ficlog, "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 */        printf("%d",i);fflush(stdout); /* print direction (parameter) i */
       fprintf(ficlog,"%d",i);fflush(ficlog);        fprintf(ficlog,"%d",i);fflush(ficlog);
 #ifdef LINMINORIGINAL  #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  #else
       linmin(p,xit,n,fret,func,&flat); /* Point p[n]. xit[n] has been loaded for direction i as input.*/        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  #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 */        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  #ifdef DEBUG
       printf("%d %.12e",i,(*fret));        printf("%d %.12e",i,(*fret));
       fprintf(ficlog,"%d %.12e",i,(*fret));        fprintf(ficlog,"%d %.12e",i,(*fret));
       for (j=1;j<=n;j++) {        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++) {        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");        printf("\n");
       fprintf(ficlog,"\n");        fprintf(ficlog,"\n");
 #endif  #endif
     } /* end loop on each direction i */      } /* 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  */      /* 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) */      /* New value of last point Pn is not computed, P(n-1) */
     for(j=1;j<=n;j++) {      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;         return; 
     } /* enough precision */       } /* enough precision */ 
     if (*iter == ITMAX*n) nrerror("powell exceeding maximum iterations.");       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];         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 */      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) {                  if (*iter <=4) {
 #else  #else
 #endif  #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); */        /* 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 */        /*  Even if f3 <f1, directest can be negative and t >0 */
       /* mu² and del² are equal when f3=f1 */        /* 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  #ifdef NRCORIGINAL
       t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)- del*SQR(fp-fptt); /* Original Numerical Recipes in C*/        t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)- del*SQR(fp-fptt); /* Original Numerical Recipes in C*/
 #else  #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 */        if (t < 0.0) { /* Then we use it for new direction */
 #else  #else
       if (directest*t < 0.0) { /* Contradiction between both tests */        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);          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,"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);          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");          fprintf(ficlog,"\n");
 #endif  #endif
       } /* end of t or directest negative */        } /* 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  #ifdef POWELLNOF3INFF1TEST
 #else  #else
       } /* end if (fptt < fp)  */        } /* end if (fptt < fp)  */
Line 3140  void powell(double p[], double **xi, int Line 4702  void powell(double p[], double **xi, int
     first++;      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(min,1,nlstate);
   free_vector(max,1,nlstate);    free_vector(max,1,nlstate);
   free_vector(meandiff,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} */    /*  0.51326036147820708, 0.48673963852179264} */
   /* If we start from prlim again, prlim tends to a constant matrix */    /* 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;    int first=0;
   double *min, *max, *meandiff, maxmax,sumnew=0.;    double *min, *max, *meandiff, maxmax,sumnew=0.;
   /* double **matprod2(); */ /* test */    /* 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.    /* 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.     * 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 sumnew=0.;
   double agefin;    double agefin;
   double k3=0.; /* constant of the w_x diagonal matrix (in order for B to sum to 1 even for death state) */    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 **out, cov[NCOVMAX+1];
   double **newm;    double **newm;
   double agexact;    double agexact;
   double agebegin, ageend;    /*double agebegin, ageend;*/
   
   /* Hstepm could be zero and should return the unit matrix */    /* Hstepm could be zero and should return the unit matrix */
   for (i=1;i<=nlstate+ndeath;i++)    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       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 **out, cov[NCOVMAX+1], **bmij();
   double **newm, ***newmm;    double **newm, ***newmm;
   double agexact;    double agexact;
   double agebegin, ageend;    /*double agebegin, ageend;*/
   double **oldm, **savm;    double **oldm, **savm;
   
   newmm=po; /* To be saved */    newmm=po; /* To be saved */
Line 4361  double func( double *x) Line 5925  double func( double *x)
 double funcone( double *x)  double funcone( double *x)
 {  {
   /* Same as func but slower because of a lot of printf and if */    /* 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 ioffset=0;
   int ipos=0,iposold=0,ncovv=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 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 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          * 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[]          * 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 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]          * 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}          *                     6, 8, 9, 10, 11}
         * TvarFind[itv]                        0      0       0          * 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       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[ncovf]]=[1]=2 [2]=3 [4]=4
         * Tvar[TvarFind[itv]]                [0]=?      ?ncovv 1 à ncovvt]          * 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}          *   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*/          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(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 */          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) */             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 */          }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 *\/ */            /* 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 */            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 */          if(ipos!=iposold){ /* Not a product or first of a product */
           cotvarvold=cotvarv;            cotvarvold=cotvarv;
Line 4568  double funcone( double *x) Line 6137  double funcone( double *x)
         }          }
         iposold=ipos;          iposold=ipos;
         cov[ioffset+ipos]=cotvarv;          cov[ioffset+ipos]=cotvarv;
           /* printf("DEBUG Product cov[ioffset+ipos=%d] \n",ioffset+ipos); */
         /* For products */          /* For products */
       }        }
       /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates single *\/ */        /* 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 []))  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 **xi;
   double fret;    /*double fret;*/
   double fretone; /* Only one call to likelihood */    /*double fretone;*/ /* Only one call to likelihood */
   /*  char filerespow[FILENAMELENGTH];*/    /*  char filerespow[FILENAMELENGTH];*/
     
     /*double * p1;*/ /* Shifted parameters from 0 instead of 1 */
 #ifdef NLOPT  #ifdef NLOPT
   int creturn;    int creturn;
   nlopt_opt opt;    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[9] = { -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL }; /\* lower bounds *\/ */
   double *lb;    double *lb;
   double minf; /* the minimum objective value, upon return */    double minf; /* the minimum objective value, upon return */
   double * p1; /* Shifted parameters from 0 instead of 1 */  
   myfunc_data dinst, *d = &dinst;    myfunc_data dinst, *d = &dinst;
 #endif  #endif
   
   
   xi=matrix(1,npar,1,npar);    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++)      for (j=1;j<=npar;j++)
       xi[i][j]=(i==j ? 1.0 : 0.0);        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_");     strcpy(filerespow,"POW_"); 
   strcat(filerespow,fileres);    strcat(filerespow,fileres);
   if((ficrespow=fopen(filerespow,"w"))==NULL) {    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);    powell(p,xi,npar,ftol,&iter,&fret,flatdir,func);
 #else  /* FLATSUP */  #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 */  #endif  /* FLATSUP */
   
 #ifdef LINMINORIGINAL  #ifdef LINMINORIGINAL
Line 5252  double hessij( double x[], double **hess Line 6839  double hessij( double x[], double **hess
       kmax=kmax+10;        kmax=kmax+10;
     }      }
     if(kmax >=10 || firstime ==1){      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);        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);        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);        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 i, m, jk, j1, bool, z1,j, iv;
   int mi; /* Effective wave */    int mi; /* Effective wave */
   int iage;    int iage;
   double agebegin, ageend;    double agebegin; /*, ageend;*/
   
   double **prop;    double **prop;
   double posprop;     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 */              if(j==0) j=1;  /* Survives at least one month after exam */
             else if(j<0){              else if(j<0){
               nberr++;                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 */                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);                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);                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;              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]);*/            /*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){            if(j<0){
             nberr++;              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;            sum=sum+j;
         }          }
Line 7589  void varprob(char optionfilefiname[], do Line 9177  void varprob(char optionfilefiname[], do
    double ***varpij;     double ***varpij;
   
    strcpy(fileresprob,"PROB_");      strcpy(fileresprob,"PROB_"); 
    strcat(fileresprob,fileres);     strcat(fileresprob,fileresu);
    if((ficresprob=fopen(fileresprob,"w"))==NULL) {     if((ficresprob=fopen(fileresprob,"w"))==NULL) {
      printf("Problem with resultfile: %s\n", fileresprob);       printf("Problem with resultfile: %s\n", fileresprob);
      fprintf(ficlog,"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 , \                    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 jprev1, double mprev1,double anprev1, double dateprev1, double dateprojd, double dateback1, \
                   double jprev2, double mprev2,double anprev2, double dateprev2, double dateprojf, double dateback2){                    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 */    /* 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 \     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 \     <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);   <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 */       /* Survival functions (period) in state j */
      for(cpt=1; cpt<=nlstate;cpt++){       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," (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);         fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
      }       }
      /* State specific survival functions (period) */       /* State specific survival functions (period) */
      for(cpt=1; cpt<=nlstate;cpt++){       for(cpt=1; cpt<=nlstate;cpt++){
        fprintf(fichtm,"<br>\n- Survival functions in state %d and in any other live state (total).\         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);   <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," (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);         fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
      }       }
      /* Period (forward stable) prevalence in each health state */       /* Period (forward stable) prevalence in each health state */
      for(cpt=1; cpt<=nlstate;cpt++){       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," (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);        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 */        /* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */
        for(cpt=1; cpt<=nlstate;cpt++){         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), \           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);  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," (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);           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 **************/  /******************* 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){  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 cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,k4=0,kf=0,kvar=0,kk=0,ipos=0,iposold=0,ij=0, ijp=0, l=0;
   int lv=0, vlv=0, kl=0;    int lv=0, vlv=0, kl=0;
   int ng=0;    int ng=0;
Line 8409  void printinggnuplot(char fileresu[], ch Line 9997  void printinggnuplot(char fileresu[], ch
     kvar=Tvar[TvarFind[kf]]; /* variable name */      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+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=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 ++) {      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,"\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));        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);               fprintf(ficgp," u %d:(",ioffset); 
             kl=0;              kl=0;
             strcpy(gplotcondition,"(");              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= 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,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(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 */                /* 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]][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++;                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++;                kl++;
               if(k <cptcoveff && cptcoveff>1)                if(k <cptcovs && cptcovs>1)
                 sprintf(gplotcondition+strlen(gplotcondition)," && ");                  sprintf(gplotcondition+strlen(gplotcondition)," && ");
             }              }
             strcpy(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{            }else{
             fprintf(ficgp,",\\\n '' ");              fprintf(ficgp,",\\\n '' ");
           }            }
           if(cptcoveff ==0){ /* No covariate */            /* if(cptcoveff ==0){ /\* No covariate *\/ */
             if(cptcovs ==0){ /* No covariate */
             ioffset=2; /* Age is in 2 */              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*/              /*# 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 */              /*#   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,"#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,"#model=1+age+%s \n",model);
     fprintf(ficgp,"# Type of graphic ng=%d\n",ng);      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(k1=1; k1 <=m; k1++)  /\* For each combination of covariate *\/ */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */      for(nres=1; nres <= nresult; nres++){ /* For each resultline */
      /* k1=nres; */       /* k1=nres; */
Line 9864  void prevforecast(char fileres[], double Line 11459  void prevforecast(char fileres[], double
   */    */
   /* double anprojd, mprojd, jprojd; */    /* double anprojd, mprojd, jprojd; */
   /* double anprojf, mprojf, jprojf; */    /* 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 agec; /* generic age */
   double agelim, ppij, yp,yp1,yp2;    double agelim, ppij;
   double *popeffectif,*popcount;    /*double *popcount;*/
   double ***p3mat;    double ***p3mat;
   /* double ***mobaverage; */    /* double ***mobaverage; */
   char fileresf[FILENAMELENGTH];    char fileresf[FILENAMELENGTH];
Line 9920  void prevforecast(char fileres[], double Line 11515  void prevforecast(char fileres[], double
   /* date2dmy(dateintmean,&jintmean,&mintmean,&aintmean); */    /* date2dmy(dateintmean,&jintmean,&mintmean,&aintmean); */
   /* date2dmy(dateprojd,&jprojd, &mprojd, &anprojd); */    /* date2dmy(dateprojd,&jprojd, &mprojd, &anprojd); */
   /* date2dmy(dateprojf,&jprojf, &mprojf, &anprojf); */    /* 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,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jintmean,mintmean,aintmean,dateintmean,dateprev1,dateprev2); 
       
   fprintf(ficresf,"#****** Routine prevforecast **\n");    fprintf(ficresf,"#****** Routine prevforecast **\n");
       
 /*            if (h==(int)(YEARM*yearp)){ */  /*            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#");      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");      fprintf(ficresf," yearproj age");
     for(j=1; j<=nlstate+ndeath;j++){       for(j=1; j<=nlstate+ndeath;j++){ 
       for(i=1; i<=nlstate;i++)                for(i=1; i<=nlstate;i++)        
Line 9968  void prevforecast(char fileres[], double Line 11569  void prevforecast(char fileres[], double
           }            }
         }          }
         fprintf(ficresf,"\n");          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,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);          fprintf(ficresf,"%.f %.f ",anprojd+yearp,agec+h*hstepm/YEARM*stepm);
                   
         for(j=1; j<=nlstate+ndeath;j++) {          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).       anback2 year of end of backprojection (same day and month as back1).
      prevacurrent and prev are prevalences.       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 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 ***p3mat;
   /* double ***mobaverage; */    /* double ***mobaverage; */
   char fileresfb[FILENAMELENGTH];    char fileresfb[FILENAMELENGTH];
Line 10062  void prevforecast(char fileres[], double Line 11665  void prevforecast(char fileres[], double
   /* if(jintmean==0) jintmean=1; */    /* if(jintmean==0) jintmean=1; */
   /* if(mintmean==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);    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);    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");    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#");      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");      fprintf(ficresfb," yearbproj age");
     for(j=1; j<=nlstate+ndeath;j++){      for(j=1; j<=nlstate+ndeath;j++){
       for(i=1; i<=nlstate;i++)        for(i=1; i<=nlstate;i++)
Line 10115  void prevforecast(char fileres[], double Line 11724  void prevforecast(char fileres[], double
           }            }
         }          }
         fprintf(ficresfb,"\n");          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);          fprintf(ficresfb,"%.f %.f ",anbackd+yearp,agec-h*hstepm/YEARM*stepm);
         for(i=1; i<=nlstate+ndeath;i++) {          for(i=1; i<=nlstate+ndeath;i++) {
           ppij=0.;ppi=0.;            ppij=0.;ppi=0.;
Line 10685  void printinggnuplotmort(char fileresu[] Line 12296  void printinggnuplotmort(char fileresu[]
   
   char dirfileres[132],optfileres[132];    char dirfileres[132],optfileres[132];
   
   int ng;    /*int ng;*/
   
   
   /*#ifdef windows */    /*#ifdef windows */
Line 10709  int readdata(char datafile[], int firsto Line 12320  int readdata(char datafile[], int firsto
   /*-------- data file ----------*/    /*-------- data file ----------*/
   FILE *fic;    FILE *fic;
   char dummy[]="                         ";    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 lstra;
   int linei, month, year,iout;    int linei, month, year,iout;
   int noffset=0; /* This is the offset if BOM data file */    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){    if (strlen(resultsav) >1){
     j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' in this resultline */      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 */      TKresult[nres]=0; /* Combination for the nresult and the model */
     return (0);      return (0);
   }    }
   if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */    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 */    for(k=1; k<=j;k++){ /* Loop on any covariate of the RESULT LINE */
     if(nbocc(resultsav,'=') >1){      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 */  /* 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 n,m;
   int  j1, k1, k11, k12, k2, k3, k4;    int  j1, k1, k11, k12, k2, k3, k4;
   char modelsav[300];    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 */      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 '+' */        j=nbocc(modelsav,'+'); /**< j=Number of '+' */
       j1=nbocc(modelsav,'*'); /**< j1=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        cptcovt= j+1; /* Number of total covariates in the model, not including
                      * cst, age and age*age                        * cst, age and age*age 
                      * V1+V1*age+ V3 + V3*V4+age*age=> 3+1=4*/                       * 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;          Tvar[k]=0; Tprod[k]=0; Tposprod[k]=0;
       }        }
       cptcovage=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 */        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          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" */                                           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);                  strcat(strb,stre);
                 strcpy(strd,strb); /* in order for strd to not be "age"  for next test (will be Vn*Vm */                  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  */              }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(stre,strb); /* save full b in stre */
               strcpy(strb,strc); /* save short c in new short b for next block strb=Vn*Vm*/                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*/                Tvardk[k][1] =m; /* m 1 for V1*/
               Tvard[k1][2] =n; /* n 4 for V4*/                Tvard[k1][2] =n; /* n 4 for V4*/
               Tvardk[k][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 */                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 */                  for (i=1; i<=lastobs;i++){/* For fixed product */
                   /* Computes the new covariate which is a product of                    /* 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);                                  /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
                                   scanf("%d",i);*/                                    scanf("%d",i);*/
       } /* end of loop + on total covariates */        } /* end of loop + on total covariates */
   
         
     } /* end if strlen(modelsave == 0) age*age might exist */      } /* end if strlen(modelsave == 0) age*age might exist */
   } /* end if strlen(model == 0) */    } /* 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  */    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);  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++){ Fixed[k]=0; Dummy[k]=0;}
   for(k=1;k<=NCOVMAX; k++){TvarFind[k]=0; TvarVind[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 */    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 */      if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */
       Fixed[k]= 0;        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 agelim;
   int hstepm;    int hstepm;
   int nhstepm;    int nhstepm;
   int h, i, i1, j, k, k4, nres=0;    int h, i, i1, j, k, nres=0;
   
   double agedeb;    double agedeb;
   double ***p3mat;    double ***p3mat;
Line 12936  int main(int argc, char *argv[]) Line 14571  int main(int argc, char *argv[])
   
   double fret;    double fret;
   double dum=0.; /* Dummy variable */    double dum=0.; /* Dummy variable */
   double ***p3mat;    /* double*** p3mat;*/
   /* double ***mobaverage; */    /* double ***mobaverage; */
   double wald;    double wald;
   
   char line[MAXLINE];    char line[MAXLINE], linetmp[MAXLINE];
   char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE];    char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE];
   
   char  modeltemp[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 pathr[MAXLINE], pathimach[MAXLINE]; 
   char *tok, *val; /* pathtot */    char *tok, *val; /* pathtot */
   /* int firstobs=1, lastobs=10; /\* nobs = lastobs-firstobs declared globally ;*\/ */    /* int firstobs=1, lastobs=10; /\* nobs = lastobs-firstobs declared globally ;*\/ */
   int c,  h , cpt, c2;    int c, h; /* c2; */
   int jl=0;    int jl=0;
   int i1, j1, jk, stepsize=0;    int i1, j1, jk, stepsize=0;
   int count=0;    int count=0;
Line 12984  int main(int argc, char *argv[]) Line 14619  int main(int argc, char *argv[])
   double ***delti3; /* Scale */    double ***delti3; /* Scale */
   double *delti; /* Scale */    double *delti; /* Scale */
   double ***eij, ***vareij;    double ***eij, ***vareij;
   double **varpl; /* Variances of prevalence limits by age */    //double **varpl; /* Variances of prevalence limits by age */
   
   double *epj, vepp;    double *epj, vepp;
   
Line 13042  int main(int argc, char *argv[]) Line 14677  int main(int argc, char *argv[])
   getcwd(pathcd, size);    getcwd(pathcd, size);
 #endif  #endif
   syscompilerinfo(0);    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){    if(argc <=1){
     printf("\nEnter the parameter file name: ");      printf("\nEnter the parameter file name: ");
     if(!fgets(pathr,FILENAMELENGTH,stdin)){      if(!fgets(pathr,FILENAMELENGTH,stdin)){
Line 13273  int main(int argc, char *argv[]) Line 14908  int main(int argc, char *argv[])
     }else      }else
       break;        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){      if (num_filled != 1){
       printf("ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line);        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);        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    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.                               * 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) */                              * 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    Tage=ivector(1,NCOVMAX); /* Gives the covariate id of covariates associated with age: V2 + V1 + age*V4 + V3*age
                          4 covariates (3 plus signs)                           4 covariates (3 plus signs)
                          Tage[1=V3*age]= 4; Tage[2=age*V4] = 3                           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     /* Calculates basic frequencies. Computes observed prevalence at single age 
                  and for any valid combination of covariates                   and for any valid combination of covariates
      and prints on file fileres'p'. */       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);                firstpass, lastpass,  stepm,  weightopt, model);
   
   fprintf(fichtm,"\n");    fprintf(fichtm,"\n");
Line 14024  Interval (in months) between two waves: Line 15670  Interval (in months) between two waves:
 #ifdef GSL  #ifdef GSL
     printf("GSL optimization\n");  fprintf(ficlog,"Powell\n");      printf("GSL optimization\n");  fprintf(ficlog,"Powell\n");
 #else  #else
     printf("Powell\n");  fprintf(ficlog,"Powell\n");      printf("Powell-mort\n");  fprintf(ficlog,"Powell-mort\n");
 #endif  #endif
     strcpy(filerespow,"POW-MORT_");       strcpy(filerespow,"POW-MORT_"); 
     strcat(filerespow,fileresu);      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(i=1; i <=NDIM; i++)
       for(j=i+1;j<=NDIM;j++)        for(j=i+1;j<=NDIM;j++)
                                 matcov[i][j]=matcov[j][i];          matcov[i][j]=matcov[j][i];
           
     printf("\nCovariance matrix\n ");      printf("\nCovariance matrix\n ");
     fprintf(ficlog,"\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 */      /* 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 */        /* 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);      precov=matrix(1,MAXRESULTLINESPONE,1,NCOVMAX+1);
     endishere=0;      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);      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#****** ");        fprintf(ficreseij,"\n#****** ");
       printf("\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 */        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) */          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 */        if(i1 != 1 && TKresult[nres]!= k) /* TKresult[nres] is the combination of this nres resultline. All the i1 combinations are not output */
         continue;          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 */        /* 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<=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 */        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_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(Tvard,1,NCOVMAX,1,2);
     free_imatrix(s,1,maxwav+1,firstobs,lastobs);      free_imatrix(s,1,maxwav+1,firstobs,lastobs);
     free_matrix(anint,1,maxwav,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);      free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);
   }  /* mle==-3 arrives here for freeing */    }  /* mle==-3 arrives here for freeing */
   /* endfree:*/    /* 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(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
   free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);    free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
   free_matrix(savms, 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(TmodelInvind,1,NCOVMAX);
   free_ivector(TmodelInvQind,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(nbcode,0,NCOVMAX,0,NCOVMAX);
   /* free_imatrix(codtab,1,100,1,10); */    /* free_imatrix(codtab,1,100,1,10); */

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


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