# LogisticBifn.py:
#	Plot a bifurcation diagram for the Logistic map.

# My attempt to make a more storage efficient bifurcation diagram

# Import modules
from numpy import *
from pylab import *

# Define the Logistic map's function 
def LogisticMap(r,x):
	return r * x * (1.0 - x)

# Setup parameter range
rlow  = 3.0
rhigh = 4.0

# Setup the plot
# It's numbered 1 and is 8 x 6 inches.
figure(1,(8,6))
# Stuff parameter range into a string via the string formating commands.
TitleString = 'Logistic map, f(x) = r x (1 - x), bifurcation diagram for r in [%g,%g]' % (rlow,rhigh)
title(TitleString)
# Label axes
xlabel('Control parameter r')
ylabel('{X_n}')

# We tell plot()'s autoscaling the range in (r,x) we want by
#	Putting dots at the corners of the desired window.
#	Otherwise, the plot constantly rescales with new data.
plot([rlow], [0.0], 'k,')
plot([rlow], [1.0], 'k,')
plot([rhigh], [0.0], 'k,')
plot([rhigh], [1.0], 'k,')

# Set the initial condition used across the different parameters
ic = 0.2
# Establish the arrays to hold the set of iterates at each parameter value
# The iterates we'll throw away
nTransients = 200
# This sets how much the attractor is filled in
nIterates = 250
# This sets how dense the bifurcation diagram will be
nSteps = 400
# Sweep the control parameter over the desired range
rInc = (rhigh-rlow)/float(nSteps)
for r in arange(rlow,rhigh+rInc,rInc):
	# Set the initial condition to the reference value
	state = ic
	# Throw away the transient iterations
	for i in range(nTransients):
		state = LogisticMap(r,state)
	# Now store the next batch of iterates
	rsweep = [ ]   # The parameter value
	x = [ ]        # The iterates
	for i in range(nIterates):
		state = LogisticMap(r,state)
		rsweep.append(r)
		x.append( state )
	plot(rsweep, x, 'k,') # Plot the list of (r,x) pairs as pixels

# Use this to save figure as a bitmap png file
#savefig('LogisticCobweb', dpi=600)

# Display plot in window
show()