MTH 337: Week 12

%pylab

Using matplotlib backend: Qt4Agg
Populating the interactive namespace from numpy and matplotlib

The "Stumble Down" Algorithm

def lake_depth(x,y):
    return (x-5)**2 + (y-3)**2 + 17

def calculate_stumbledown(myf, p, stepsize, nsteps):
        stumble = empty((nsteps+1, 2))
        successes = empty(nsteps+1, dtype=bool)
        stepsizes = empty(nsteps+1)
        successes[0] = True
        stumble[0] = p
        stepsizes[0] = stepsize
        counter = 0
        for i in xrange(1, nsteps+1):
            step = stepsize*(random.rand(2)-0.5)
            q = p + step
            if myf(*q) < myf(*p):
                stumble[i] = q
                successes[i] = True
                p = q
                counter = 0
            else:
                stumble[i] = q
                successes[i] = False
                counter += 1
                if counter >=10:
                    stepsize = float(stepsize)*(1/2.)
                    counter = 0
                    if stepsize < 10e-8:
                        print i
                        break
            stepsizes[i] = stepsize
        return stumble, successes, stepsizes

Animating the "Stumble Down" Algorithm

from matplotlib import animation

steps = 120
h = 1.5
p = (0., 0.)
fig, ax = subplots(figsize=(6,6))
ax.set_aspect('equal')
point, = plot([],[],'go', ms=3)
line, = plot([],[],'k')
line2, = plot([],[],'r')
rect = gca().add_patch(Rectangle((0,0), 0,0, fill=False, ec='y'))
stumble, successes, stepsizes = calculate_stumbledown(lake_depth, p, h, steps)
max_steps = ones(2)*stepsizes[0]
maxx, maxy = amax(stumble,axis=0) + max_steps
minx, miny = amin(stumble,axis=0) - max_steps
xlim(minx,maxx)
ylim(miny,maxy)

def init(): 
    point.set_data([],[])
    x = linspace(minx,maxx,100)
    y = linspace(miny,maxy,100)
    X,Y = meshgrid(x,y)
    Z = lake_depth(X,Y)
    contourf(x,y,Z,cmap='gist_rainbow', levels = arange(100))
    return line, line2, point, rect,

def animate(i):
    x,y = stumble[i] 
    path = stumble[:i+1][successes[:i+1]]
    px, py = path[-1]
    point.set_data(px,py) 
    line.set_data(path[:,0], path[:,1])
    line2.set_data([px,x], [py,y])
    stepsize = stepsizes[i]
    rect.set_xy((px-stepsize, py-stepsize))
    rect.set_height(2*stepsize)
    rect.set_width(2*stepsize)
    return line, line2, point, rect,

animation.FuncAnimation(fig, animate, init_func = init, frames = steps, interval = 250, blit = False, repeat = True)

<matplotlib.animation.FuncAnimation at 0x7f7a5ce94ed0>