#! /usr/bin/env python3
#
def quartic_gradient1 ( ):

#*****************************************************************************80
#
## quartic_gradient1() applies gradient descent to the quartic function.
#
#  Discussion:
#
#    f (x) = 2 * x**4 - 4 * x**2 + x + 20
#    f'(x) = 8 * x**3 - 8 * x + 1
#
#    We seek a minimizer of f(x) in the range -2 <= x <= 2.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    22 September 2019
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  print ( '' )
  print ( 'quartic_gradient1:' )
  print ( '  Seek minimizer of a quartic scalar function f(x)' )
#
#  Plot the function.
#
  xmin = -2.0
  xmax = +2.0;

  xp = np.linspace ( xmin, xmax, 101 )
  fp = quartic ( xp )
  plt.plot ( xp, fp, linewidth = 3 )
  plt.plot ( [ xmin, xmax ], [ 0, 0 ], 'k', linewidth = 3 )
  plt.grid ( True )
  plt.xlabel ( '<-- X -->' )
  plt.ylabel ( '<-- F(X) -->' )
  plt.title ( 'Quartic function', fontsize = 16 )
  filename = 'quartic.png'
  plt.savefig ( filename )
  print ( '' )
  print ( '  Graphics saved as "%s"' % ( filename ) )
  plt.show ( )
  plt.close ( )
#
#  Choose a random initial guess.
#
  for test in range ( 0, 2 ):
    r = 0.05

    it = 0
    i2 = 0
    terminate = False

    print ( "" )
    print ( "   it    x    f(x)    f'(x)" )
    while ( True ):

      if ( terminate ):
        break
      elif ( it == 0 ):
        if ( test == 0 ):
          x0 = -1.6
        else:
          x0 = 1.8
        x = x0
      elif ( it < 1000 ):

        beta = 1.0
        while ( True ):
          xn = x - beta * r * quartic_df ( x )
          it = it + 1
          if ( quartic ( xn ) < quartic ( x ) ):
            i2 = i2 + 1
            x = xn
            break
          beta = beta * 0.5
          if ( beta < 1.0E-10 ):
            terminate = True
            break

      else:
        break
      print ( '  %3d  %12.8g  %12.8g  %12.8g' % ( it, x, quartic ( x ), quartic_df ( x ) ) )
      it = it + 1

    print ( '' )
    print ( '  ', it, 'gradient descent steps were tried.' )
    print ( '  ', i2, 'steps were successful.' )
#
#  Plot the function with initial point and local minimizer
#
    xvec = np.linspace ( xmin, xmax, 101 )
    fxvec = quartic ( xvec )
    plt.plot ( xvec, fxvec, linewidth = 3 )
    plt.plot ( [ x0, x0 ], [ 0, quartic ( x0 ) ], 'b', linewidth = 3 )
    plt.plot ( x0,  0.0, 'bo', markersize = 10 )
    plt.plot ( x0,  quartic ( x0 ), 'bo', markersize = 10 )
    plt.plot ( [ x, x ], [ 0, quartic ( x ) ], 'r', linewidth = 3 )
    plt.plot ( x,  0.0, 'ro', markersize = 10 )
    plt.plot ( x,  quartic ( x ), 'ro', markersize = 10 )
    plt.plot ( [ xmin, xmax ], [ 0, 0 ], 'k', linewidth = 3 )
    plt.grid ( True )
    plt.xlabel ( '<-- X -->' )
    plt.ylabel ( '<-- F(X) -->' )
    plt.title ( 'Initial point in blue, local minimizer in red', fontsize = 16 )
    filename = 'quartic_minimizer_x' + str(test) + '.png'
    plt.savefig ( filename )
    print ( '' )
    print ( '  Graphics saved as "%s"' % ( filename ) )
    plt.show ( )
    plt.close ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'quartic_gradient1:' )
  print ( '  Normal end of execution.' )
  return

def quartic ( x ):

#*****************************************************************************80
#
## quartic() evaluates the quartic function.
#
  f = 2.0 * x**4 - 4.0 * x**2 + x + 20.0

  return f

def quartic_df ( x ):

#*****************************************************************************80
#
## quartic_df() evaluates the derivative of the quartic function.
#
  df = 8.0 * x**3 - 8.0 * x + 1.0

  return df

if ( __name__ == '__main__' ):
  quartic_gradient1 ( )