[BACK]Return to imachprax.c CVS log [TXT] [DIR] Up to [Local Repository] / imach / src

Diff for /imach/src/imachprax.c between versions 1.5 and 1.6

version 1.5, 2023/10/09 09:10:01 version 1.6, 2024/04/24 21:10:29
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     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    Revision 1.5  2023/10/09 09:10:01  brouard
   Summary: trying to reconsider    Summary: trying to reconsider
   
Line 1436  int *wav; /* Number of waves for this in Line 1439  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 2623  void linmin(double p[], double xi[], int Line 2627  void linmin(double p[], double xi[], int
   free_vector(pcom,1,n);     free_vector(pcom,1,n); 
 }   } 
   
 /**** praxis ****/  /**** praxis gegen ****/
 # include <float.h>  
   
 void transpose_in_place ( int n, double **a )  /* 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>
   Purpose:  #include <stdlib.h>
   #include <float.h> /* for DBL_EPSILON */
   /* #include "machine.h" */
   
    TRANSPOSE_IN_PLACE transposes a square matrix in place.  
   Licensing:  
     This code is distributed under the GNU LGPL license.  
   
     Input, int N, the number of rows and columns of the matrix A.  /* 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;
   
     Input/output, double A[N*N], the matrix to be transposed.  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 i;
    
     int i2hi;
     int ihi;
     int ilo;
     int i2lo;
     int jlo=1;
   int j;    int j;
   double t;    int j2hi;
     int jhi;
   /* for ( j = 0; j < n; j++ ){ */    int j2lo;
   /*   for ( i = 0; i < j; i++ ) { */    ilo=1;
   for ( j = 1; j <= n; j++ ){    ihi=n;
     for ( i = 1; i < j; i++ ) {    jlo=1;
       /* t        = a[i+j*n]; */    jhi=n;
       /* a[i+j*n] = a[j+i*n]; */    
       /* a[j+i*n] = t; */    printf ("\n" );
       t        = a[i][j];    printf ("%s\n", s );
       a[i][j] = a[j][i];    for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX )
       a[j][i] = t;    {
       j2hi = j2lo + INCX - 1;
       if ( n < j2hi )
       {
         j2hi = n;
       }
       if ( jhi < j2hi )
       {
         j2hi = jhi;
     }      }
   }  
   return;  
 }  
   
 double pythag( double x, double y )      /* 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" );
 /*  /*
   Purpose:    Determine the range of the rows in this strip.
     R8_HYPOT returns the value of sqrt ( X^2 + Y^2 ).  
   Licensing:  
     This code is distributed under the GNU LGPL license.  
   Modified:  
     26 March 2012  
   Author:  
     John Burkardt  
   Parameters:  
     Input, double X, Y, the arguments.  
     Output, double R8_HYPOT, the value of sqrt ( X^2 + Y^2 ).  
 */  */
 {      if ( 1 < ilo ){
   double a;        i2lo = ilo;
   double b;      }else{
   double value;        i2lo = 1;
       }
   if ( fabs ( x ) < fabs ( y ) )  {      if ( m < ihi ){
     a = fabs ( y );        i2hi = m;
     b = fabs ( x );      }else{
   }  else  {        i2hi = ihi;
     a = fabs ( x );      }
     b = fabs ( y );  
   }      for ( i = i2lo; i <= i2hi; i++ ){
 /*  /*
   A contains the larger value.    Print out (up to) 5 entries in row I, that lie in the current strip.
 */  */
   if ( a == 0.0 )  {        /* fprintf ( ficlog, "%5d:", i - 1 ); */
     value = 0.0;        /* printf ("%5d:", i - 1 ); */
   }  else  {        printf ("%5d:", i );
     value = a * sqrt ( 1.0 + ( b / a ) * ( b / a ) );        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" );
       }
   }    }
   return value;   
      /* 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 svsort ( int n, double d[], double **v )   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 */
 /*  {
   Purpose:   
   
     SVSORT descending sorts D and adjusts the corresponding columns of V.     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.;
   
   Discussion:     /* printf("\n"); */
     A simple bubble sort is used on D.     /* printf("... chi square reduced to ... %20.10e\n", fx); */
     In our application, D contains singular values, and the columns of V are     /* printf("... after %u function calls ...\n", nf); */
     the corresponding right singular vectors.     /* printf("... including %u linear searches ...\n", nl); */
   Author:     /* printf("%10d    %10d%14.7g",nl, nf, fx); */
     Original FORTRAN77 version by Richard Brent.    printf ( "\n" );
     Richard Brent,    printf ( "  Linear searches      %d", nl );
     Algorithms for Minimization with Derivatives,    /* printf ( "  Linear searches      %d\n", nl ); */
     Prentice Hall, 1973,    /* printf ( "  Function evaluations %d\n", nf ); */
     Reprinted by Dover, 2002.    /* printf ( "  Function value FX = %g\n", fx ); */
     printf ( "  Function evaluations %d", nf );
   Parameters:    printf ( "  Function value FX = %.12lf\n", fx );
     Input, int N, the length of D, and the order of V.  #ifdef DEBUGPRAX
     Input/output, double D[N], the vector to be sorted.       printf("n=%d prin=%d\n",n,prin);
     On output, the entries of D are in descending order.  #endif
      if(fx <= fmin) printf(" UNDEFINED "); else  printf("%14.7g",log(fx-fmin));
     Input/output, double V[N,N], an N by N array to be adjusted      if ( n <= 4 || 2 < prin )
     as D is sorted.  In particular, if the value that was in D(I) on input is     {
     moved to D(J) on output, then the input column V(*,I) is moved to       /* for(i=1;i<=n;i++)printf("%14.7g",x[i-1]); */
     the output column V(*,J).       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]; */
 {  {
   int i, j1, j2, j3;     /* init global extern variables and parameters */
   double t;     /* 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();
   
   for (j1 = 1; j1 < n; j1++) {  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++) {
     /*      /*
      * Find J3, the index of the largest entry in D[J1:N-1].        The inner loop starts here.
      * MAXLOC apparently requires its output to be an array.  
     */      */
     j3 = j1;  #ifdef DEBUGPRAX
     for (j2 = j1+1; j2 < n; j2++) {        printf("      The inner loop  here from k=%d to n=%d.\n",k,n);
       if (d[j3] < d[j2]) {        /* fprintf(ficlog,"      The inner loop  here from k=%d to n=%d.\n",k,n); */
         j3 = j2;  #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. */
     /*      /*
      * If J1 != J3, swap D[J1] and D[J3], and columns J1 and J3 of V.        Unscale the axes.
     */      */
     if (j1 != j3) {     if (scbd > 1.0) {
       t     = d[j1];  #ifdef DEBUGPRAX
       d[j1] = d[j3];        printf(" Unscale the axes.\n");
       d[j3] = t;  #endif
       for (i = 1; i <= n; i++) {        /* for (i=0; i<n; i++) { */
         t = v[i][j1];        for (i=1; i<=n; i++) {
         v[i][j1] = v[i][j3];            s = z[i];
         v[i][j3] = t;            /* for (j=0; j<n; j++) */
       } /* end i */            for (j=1; j<=n; j++)
     } /* end j1 != j3 */                v[i][j] *= s;
   } /* end j1 */        }
   return;        /* 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 */
   
 /* void svdcmp(double **a, int m, int n, double w[], double **v)  */  fret:
 void svdminfit(double **a, int m, int n, double w[])      if (prin > 0) {
 {           vecprint("\n  X:", x, n);
   /* From numerical recipes */           /* printf("\n... ChiSq reduced to %20.10e ...\n", fx); */
   /* Given a matrix a[1..m][1..n], this routine computes its singular value */           /* printf("... after %20u function calls.\n", nf); */
   /*   decomposition, A = U ·W ·V T . The matrix U replaces a on output. */     }
   /*   The diagonal matrix of singular values W is out- put as a vector w[1..n]. */     free_vector(d, 1, n);
   /*   The matrix V (not the transpose V T ) is output as v[1..n][1..n]. */     free_vector(y, 1, n);
        free_vector(z, 1, n);
   /* But in fact from Golub 1970 Algol60 */     free_vector(q0, 1, n);
      free_vector(q1, 1, n);
 /*   Computation of the singular values and complete orthogonal decom- */     free_matrix(v, 1, n, 1, n);
 /* position of a real rectangular matrix A, */     /*   double *d, *y, *z, */
 /* A = U diag (q) V^T, U^T U = V^T V =I , */     /* *q0, *q1, **v; */
 /*     where the arrays a[1:m, 1:n], u[1:m, 1 :n], v[1:n, 1:n], q[1:n] re- */     free_vector(tflin, 1, n);
 /* present A, U, V, q respectively. The actual parameters corresponding */     /* double *tflin; /\* used in flin: return (*fun)(tflin, n); *\/ */
 /* to a, u, v may all be identical unless withu=withv = true . In this */     free_vector(e, 1, n);
 /* case, the actual parameters corresponding to u and v must differ. */     /* double *e; /\* used in minfit, don't konw how to free memory and thus made global *\/ */
 /* m >= n is assumed; */     
        return(fx);
   /* Simplified (as in praxis) in order to output only V (replacing A), w is the diagonal matrix q  */  }
   
     
   double pythag(double a, double b);   
   int flag,i,its,j,jj,k,l,nm;   
   double anorm,c,f,g,h,s,scale,x,y,z,*rv1;   
     
   rv1=vector(1,n);  
   /* Householder's reduction to bidiagonal form; */  
   g=scale=anorm=0.0;   
   for (i=1;i<=n;i++) {   
     l=i+1;   
     rv1[i]=scale*g;   
     g=s=scale=0.0;   
     if (i <= m) {   
       for (k=i;k<=m;k++) scale += fabs(a[k][i]);   
       if (scale) {   
         for (k=i;k<=m;k++) {   
           a[k][i] /= scale;   
           s += a[k][i]*a[k][i];   
         }   
         f=a[i][i];   
         g = -SIGN(sqrt(s),f);   
         h=f*g-s;   
         a[i][i]=f-g;   
         for (j=l;j<=n;j++) {   
           for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j];   
           f=s/h;   
           for (k=i;k<=m;k++) a[k][j] += f*a[k][i];   
         }   
         for (k=i;k<=m;k++) a[k][i] *= scale;   
       }   
     }   
     w[i]=scale *g; /* p 411*/  
     g=s=scale=0.0;   
     if (i <= m && i != n) {   
       for (k=l;k<=n;k++) scale += fabs(a[i][k]);   
       if (scale) {   
         for (k=l;k<=n;k++) {   
           a[i][k] /= scale;   
           s += a[i][k]*a[i][k];   
         }   
         f=a[i][l];   
         g = -SIGN(sqrt(s),f);   
         h=f*g-s;   
         a[i][l]=f-g;   
         for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;   
         for (j=l;j<=m;j++) {   
           for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k];   
           for (k=l;k<=n;k++) a[j][k] += s*rv1[k];   
         }   
         for (k=l;k<=n;k++) a[i][k] *= scale;   
       }   
     }  
     /* y : = abs(q[i]) +abs(e[i]); if y > x then x : = y */  
     /* anorm=DMAX(anorm,(fabs(w[i])+fabs(rv1[i])));  */  
     anorm=FMAX(anorm,(fabs(w[i])+fabs(rv1[i])));   
   }  
   /* Comment accumulation of right-hand transformations */  
   for (i=n;i>=1;i--) {   
     if (i < n) {   
       if (g) {   
         for (j=l;j<=n;j++) a[j][i]=(a[i][j]/a[i][l])/g; /* Double division to avoid possible underflow. */  
         /* for (j=l;j<=n;j++) v[j][i]=(a[i][j]/a[i][l])/g;  */  
         for (j=l;j<=n;j++) {   
           /* for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j];  */  
            for (s=0.0,k=l;k<=n;k++) s += a[i][k]*a[k][j];  
           /* for (k=l;k<=n;k++) v[k][j] += s*v[k][i];  */  
           for (k=l;k<=n;k++) a[k][j] += s*a[k][i];   
         }   
       }   
       /* for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;  */  
       for (j=l;j<=n;j++) a[i][j]=a[j][i]=0.0;   
     }   
     /* v[i][i]=1.0;  */  
     a[i][i]=1.0;   
     g=rv1[i];   
     l=i;   
   }   
  /* Comment accumulation of left-hand transformations; */  
  /* for (i=IMIN(m,n);i>=1;i--) {  */  
  for (i=FMIN(m,n);i>=1;i--) {   
     l=i+1;   
     g=w[i];   
     for (j=l;j<=n;j++) a[i][j]=0.0;   
     if (g) {   
       g=1.0/g;   
       for (j=l;j<=n;j++) {   
         for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j];   
         f=(s/a[i][i])*g;   
         for (k=i;k<=m;k++) a[k][j] += f*a[k][i];   
       }   
       for (j=i;j<=m;j++) a[j][i] *= g;   
     } else for (j=i;j<=m;j++) a[j][i]=0.0;   
     ++a[i][i];   
   }   
 /* comment diagonalization of the bidiagonal form; */  
  for (k=n;k>=1;k--) {   
     for (its=1;its<=30;its++) { /* Loop over singular values, and over allowed iterations. */  
       flag=1;   
       for (l=k;l>=1;l--) { /* Test for splitting. */  
         nm=l-1;            /* Note that rv1[1] is always zero. */  
         if ((double)(fabs(rv1[l])+anorm) == anorm) {   
           flag=0;   
           break;   
         }   
         if ((double)(fabs(w[nm])+anorm) == anorm) break;   
       }   
       if (flag) {   
         c=0.0;    /* Cancellation of rv1[l], if l > 1. */  
         s=1.0;   
         for (i=l;i<=k;i++) {   
           f=s*rv1[i];   
           rv1[i]=c*rv1[i];   
           if ((double)(fabs(f)+anorm) == anorm) break;   
           g=w[i];   
           h=pythag(f,g);   
           w[i]=h;   
           h=1.0/h;   
           c=g*h;   
           s = -f*h;   
           for (j=1;j<=m;j++) {   
             y=a[j][nm];   
             z=a[j][i];   
             a[j][nm]=y*c+z*s;   
             a[j][i]=z*c-y*s;   
           }   
         }   
       }   
       z=w[k];   
       if (l == k) { /* Convergence */  
         if (z < 0.0) { /* Singular value is made nonnegative. */     
           w[k] = -z;     
           /* for (j=1;j<=n;j++) v[j][k] = -v[j][k];  */  
           for (j=1;j<=n;j++) a[j][k] = -a[j][k];   
         }   
         break;   
       }   
       if (its == 30) nrerror("no convergence in 30 svdcmp iterations");   
       x=w[l];  /* shift from bottom 2 x 2 minor; */  
       nm=k-1;   
       y=w[nm];   
       g=rv1[nm];   
       h=rv1[k];   
       f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);   
       g=pythag(f,1.0);   
       /* g=dpythag(f,1.0);  */  
       f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;   
       c=s=1.0;  /* next QR transformation  */  
       for (j=l;j<=nm;j++) {   
         i=j+1;   
         g=rv1[i];   
         y=w[i];   
         h=s*g;   
         g=c*g;   
         /* z=dpythag(f,h);  */  
         z=pythag(f,h);   
         rv1[j]=z;   
         c=f/z;   
         s=h/z;   
         f=x*c+g*s;   
         g = g*c-x*s;   
         h=y*s;   
         y *= c;   
 /* if withv then for j:= 1 step 1 until n do */  
         for (jj=1;jj<=n;jj++) {   
           /* x=v[jj][j];  */  
           /* z=v[jj][i];  */  
           x=a[jj][j];   
           z=a[jj][i];   
           /* v[jj][j]=x*c+z*s;  */  
           /* v[jj][i]=z*c-x*s;  */  
           a[jj][j]=x*c+z*s;   
           a[jj][i]=z*c-x*s;   
         }   
         /* z=dpythag(f,h);  */  
         z=pythag(f,h);   
         w[j]=z;   
         if (z) {   
           z=1.0/z;   
           c=f*z;   
           s=h*z;   
         }   
         f=c*g+s*y;   
         x=c*y-s*g;  
         /* if withu then for j:=1 step 1 until m do */  
         /* for (jj=1;jj<=m;jj++) {  */  
         /*   y=a[jj][j];  */  
         /*   z=a[jj][i];  */  
         /*   a[jj][j]=y*c+z*s;  */  
         /*   a[jj][i]=z*c-y*s;  */  
         /* }  */  
       }   
       rv1[l]=0.0;   
       rv1[k]=f;   
       w[k]=x;   
     }   
   }   
   free_vector(rv1,1,n);   
 }   
   
 /* end praxis */  /* end praxis gegen */
   
 /*************** powell ************************/  /*************** powell ************************/
 /*  /*
Line 3106  void powell(double p[], double **xi, int Line 4257  void powell(double p[], double **xi, int
       }        }
     }      }
     for (i=1;i<=n;i++) { /* For each direction i, maximisation after loading directions */      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 3115  void powell(double p[], double **xi, int Line 4267  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); /* New point i minimizing in direction xit i has coordinates p[j].*/        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.*/        /* xit[j] gives the n coordinates of direction i as input.*/
       /* *fret gives the maximum value on direction xit */        /* *fret gives the maximum value on direction xit */
 #else  #else
Line 3149  void powell(double p[], double **xi, int Line 4301  void powell(double p[], double **xi, int
     } /* end loop on each direction i */      } /* end loop on each direction i */
     /* Convergence test will use last linmin estimation (fret) and compare to 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) */
     for(j=1;j<=n;j++) {      for(j=1;j<=n;j++) {
       if(flatdir[j] >0){        if(flatdir[j] >0){
         printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);          printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
Line 3171  void powell(double p[], double **xi, int Line 4324  void powell(double p[], double **xi, int
       /* the scales of the directions and the directions, because the are reset to canonical directions */        /* the scales of the directions and the directions, because the are reset to canonical directions */
       /* Thus the first calls to linmin will give new points and better maximizations until fp-(*fret) is */        /* Thus the first calls to linmin will give new points and better maximizations until fp-(*fret) is */
       /* under the tolerance value. If the tolerance is very small 1.e-9, it could last long.  */        /* under the tolerance value. If the tolerance is very small 1.e-9, it could last long.  */
   #ifdef DEBUG
         int k[2],l;
         k[0]=1;
         k[1]=-1;
         printf("Max: %.12e",(*func)(p));
         fprintf(ficlog,"Max: %.12e",(*func)(p));
         for (j=1;j<=n;j++) {
           printf(" %.12e",p[j]);
           fprintf(ficlog," %.12e",p[j]);
         }
         printf("\n");
         fprintf(ficlog,"\n");
         for(l=0;l<=1;l++) {
           for (j=1;j<=n;j++) {
             ptt[j]=p[j]+(p[j]-pt[j])*k[l];
             printf("l=%d j=%d ptt=%.12e, xits=%.12e, p=%.12e, xit=%.12e", l,j,ptt[j],xits[j],p[j],xit[j]);
             fprintf(ficlog,"l=%d j=%d ptt=%.12e, xits=%.12e, p=%.12e, xit=%.12e", l,j,ptt[j],xits[j],p[j],xit[j]);
           }
           printf("func(ptt)=%.12e, deriv=%.12e\n",(*func)(ptt),(ptt[j]-p[j])/((*func)(ptt)-(*func)(p)));
           fprintf(ficlog,"func(ptt)=%.12e, deriv=%.12e\n",(*func)(ptt),(ptt[j]-p[j])/((*func)(ptt)-(*func)(p)));
         }
   #endif
   
       free_vector(xit,1,n);         free_vector(xit,1,n); 
       free_vector(xits,1,n);         free_vector(xits,1,n); 
Line 3179  void powell(double p[], double **xi, int Line 4354  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)=2Pn-P0 */      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]; /* Coordinate j of last direction xi_n=P_n-_0 */        xit[j]=p[j]-pt[j]; /* Coordinate j of last direction xi_n=P_n-P_0 */
   #ifdef DEBUG
       printf("\n %d xit=%12.7g p=%12.7g pt=%12.7g ",j,xit[j],p[j],pt[j]);        printf("\n %d xit=%12.7g p=%12.7g pt=%12.7g ",j,xit[j],p[j],pt[j]);
       pt[j]=p[j];   #endif
         pt[j]=p[j]; /* New P0 is Pn */
     }      }
   #ifdef DEBUG
     printf("\n");      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 3217  void powell(double p[], double **xi, int Line 4396  void powell(double p[], double **xi, int
       t= t- del*SQR(fp-fptt);        t= t- del*SQR(fp-fptt);
 #endif  #endif
       directest = fp-2.0*(*fret)+fptt - 2.0 * del; /* If delta was big enough we change it for a new direction */        directest = fp-2.0*(*fret)+fptt - 2.0 * del; /* If delta was big enough we change it for a new direction */
       printf(" t=%g, directest=%g\n",t, directest);  #ifdef DEBUG
 #ifdef POWELLNOTTRUECONJUGATE   /* Searching for IBIG and testing for replacement */          printf("t1= %.12lf, t2= %.12lf, t=%.12lf  directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest);
         fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest);
         printf("t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt),
                (fp-(*fret)-del)*(fp-(*fret)-del),(fp-fptt)*(fp-fptt));
         fprintf(ficlog,"t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt),
                (fp-(*fret)-del)*(fp-(*fret)-del),(fp-fptt)*(fp-fptt));
         printf("tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t);
         fprintf(ficlog, "tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t);
   #endif
 #ifdef POWELLORIGINAL  #ifdef POWELLORIGINAL
       if (t < 0.0) { /* Then we use it for new direction */        if (t < 0.0) { /* Then we use it for new direction */
 #else  /* Not POWELLOriginal but Brouard's */  #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);
       }         } 
       if (directest < 0.0) { /* Then we use (P0, Pn) for new direction Xi_n or Xi_iBig */        if (directest < 0.0) { /* Then we use it for new direction */
 #endif /* end POWELLOriginal */  #endif
 #endif  /* POWELLNOTTRUECONJUGATE  else means systematic replacement by new direction P_0P_n */  #ifdef DEBUGLINMIN
           printf("Before linmin in direction P%d-P0\n",n);
           for (j=1;j<=n;j++) {
             printf(" Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
             fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
             if(j % ncovmodel == 0){
               printf("\n");
               fprintf(ficlog,"\n");
             }
           }
   #endif
 #ifdef LINMINORIGINAL  #ifdef LINMINORIGINAL
         /* xit[j]=p[j]-pt[j] */          linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
         printf("\n Computes min on P_0, P_n direction iter=%d Bigter=%d\n",*iter,Bigter);  #else
         linmin(p,xit,n,fret,func); /* computes minimum on P_0,P_n direction: changes p and rescales xit.*/  
 #else  /* NOT LINMINORIGINAL but with searching for flat directions */  
         printf("\n Flat Computes min on P_0, P_n direction iter=%d Bigter=%d\n",*iter,Bigter);  
         linmin(p,xit,n,fret,func,&flat); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/          linmin(p,xit,n,fret,func,&flat); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
         flatdir[i]=flat; /* Function is vanishing in that direction i */          flatdir[i]=flat; /* Function is vanishing in that direction i */
 #endif  #endif
                   
 #ifdef POWELLNOTTRUECONJUGATE  #ifdef DEBUGLINMIN
 #else  
 #ifdef POWELLORIGINCONJUGATE  
         for (j=1;j<=n;j++) {           for (j=1;j<=n;j++) { 
           xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */            printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
           xi[j][n]=xit[j];      /* and this nth direction by the by the average p_0 p_n */            fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
         }            if(j % ncovmodel == 0){
 #else              printf("\n");
         for (i=1;i<=n-1;i++) {               fprintf(ficlog,"\n");
           for (j=1;j<=n;j++) {   
             xi[j][i]=xi[j][i+1]; /* Standard method of conjugate directions, not Powell who changes the nth direction by p0 pn . */  
           }            }
         }          }
   #endif
         for (j=1;j<=n;j++) {           for (j=1;j<=n;j++) { 
             xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */
           xi[j][n]=xit[j];      /* and this nth direction by the by the average p_0 p_n */            xi[j][n]=xit[j];      /* and this nth direction by the by the average p_0 p_n */
         }          }
 #endif  /* POWELLORIGINCONJUGATE*/  
 #endif  /*POWELLNOTTRUECONJUGATE*/  
         printf(" Standard method of conjugate directions\n");  
         printf("\n#A Before prax Bigter=%d model=  1      +     age ", Bigter);  
         for(j=1;j<=n;j++){  
           printf("%d \n",j);  
           for(i=1;i<=n;i++){  
             printf(" %f",xi[j][i]);  
           }  
         }  
         printf("\n");  
   
 #ifdef NOTMINFIT  
 #else   
         if(*iter >n){  
         /* if(Bigter >n){ */  
           printf("\n#Before prax Bigter=%d model=  1      +     age ", Bigter);  
           printf("\n");  
           for(j=1;j<=n;j++){  
             printf("%d \n",j);  
             for(i=1;i<=n;i++){  
               printf(" %f",xi[j][i]);  
             }  
           }  
           printf("\n");  
           /*  
            * Calculate a new set of orthogonal directions before repeating  
            * the main loop.  
            * Transpose V for SVD (minfit) (because minfit returns the right V in ULV=A):  
            */  
           printf(" Bigter=%d Calculate a new set of orthogonal directions before repeating  the main loop.\n  Transpose V for MINFIT:...\n",Bigter);  
           transpose_in_place ( n, xi );  
           /*  
             SVD/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.  
           */  
           printf(" SVDMINFIT 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");  
           double *d; /* eigenvalues of principal directions */  
           d=vector(1,n);  
             
             
           svdminfit (xi, n, n, d ); /* In Brent's notation find d such that V=Q Diagonal(d) R, and Lambda=d^(-1/2) */  
           
           printf("\n#After prax model=  1      +     age ");  
           fprintf(ficlog,"\n#model=  1      +     age ");  
             
           if(nagesqr==1){  
             printf("  + age*age  ");  
             fprintf(ficlog,"  + age*age  ");  
           }  
           for(j=1;j <=ncovmodel-2;j++){  
             if(Typevar[j]==0) {  
               printf("  +      V%d  ",Tvar[j]);  
               fprintf(ficlog,"  +      V%d  ",Tvar[j]);  
             }else if(Typevar[j]==1) {  
               printf("  +    V%d*age ",Tvar[j]);  
               fprintf(ficlog,"  +    V%d*age ",Tvar[j]);  
             }else if(Typevar[j]==2) {  
               printf("  +    V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);  
               fprintf(ficlog,"  +    V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);  
             }else if(Typevar[j]==3) {  
               printf("  +    V%d*V%d*age ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);  
               fprintf(ficlog,"  +    V%d*V%d*age ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);  
             }  
           }  
           printf("\n");  
           fprintf(ficlog,"\n");  
           for(i=1,jk=1; i <=nlstate; i++){  
             for(k=1; k <=(nlstate+ndeath); k++){  
               if (k != i) {  
                 printf("%d%d ",i,k);  
                 fprintf(ficlog,"%d%d ",i,k);  
                 for(j=1; j <=ncovmodel; j++){  
                   printf("%12.7f ",p[jk]);  
                   fprintf(ficlog,"%12.7f ",p[jk]);  
                   jk++;   
                 }  
                 printf("\n");  
                 fprintf(ficlog,"\n");  
               }  
             }  
           }  
           /* minfit ( n, vsmall, v, d ); */  
           /* v is overwritten with R. */  
           /*  
             Heuristic numbers:  
             If the axes may be badly scaled (which is to be avoided if  
             possible), then set SCBD = 10.  Otherwise set SCBD = 1.  
               
             If the problem is known to be ill-conditioned, initialize ILLC = true.  
             KTM is the number of iterations without improvement before the  
             algorithm terminates.  KTM = 4 is very cautious; usually KTM = 1  
             is satisfactory.  
           */  
           double machep, small;  
           double dmin;  
           int illc=0; /*    Local, int ILLC, is TRUE if the system is ill-conditioned. */  
           machep = DBL_EPSILON;  
           small = machep * machep;  
           /* m2 = dsqrt(machep); */  
             
           /*  
            * Sort the eigenvalues and eigenvectors.  
            */  
           printf(" Sort the eigenvalues and eigenvectors....\n");  
           svsort ( n, d, xi );  
           printf("Sorted Eigenvalues:\n");  
           for(i=1; i<=n;i++){  
             printf(" d[%d]=%g",i,d[i]);  
           }  
           printf("\n");  
           /*  
            * Determine the smallest eigenvalue.  
            */  
           printf("  Determine the smallest eigenvalue.\n");  
           dmin = fmax ( d[n], small );  
           /*  
            * The ratio of the smallest to largest eigenvalue determines whether  
            * the system is ill conditioned.  
            */  
           if ( dmin < sqrt(machep) * d[1] ) {  /* m2*d[0] */  
             illc = 1;  
           } else  {  
             illc = 0;  
           }  
           printf("  The ratio of the smallest to largest eigenvalue determines whether\n \  
   the system is ill conditioned=%d . dmin=%.12lf < m2=%.12lf * d[1]=%.12lf \n",illc, dmin,sqrt(machep), d[1]);  
           /*    if ( 1.0 < scbd ) { */  
           /*      r8vec_print ( n, z, "  The scale factors:" ); */  
           /*    }  */  
           /*    r8vec_print ( n, d, "  Principal values of the quadratic form:" ); */  
           /* } */  
           /* if ( 3 < prin ) { */  
           /*    r8mat_print ( n, n, v, "  The principal axes:" ); */  
           /* } */  
           free_vector(d,1,n);     
           /*  
            * The main loop ends here.  
            */  
             
           /* if ( 0 < prin ) */  
           /* { */  
           /*   r8vec_print ( n, x, "  X:" ); */  
           /* } */  
         }         
 #endif  /* NOTMINFIT */  
 #ifdef LINMINORIGINAL  #ifdef LINMINORIGINAL
 #else  #else
         for (j=1, flatd=0;j<=n;j++) {          for (j=1, flatd=0;j<=n;j++) {
Line 3429  void powell(double p[], double **xi, int Line 4473  void powell(double p[], double **xi, int
           free_vector(pt,1,n);             free_vector(pt,1,n); 
           return;            return;
 #endif  #endif
         }  /* endif(flatd >0) */          }
 #endif  #endif
         printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);          printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
         fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);          fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
   /* The minimization in direction $\xi_1$ gives $P_1$. From $P_1$ minimization in direction $\xi_2$ gives */          
   /* $P_2$. Minimization of line $P_2-P_1$ gives new starting point $P^{(1)}_0$ and direction $\xi_1$ is dropped and replaced by second */  #ifdef DEBUG
   /* direction $\xi_1^{(1)}=\xi_2$. Also second direction is replaced by new direction $\xi^{(1)}_2=P_2-P_0$. */          printf("Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
           fprintf(ficlog,"Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
   /* At the second iteration, starting from $P_0^{(1)}$, minimization amongst $\xi^{(1)}_1$ gives point $P^{(1)}_1$. */          for(j=1;j<=n;j++){
   /* Minimization amongst $\xi^{(1)}_2=(P_2-P_0)$ gives point $P^{(1)}_2$.  As $P^{(2)}_1$ and */            printf(" %lf",xit[j]);
   /* $P^{(1)}_0$ are minimizing in the same direction $P^{(1)}_2 - P^{(1)}_1= P_2-P_0$, directions $P^{(1)}_2-P^{(1)}_0$ */            fprintf(ficlog," %lf",xit[j]);
   /* and $P_2-P_0$ (parallel to $\xi$ and $\xi^c$) are conjugate.  } */          }
 #ifdef POWELLNOTTRUECONJUGATE            printf("\n");
       } /* end of t or directest negative */          fprintf(ficlog,"\n");
 #endif  #endif
         } /* end of t or directest negative */
         printf(" Directest is positive, P_n-P_0 does not increase the conjugacy. n=%d\n",n);
         fprintf(ficlog," Directest is positive, P_n-P_0 does not increase the conjugacy. n=%d\n",n);
 #ifdef POWELLNOF3INFF1TEST  #ifdef POWELLNOF3INFF1TEST
 #else  #else
       } /* end if (fptt < fp)  */        } /* end if (fptt < fp)  */
 #endif  #endif
 #ifdef POWELLORIGINCONJUGATE  
     } /* if (t < 0.0) */  
 #endif  
 #ifdef NODIRECTIONCHANGEDUNTILNITER  /* No change in drections until some iterations are done */  #ifdef NODIRECTIONCHANGEDUNTILNITER  /* No change in drections until some iterations are done */
     } /*NODIRECTIONCHANGEDUNTILNITER  No change in drections until some iterations are done */      } /*NODIRECTIONCHANGEDUNTILNITER  No change in drections until some iterations are done */
 #else  #else
Line 3655  void powell(double p[], double **xi, int Line 4699  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 3690  void powell(double p[], double **xi, int Line 4736  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 3957  double **pmij(double **ps, double *cov, Line 5003  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 4172  double ***hpxij(double ***po, int nhstep Line 5218  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 4353  double ***hbxij(double ***po, int nhstep Line 5399  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 4876  double func( double *x) Line 5922  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 5426  void likelione(FILE *ficres,double p[], Line 6472  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 */    /*double * p1;*/ /* Shifted parameters from 0 instead of 1 */
 #ifdef NLOPT  #ifdef NLOPT
   int creturn;    int creturn;
   nlopt_opt opt;    nlopt_opt opt;
Line 5512  void mlikeli(FILE *ficres,double p[], in Line 6558  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 ); */  /*   praxis ( t0, h0, n, prin, x, beale_f ); */
 /*   int prin=4; */    int prin=1;
 /*   double h0=0.25; */    double h0=0.25;
     double macheps;
     double fmin;
     macheps=pow(16.0,-13.0);
 /* #include "praxis.h" */  /* #include "praxis.h" */
   /* Be careful that praxis start at x[0] and powell start at p[1] */    /* Be careful that praxis start at x[0] and powell start at p[1] */
    /* praxis ( ftol, h0, npar, prin, p, func ); */     /* praxis ( ftol, h0, npar, prin, p, func ); */
 /* p1= (p+1); /\*  p *(p+1)@8 and p *(p1)@8 are equal p1[0]=p[1] *\/ */  /* p1= (p+1); */ /*  p *(p+1)@8 and p *(p1)@8 are equal p1[0]=p[1] */
 /* printf("Praxis \n"); */  printf("Praxis Gegenfurtner \n");
 /* fprintf(ficlog, "Praxis \n");fflush(ficlog); */  fprintf(ficlog, "Praxis  Gegenfurtner\n");fflush(ficlog);
 /* praxis ( ftol, h0, npar, prin, p1, func ); */  /* praxis ( ftol, h0, npar, prin, p1, func ); */
 /* printf("End Praxis\n"); */    /* 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 6654  void prevalence(double ***probs, double Line 7705  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 8541  void printinghtml(char fileresu[], char Line 9592  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 10405  void prevforecast(char fileres[], double Line 11456  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 10556  void prevforecast(char fileres[], double Line 11607  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 11242  void printinggnuplotmort(char fileresu[] Line 12293  void printinggnuplotmort(char fileresu[]
   
   char dirfileres[132],optfileres[132];    char dirfileres[132],optfileres[132];
   
   int ng;    /*int ng;*/
   
   
   /*#ifdef windows */    /*#ifdef windows */
Line 11266  int readdata(char datafile[], int firsto Line 12317  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 11883  int decodemodel( char model[], int lasto Line 12934  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 13311  int hPijx(double *p, int bage, int fage) Line 14362  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 13517  int main(int argc, char *argv[]) Line 14568  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;
   
Line 13530  int main(int argc, char *argv[]) Line 14581  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 13565  int main(int argc, char *argv[]) Line 14616  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 13623  int main(int argc, char *argv[]) Line 14674  int main(int argc, char *argv[])
   getcwd(pathcd, size);    getcwd(pathcd, size);
 #endif  #endif
   syscompilerinfo(0);    syscompilerinfo(0);
   printf("\nIMaCh prax version minfit %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 14525  This file: <a href=\"%s\">%s</a></br>Tit Line 15576  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 15835  Please run with mle=-1 to get a correct Line 16886  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.5  
changed lines
  Added in v.1.6

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