"""Basic functionality for iterated maps"""

import scipy
import pylab

# ***** Shared software, 
# ***** used by ChaosLyapunov, InvariantMeasure, FractalDimension, and 
# ***** PeriodDoubling exercises.

def f(x, mu):
    """
    Logistic map f(x) = 4 mu x (1-x), which folds the unit interval (0,1)
    into itself.
    """
    return 4*mu*x*(1.-x)


def Iterate(g, x0, N, mu):
    """
    Iterate the function g(x,mu) N times, starting at x=x0.
    Return g(g(...(g(x))...)). Used to find a point on the attractor
    starting from some arbitrary point x0.

    Calling Iterate for the Feigenbaum map at mu=0.9 would look like
        Iterate(f, 0.1, 1000, 0.9)

    We'll later be using Iterate to study the sine map 
        fsin(x,B) = B sin(pi x)
    so passing in the function and arguments will be necessary for 
    comparing the logistic map f to fsin.
    
    Inside Iterate you'll want to apply g(x0, mu).
    """
    for i in range(N):
        x0 = g(x0, mu)
    return x0


def IterateList(g, x0, N, mu):
    """
    Iterate the function g(x, mu) N-1 times, starting at x0, so that the
    full trajectory contains N points.
    Returns the entire list 
    (x, g(x), g(g(x)), ... g(g(...(g(x))...))). 

    Can be used to explore the dynamics starting from an arbitrary point 
    x0, or to explore the attractor starting from a point x0 on the 
    attractor (say, initialized using Iterate).

    For example, you can use Iterate to find a point xAttractor on the 
    attractor and IterateList to create a long series of points attractorXs
    (thousands, or even millions long, if you're in the chaotic region), 
    and then use
        pylab.hist(attractorXs, bins=500, normed=1)
        pylab.show()
    to see the density of points.
    """
    xs = [x0]
    for i in range(N-1):
        x0 = g(x0, mu)
        xs.append(x0)
    return xs


def BifurcationDiagram(g, x0, nTransient, nCycle, muArray, showPlot=True):
    """
    For each parameter value mu in muArray,
    iterate g nTransient times to find a point on the attractor, and then
    make a list nCycle long to explore the attractor.

    To generate muArray, it's convenient to use scipy.linspace: for example,
    BifurcationDiagram(f, 0.1, 500, 128, scipy.linspace(0.8, 1.0, 200))

    pylab.plot allows one to plot an entire array of abscissa-values versus an
    array of ordinate-values of the same shape. Our vertical axis (ordinate)
    is a list of arrays of attractor points of length nCycle (created by 
    IterateList after Iterating), one list per value of mu in muArray.
    Our horizontal axis (abscissa) should thus be a list of arrays 
        [mu, mu, mu, ...] = [mu]*nCycle = mu*scipy.ones(nCycle) 
    of length nCycle.
    Use
        pylab.plot(muMatrix, xMatrix, 'k,')
        pylab.show()
    to visualize the resulting bifurcation diagram, where 'k,' denotes 
    black pixels.
    """
    muMatrix = []
    xMatrix = []
    for mu in muArray:
        muMatrix.append([mu]*nCycle)
        xMatrix.append(IterateList(g, Iterate(g, x0, nTransient, mu),
                              nCycle, mu))
    pylab.plot(muMatrix, xMatrix, 'k,')
    if showPlot: pylab.show()

def demo():
    """Demonstrates solution for exercise: example of usage"""
    print("IterateLogistic Demo")
    print("  Creating Bifurcation Diagram")
    BifurcationDiagram(f, 0.1, 500, 500, scipy.linspace(0., 1.00001, 200))
    print("  Creating Bifurcation Diagram (zoom)")
    BifurcationDiagram(f, 0.1, 500, 500, scipy.linspace(0.85, 1.00001, 300))
    print("  Creating Bifurcation Diagram (zoom again)")
    BifurcationDiagram(f, 0.1, 500, 500, scipy.linspace(0.96, 0.97, 400))

if __name__=="__main__":
   demo()

# Copyright (C) Cornell University
# All rights reserved.
# Apache License, Version 2.0