# HenonMap.py:
#   Estimate the Henon map's Lyapunov Characteristic Exponent.
#   Plot out iterates of Henon's 2D map, listing the LCE estimate.
#
# The Henon Map is given by
#     x_n+1 = f( x_n , y_n )
#     y_n+1 = g( x_n , y_n )
# with
#     f(x,y) = y + 1 - a x^2
#     g(x,y) = b x
#
# The state space is R^2
# and the control parameters range in
#       a in [0,2]
#       b in [0,1]
# Henon's original parameters for a chaotic attractor:
#       (a,b) = (1.4,0.3)
#
# It has Jacobian (d = partial derivative)
#     df/dx (x,y) = - 2 a x
#     df/dy (x,y) = 1
#     dg/dx (x,y) = b
#     dg/dy (x,y) = 0
# and area contraction given by
#     Determinant(Jacobian(x,y)) = -b
# 
# We look at the symmetric part of J: S = (J + J*)/2, * = transpose
#          [ -2 a x       (1+b)/2 ]
#     S =  [                      ]
#          [ (1+b)/2         0    ]
# and get the two eigenvalues explicitly
#     Eigen1 = - a x - sqrt(a^2 x^2 + (1 + b)^2 )
#     Eigen2 = - a x + sqrt(a^2 x^2 + (1 + b)^2 )
# Largest LCE = average of log(abs(max{Eigen1,Eigen2}(x_n,y_n))
#    along a trajectory
# Smallest LCE = average of log(abs(min{Eigen1,Eigen2}(x_n,y_n))
#    along a trajectory

# Note: Return is a list: the next (x,y)
def HenonMap(a,b,x,y):
	return y + 1.0 - a *x*x, b * x

# Return is a 2 x 2 array
def HenonMapJacobian(a,b,x,y):
	return array( [ [-2.0*a*x,1.0] , [b,0.0] ] )

# The area contraction rate
def HenonMapDetJac(a,b,x,y):
	return -b

# The maximum eigenvalue of the symmetric part of the Jacobian
#  These are all that is used for the LCE below
#  We could estimate the eigenvalue numerically, but in 2D we can
#     work them out directly.
def Eigen1(a,b,x,y):
	return - a * x - sqrt( a*a*x*x + (1.0 + b) * (1.0 + b) )

# The minimum eigenvalue of the symmetric part of the Jacobian
def Eigen2(a,b,x,y):
	return - a * x + sqrt( a*a*x*x + (1.0 + b) * (1.0 + b) )

# Import plotting routines
# Use this on Linux/Unix/OS X
from pylab import *
# Use this on Windows
#from matplotlib.matlab import *

# Simulation parameters
#
# Control parameters:
a = 1.4
b = 0.3
# The number of iterations to throw away
nTransients = 100
# The number of iterations to generate
nIterates = 1000

# Initial condition for iterates to plot
xtemp = 0.1
ytemp = 0.3
for n in range(0,nTransients):
	xtemp, ytemp = HenonMap(a,b,xtemp,ytemp)
# Set up arrays of iterates (x_n,y_n) and set the initital condition
x = [xtemp]
y = [ytemp]
# The main loop that generates iterates and stores them
for n in range(0,nIterates):
	# at each iteration calculate (x_n+1,y_n+1)
	xtemp, ytemp = HenonMap(a,b,x[n],y[n])
	# and append to lists x and y
	x.append( xtemp )
	y.append( ytemp )

# Estimate the LCEs
# The number of iterations to throw away
nTransients = 200
# The number of iterations to generate
nIterates = 100000
# Initial condition
xState = 0.1
yState = 0.3
# An array to hold the Jacobian
for n in range(0,nTransients):
	xState, yState = HenonMap(a,b,xState,yState)
maxLCE = 0.0
minLCE = 0.0
# Begin estimation
for n in range(0,nIterates):
  # Get next state
	xState, yState = HenonMap(a,b,xState,yState)
  # Evaluate stability there
	e1 = abs(Eigen1(a,b,xState,yState))
	e2 = abs(Eigen2(a,b,xState,yState))
	if e1 > e2:
		maxLCE = maxLCE + log(e1)
		minLCE = minLCE + log(e2)
	else:
		maxLCE = maxLCE + log(e2)
		minLCE = minLCE + log(e1)
#	print e1,e2

# Convert to per-iterate LCEs and to log_2
maxLCE = maxLCE / float(nIterates) / log(2.)
minLCE = minLCE / float(nIterates) / log(2.)
print maxLCE, minLCE
# Setup the plot
xlabel('x_n') # set x-axis label
ylabel('y_n') # set y-axis label
# Set plot title
LCEString = '(%g,%g)' % (maxLCE,minLCE)
PString   = '(%g,%g)' % (a,b)
title('Henon LCEs = ' + LCEString + ' at (a,b) = ' + PString)
# Plot the time series
plot(x,y, 'r,')
# Try to scale the plotting box (This doesn't give the desired control!)
axis([-2.0, 2.0, -1.0, 1.0])

# Use command below to save figure
#savefig('HenonMapIterates', dpi=600)

# Display the plot in a window
show()