%matplotlib inline
import pylab as pl
from IPython.html.widgets import interact

from IPython.display import Image
Image(url="http://kristinruffin.files.wordpress.com/2014/03/jenga-falling-for-web.jpg")

def generate_brick_sequence(width, height, n):
    return [(width, height) for i in range(n)]

generate_brick_sequence(3, 4, 5)

from matplotlib.patches import Rectangle

def plot_brick_sequence(brick_sequence, shifts):
    ax = pl.gca()
    total_height = 0.
    for brick, shift in zip(brick_sequence, shifts):
        ax.add_patch(Rectangle((shift, total_height), brick[0], brick[1]))
        total_height += brick[1]

bricks = [(1, 2), (2, 3), (3, 4)]
shifts = [0, 0.5, 2]
plot_brick_sequence(bricks, shifts)
pl.xlim(0, 10)
pl.ylim(0, 10)

from IPython.display import SVG
SVG(filename='files/forces.svg')

SVG(filename='files/forces_equilibrium.svg')

def compute_horizontal_center_of_mass(brick_sequence, shifts):
    compute_mass = lambda b: b[0] * b[1]
    sum_xi_mi = sum([(shifts[i] + brick_sequence[i][0]/2.) * compute_mass(brick_sequence[i]) for i in range(len(brick_sequence))])
    sum_mi = sum([compute_mass(brick_sequence[i]) for i in range(len(brick_sequence))])
    return sum_xi_mi / sum_mi                                   

def compute_brick_shifts(brick_sequence):
    shifts = []
    reverse_brick_sequence = brick_sequence[::-1]
    for index, brick in enumerate(reverse_brick_sequence):
        if index == 0:
            shifts.append(0.)
        else:
            previous_center_of_mass = compute_horizontal_center_of_mass(reverse_brick_sequence[:index], shifts)
            shifts.append(previous_center_of_mass - brick[0])
    min_shift = min(shifts)
    return [s - min_shift for s in shifts[::-1]]

def plot_nbricks(l, n):
    bricks = generate_brick_sequence(l, 1, n)
    plot_brick_sequence(bricks, compute_brick_shifts(bricks))
    pl.xlim(0, n)
    pl.ylim(0, n)
    pl.grid(True)

interact(plot_nbricks,
            l=(0.1, 4, 0.1),
            n=(2, 20))

pl.figure(figsize=(10, 4))
format_plot = lambda : pl.xlim(0, 8) and pl.ylim(0, 40) and pl.grid(True);

pl.subplot(131)
plot_nbricks(2, 4)
format_plot()

pl.subplot(132)
plot_nbricks(2, 13)
format_plot()

pl.subplot(133)
plot_nbricks(2, 40)
format_plot()

l = 1
number_of_bricks = pl.linspace(1, 200, 50, pl.int32)
distances = pl.zeros_like(number_of_bricks)
for index, n in enumerate(number_of_bricks):
    bricks = generate_brick_sequence(l, 1, int(n))
    shifts = compute_brick_shifts(bricks)
    distances[index] = max(shifts) + l/2.

pl.plot(number_of_bricks, distances)
pl.title("distance reached as a function of number of {0} width bricks".format(l))
pl.xlabel('brick unit width (a. u.)')
pl.ylabel('distance covered (a. u.)')
pl.grid(True)

brick_widths = pl.linspace(0.1, 10)
distances = pl.zeros_like(brick_widths)
for l in brick_widths:
    bricks = generate_brick_sequence(l, 1, 100)
    distances[brick_widths == l] = max(compute_brick_shifts(bricks)) + l/2.

pl.plot(brick_widths, distances)
pl.xlabel('brick width (distance units)')
pl.ylabel('distance units')
pl.title('distance covered by stacking 100 bricks as a function of brick width')
pl.grid(True)

(distances[-1] - distances[0])/(10 - 0.1)

def generate_random_brick_sequence(n):
    return [(pl.rand(1), 0.2) for i in range(n)]

def plot_random_bricks(n):
    bricks = generate_random_brick_sequence(n)
    shifts =  compute_brick_shifts(bricks)
    plot_brick_sequence(bricks, shifts)

pl.figure(figsize=(11, 4))
format_plot = lambda : pl.xlim(0, 2.5) and pl.ylim(0, 9) and pl.grid(True);

pl.subplot(131)
plot_random_bricks(40)
format_plot()

pl.subplot(132)
plot_random_bricks(40)
format_plot()

pl.subplot(133)
plot_random_bricks(40)
format_plot()

SVG(filename='files/forces_analytical_formula.svg')

def analytical_shift(n):
    if n==1:
        return 0
    else:
        return 1./(n-1) * ((n-1) * analytical_shift(n-1) + 1/2.)

print [analytical_shift(i) for i in range(1, 10)]

bricks = generate_brick_sequence(1, 1, 10)
shifts = compute_brick_shifts(bricks)
shifts = shifts[::-1]
print [elem - shifts[0] for elem in shifts]