#! /usr/bin/env python
#
def mckinnon ( m, x ):

#*****************************************************************************80
#
## MCKINNON computes the McKinnon function.
#
#  Discussion:
#
#    This function has a global minimizer:
#
#      X* = ( 0.0, -0.5 ), F(X*) = -0.25.
#
#    There are three parameters, TAU, THETA and PHI.
#
#    1 < TAU, then F is strictly convex.
#             and F has continuous first derivatives.
#    2 < TAU, then F has continuous second derivatives.
#    3 < TAU, then F has continuous third derivatives.
#
#    However, this function can cause the Nelder-Mead optimization
#    algorithm to "converge" to a point which is not the minimizer
#    of the function F.
#
#    Sample parameter values which cause problems for Nelder-Mead 
#    include:
#
#      PHI = 10,  TAU = 1, THETA = 15
#      PHI = 60,  TAU = 2, THETA =  6
#      PHI = 400, TAU = 3, THETA =  6
#
#    To get the bad behavior, we also assume the initial simplex has the form
#
#      X1 = (0,0),
#      X2 = (1,1),
#      X3 = (A,B), 
#
#    where 
#
#      A = (1+sqrt(33))/8 =  0.84307...
#      B = (1-sqrt(33))/8 = -0.59307...
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    23 January 2016
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Ken McKinnon,
#    Convergence of the Nelder-Mead simplex method to a nonstationary point,
#    SIAM Journal on Optimization,
#    Volume 9, Number 1, 1998, pages 148-158.
#
#  Parameters:
#
#    Input, integer M, the number of variables.
#
#    Input, real X(M), the argument of the function.
#
#    Output, real F, the value of the function at X.
#
  global phi
  global tau
  global theta

  if ( x[0] <= 0.0 ):
    f = theta * phi * abs ( x[0] ) ** tau + x[1] * ( 1.0 + x[1] )
  else:
    f = theta       *       x[0] ** tau   + x[1] * ( 1.0 + x[1] )

  return f

def mckinnon_test ( ):

#*****************************************************************************80
#
## MCKINNON_TEST works with the McKinnon function.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    23 January 2016
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from compass_search import compass_search
  from r8vec_print import r8vec_print

  global phi
  global tau
  global theta

  print ( '' )
  print ( 'MCKINNON_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Test COMPASS_SEARCH with the McKinnon function.' )
  m = 2
  delta_tol = 0.00001
  delta = 0.3
  k_max = 20000
#
#  Test 1
#
  a = ( 1.0 + np.sqrt ( 33.0 ) ) / 8.0
  b = ( 1.0 - np.sqrt ( 33.0 ) ) / 8.0

  phi = 10.0
  tau = 1.0
  theta = 15.0

  x = np.array ( [ a, b ] )
  r8vec_print ( m, x, '  Initial point X0:' )
  print ( '  PHI = %g, TAU = %g, THETA = %g' % ( phi, tau, theta ) )
  print ( '' )
  print ( '  F(X0) = %g' % ( mckinnon ( m, x ) ) )

  x, fx, k = compass_search ( mckinnon, m, x, delta_tol, delta, k_max )
  r8vec_print ( m, x, '  Estimated minimizer X1:' )
  print ( '' )
  print ( '  F(X1) = %g, number of steps = %d' % ( fx, k ) )

  x = np.array ( [ 0.0, -0.5 ] )
  r8vec_print ( m, x, '  Correct minimizer X*:' )
  print ( '' )
  print ( '  F(X*) = %g' % ( mckinnon ( m, x ) ) )
#
#  Test 2
#
  a = ( 1.0 + np.sqrt ( 33.0 ) ) / 8.0
  b = ( 1.0 - np.sqrt ( 33.0 ) ) / 8.0

  phi = 60.0
  tau = 2.0
  theta = 6.0

  x = np.array ( [ a, b ] )
  r8vec_print ( m, x, '  Initial point X0:' )
  print ( '  PHI = %g, TAU = %g, THETA = %g' % ( phi, tau, theta ) )
  print ( '' )
  print ( '  F(X0) = %g' % ( mckinnon ( m, x ) ) )

  x, fx, k = compass_search ( mckinnon, m, x, delta_tol, delta, k_max )
  r8vec_print ( m, x, '  Estimated minimizer X1:' )
  print ( '' )
  print ( '  F(X1) = %g, number of steps = %d' % ( fx, k ) )

  x = np.array ( [ 0.0, -0.5 ] )
  r8vec_print ( m, x, '  Correct minimizer X*:' )
  print ( '' )
  print ( '  F(X*) = %g' % ( mckinnon ( m, x ) ) )
#
#  Test 3
#
  a = ( 1.0 + np.sqrt ( 33.0 ) ) / 8.0
  b = ( 1.0 - np.sqrt ( 33.0 ) ) / 8.0

  phi = 4000.0
  tau = 3.0
  theta = 6.0

  x = np.array ( [ a, b ] )
  r8vec_print ( m, x, '  Initial point X0:' )
  print ( '  PHI = %g, TAU = %g, THETA = %g' % ( phi, tau, theta ) )
  print ( '' )
  print ( '  F(X0) = %g' % ( mckinnon ( m, x ) ) )

  x, fx, k = compass_search ( mckinnon, m, x, delta_tol, delta, k_max )
  r8vec_print ( m, x, '  Estimated minimizer X1:' )
  print ( '' )
  print ( '  F(X1) = %g, number of steps = %d' % ( fx, k ) )

  x = np.array ( [ 0.0, -0.5 ] )
  r8vec_print ( m, x, '  Correct minimizer X*:' )
  print ( '' )
  print ( '  F(X*) = %g' % ( mckinnon ( m, x ) ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'MCKINNON_TEST' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  mckinnon_test ( )
  timestamp ( )