function [ value, x ] = praxis ( t0, h0, n, prin, x, f )

%*****************************************************************************80
%
%% praxis() seeks an N-dimensional minimizer X of a scalar function F(X).
%
%  Discussion:
%
%    PRAXIS returns the minimum of the function F(X,N) of N variables
%    using the principal axis method.  The gradient of the function is
%    not required.
%
%    The approximating quadratic form is
%
%      Q(x') = F(x,n) + (1/2) * (x'-x)' * A * (x'-x)
%
%    where X is the best estimate of the minimum and
%
%      A = inverse(V') * D * inverse(V)
%
%    V(*,*) is the matrix of search directions;
%    D(*) is the array of second differences.
%
%    If F(X) has continuous second derivatives near X0, then A will tend
%    to the hessian of F at X0 as X approaches X0.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    31 July 2016
%
%  Author:
%
%    Original FORTRAN77 version by Richard Brent.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Richard Brent,
%    Algorithms for Minimization with Derivatives,
%    Prentice Hall, 1973,
%    Reprinted by Dover, 2002.
%
%  Input:
%
%    real T0, is a tolerance.  PRAXIS attempts to return
%    praxis = f(x) such that if X0 is the true local minimum near X, then
%    norm ( x - x0 ) < T0 + sqrt ( EPSILON ( X ) ) * norm ( X ),
%    where EPSILON ( X ) is the machine precision for X.
%
%    real H0, is the maximum step size.  H0 should be
%    set to about the maximum distance from the initial guess to the minimum.
%    If H0 is set too large or too small, the initial rate of
%    convergence may be slow.
%
%    integer N, the number of variables.
%
%    integer PRIN, controls printing intermediate results.
%    0, nothing is printed.
%    1, F is printed after every n+1 or n+2 linear minimizations.
%       final X is printed, but intermediate X is printed only
%       if N is at most 4.
%    2, the scale factors and the principal values of the approximating
%       quadratic form are also printed.
%    3, X is also printed after every few linear minimizations.
%    4, the principal vectors of the approximating quadratic form are
%       also printed.
%
%    real X(N),  a guess of the point of minimum.
%
%    real F ( X, N ), the function to be minimized.
%
%  Output:
%
%    real VALUE, the function value at the minimizer.
%
%    real X(N), the estimated point of minimum.
%
%  Local:
%
%    real DMIN, an estimate for the smallest eigenvalue.
%
%    real FX, the value of F(X,N).
%
%    logical ILLC, is TRUE if the system is ill-conditioned.
%
%    real LDT, the length of the step.
%
%    integer NF, the number of function evaluations.
%
%    integer NL, the number of linear searches.
%
  x = x(:);
%
%  Initialization.
%
  machep = eps;
  small = machep * machep;
  vsmall = small * small;
  large = 1.0 / small;
  vlarge = 1.0 / vsmall;
  m2 = sqrt ( machep );
  m4 = sqrt ( m2 );
%
%  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.
%
  scbd = 1.0;
  illc = false;
  ktm = 1;

  if ( illc )
    ldfac = 0.1;
  else
    ldfac = 0.01;
  end

  kt = 0;
  nl = 0;
  nf = 1;
  fx = f ( x, n );
  qf1 = fx;
  t = small + abs ( t0 );
  t2 = t;
  dmin = small;
  h = h0;
  h = max ( h, 100.0 * t );
  ldt = h;
%
%  The initial set of search directions V is the identity matrix.
%
  v = zeros ( n, n );
  for i = 1 : n
    v(i,i) = 1.0;
  end

  d = zeros ( n, 1 );
  y = zeros ( n, 1 );
  z = zeros ( n, 1 );
  qa = 0.0;
  qb = 0.0;
  qc = 0.0;
  qd0 = 0.0;
  qd1 = 0.0;
  q0 = x;
  q1 = x;

  if ( 0 < prin )
    print2 ( n, x, prin, fx, nf, nl );
  end
%
%  The main loop starts here.
%
  while ( true )

    sf = d(1);
    d(1) = 0.0;
    s = 0.0;
%
%  Minimize along the first direction V(*,1).
%
    jsearch = 1;
    nits = 2;
    d2 = d(1);
    value = fx;
    fk = false;

    [ d2, s, value, x, nl, nf, fx, qa, qb, qc ] = minny ( n, jsearch, nits, ...
      d2, s, value, fk, f, x, t, h, v, q0, q1, nl, nf, dmin, ldt, ...
      fx, qa, qb, qc, qd0, qd1 );

    d(1) = d2;

    if ( s <= 0.0 )
      v(1:n,1) = - v(1:n,1);
    end

    if ( sf <= 0.9 * d(1) || d(1) <= 0.9 * sf )
      d(2:n) = 0.0;
    end
%
%  The inner loop starts here.
%
    for k = 2 : n

      y = x;

      sf = fx;

      if ( 0 < kt )
        illc = true;
      end

      while ( true )

        kl = k;
        df = 0.0;
%
%  A random step follows, to avoid resolution valleys.
%
        if ( illc )
          for i = 1 : n
            r = rand ( 1, 1 );
            s = ( 0.1 * ldt + t2 * 10.0 ^ kt ) * ( r - 0.5 );
            z(i) = s;
            x(1:n) = x(1:n) + s * v(1:n,i);
          end

          fx = f ( x, n );
          nf = nf + 1;

        end
%
%  Minimize along the "non-conjugate" directions V(*,K),...,V(*,N).
%
        for k2 = k : n

          sl = fx;

          jsearch = k2;
          nits = 2;
          d2 = d(k2);
          s = 0.0;
          value = fx;
          fk = false;

          [ d2, s, value, x, nl, nf, fx, qa, qb, qc ] = minny ( n, jsearch, nits, ...
            d2, s, value, fk, f, x, t, h, v, q0, q1, nl, nf, dmin, ldt, ...
            fx, qa, qb, qc, qd0, qd1 );

          d(k2) = d2;

          if ( illc )
            s = d(k2) * ( ( s + z(k2) ) ^ 2 );
          else
            s = sl - fx;
          end

          if ( df <= s )
            df = s;
            kl = k2;
          end

        end
%
%  If there was not much improvement on the first try, set
%  ILLC = true and start the inner loop again.
%
        if ( illc )
          break
        end

        if ( abs ( 100.0 * machep * fx ) <= df )
          break
        end

        illc = true;

      end

      if ( k == 2 && 1 < prin )
        r8vec_print ( n, d, '  The second difference array:' );
      end
%
%  Minimize along the "conjugate" directions V(*,1),...,V(*,K-1).
%
      for k2 = 1 : k - 1

        jsearch = k2;
        nits = 2;
        d2 = d(k2);
        s = 0.0;
        value = fx;
        fk = false;

        [ d2, s, value, x, nl, nf, fx, qa, qb, qc ] = minny ( n, jsearch, nits, ...
          d2, s, value, fk, f, x, t, h, v, q0, q1, nl, nf, dmin, ldt, ...
          fx, qa, qb, qc, qd0, qd1 );

        d(k2) = d2;

      end

      f1 = fx;
      fx = sf;
      lds = 0.0;

      for i = 1 : n
        sl = x(i);
        x(i) = y(i);
        sl = sl - y(i);
        y(i) = sl;
        lds = lds + sl * sl;
      end

      lds = sqrt ( lds );
%
%  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 ( small < lds )

        for j = kl - 1 : -1 : k
          v(1:n,j+1) = v(1:n,j);
          d(j+1) = d(j);
        end

        d(k) = 0.0;

        v(1:n,k) = y(1:n) / lds;
%
%  Minimize along the new "conjugate" direction V(*,k), which is
%  the normalized vector:  (new x) - (old x).
%
        jsearch = k;
        nits = 4;
        d2 = d(k);
        value = f1;
        fk = true;

        [ d2, lds, value, x, nl, nf, fx, qa, qb, qc ] = minny ( n, jsearch, nits, ...
          d2, lds, value, fk, f, x, t, h, v, q0, q1, nl, nf, dmin, ldt, ...
          fx, qa, qb, qc, qd0, qd1 );

        d(k) = d2;

        if ( lds <= 0.0 )
          lds = - lds;
          v(1:n,k) = - v(1:n,k);
        end

      end

      ldt = ldfac * ldt;
      ldt = max ( ldt, lds );

      if ( 0 < prin )
        print2 ( n, x, prin, fx, nf, nl );
      end

      t2 = m2 * norm ( x ) + t;
%
%  See whether the length of the step taken since starting the
%  inner loop exceeds half the tolerance.
%
      if ( 0.5 * t2 < ldt )
        kt = - 1;
      end

      kt = kt + 1;

      if ( ktm < kt )

        if ( 0 < prin )
          r8vec_print ( n, x, '  X:' );
        end

        value = fx;

        return

      end

    end
%
%  The inner loop ends here.
%
%  Try quadratic extrapolation in case we are in a curved valley.
%
    [ x, q0, q1, nl, nf, fx, qf1, qa, qb, qc, qd0, qd1 ] = quad ( ...
      n, f, x, t, h, v, q0, q1, nl, nf, dmin, ldt, fx, qf1, qa, qb, qc, qd0, qd1 );

    d(1:n) = 1.0 ./ sqrt ( d(1:n) );

    dn = max ( d(1:n) );

    if ( 3 < prin )
      r8mat_print ( n, n, v, '  The new direction vectors:' );
    end

    for j = 1 : n
      v(1:n,j) = ( d(j) / dn ) * v(1:n,j);
    end
%
%  Scale the axes to try to reduce the condition number.
%
    if ( 1.0 < scbd )

      for i = 1 : n
        z(i) = max ( m4, norm ( v(i,1:n) ) );
      end

      s = min ( z(1:n) );

      for i = 1 : n

        sl = s / z(i);
        z(i) = 1.0 / sl;

        if ( scbd < z(i) )
          sl = 1.0 / scbd;
          z(i) = scbd;
        end

        v(i,1:n) = sl * v(i,1:n);

      end

    end
%
%  Calculate a new set of orthogonal directions before repeating
%  the main loop.
%
%  Transpose V for MINFIT:
%
    v = v';
%
%  Call MINFIT to find the singular value decomposition of V.
%
%  This gives the principal values and principal directions of the
%  approximating quadratic form without squaring the condition number.
%
    [ v, d ] = minfit ( n, vsmall, v );
%
%  Unscale the axes.
%
    if ( 1.0 < scbd )

      for i = 1 : n
        v(i,1:n) = z(i) * v(i,1:n);
      end

      for i = 1 : n

        s = norm ( v(1:n,i) );
        d(i) = s * d(i);
        v(1:n,i) = v(1:n,i) / s;

      end

    end

    for i = 1 : n

      dni = dn * d(i);

      if ( large < dni )
        d(i) = vsmall;
      elseif ( dni < small )
        d(i) = vlarge;
      else
        d(i) = 1.0 / dni ^ 2;
      end

    end
%
%  Sort the singular values and singular vectors.
%
    [ d, v ] = svsort ( n, d, v );
%
%  Determine the smallest eigenvalue.
%
    dmin = max ( d(n), small );
%
%  The ratio of the smallest to largest eigenvalue determines whether
%  the system is ill conditioned.
%
    if ( dmin < m2 * d(1) )
      illc = true;
    else
      illc = false;
    end

    if ( 1 < prin )

      if ( 1.0 < scbd )
        r8vec_print ( n, z, '  The scale factors:' );
      end

      r8vec_print ( n, d, '  Principal values of the quadratic form:' );

    end

    if ( 3 < prin )
      r8mat_print ( n, n, v, '  The principal axes:' );
    end
%
%  The main loop ends here.
%
  end

  if ( 0 < prin )
    r8vec_print ( n, x, '  X:' );
  end

  value = fx;

  return
end
