#%qtconsole

from IPython.core.display import HTML
def css_styling():
    styles = open("styles/custom.css", "r").read()
    return HTML(styles)
css_styling()

import json
s = json.load( open("styles/bmh_matplotlibrc.json") )
matplotlib.rcParams.update(s)

from __future__ import  division
from matplotlib.lines import Line2D
from matplotlib.patches import Rectangle

def needle_gen(L=0.5):
    '''Drops needle on unit square. Return dictionary describing needle position and 
    orientation. The `cut` key in the dictionary is True when the needle intersects 
    an edge of the unit square and is False otherwise.
    
    L := length of needle (default = 0.5)
    '''
    assert 0<L<1
    xc,yc = rand(2) # uniform random in [0,1] for each coordinate dimension
    ang = rand()*pi # uniform random angle in [0,pi]
    x0,y0 = xc-L/2*cos(ang),yc-L/2*sin(ang) # coordinates for one end of needle 
    xe,ye = xc+L/2*cos(ang),yc+L/2*sin(ang) # coordinates for the other end
    not_cut = (0<x0<1) and (0<y0<1) and (0<xe<1) and (0<ye<1)
    return dict(x0=x0,y0=y0,xc=xc,yc=yc,xe=xe,ye=ye,ang=ang,cut=( not not_cut ))


def draw_needle(ax,d):
    '''
    Draw needle symbol on given axis.
    
    If red, then the line cuts an edge of the unit square. If blue, then it does not. 
    One end of the needle is marked with a `o` symbol and the other end is unmarked. 
    
    ax := matplotlib axes to draw needle symbol
    d := dictionary output from needle_gen
    '''
    for i,j in d.iteritems(): exec('%s=%r'%(i,j)) # dump dictionary in local namespace
    if cut: 
        line= [Line2D([x0,xe],[y0,ye],color='r',alpha=0.3), 
               Line2D([xc],[yc],marker='o',color='r')] 
    else:   
        line= [Line2D([x0,xe],[y0,ye],alpha=0.3), 
               Line2D([xc],[yc],marker='o',alpha=0.3)] 
    for i in line: 
        ax.add_line(i)

L=0.5

samples = [needle_gen() for i in range(500)]

def draw_sim(samples,L=0.5):
    'Draw simulation results. Package this plot for reuse later'
    fig,ax=subplots()
    fig.set_size_inches(5,5)
    
    red_count=0
    n=len(samples)
    for k in samples:
        red_count+=k['cut']
        draw_needle(ax,k)
        
    ax.set_aspect(1) # set aspect ratio to 1
    ax.add_patch(Rectangle((0,0),1,1,alpha=0.3,color='y')) # unit square background is light yellow
    ax.set_xlabel('x-direction')
    ax.set_ylabel('y-direction')
    ax.axis([-.1,1.1,-.1,1.1]) # add some space around the unit square
    ax.set_title(r'$\mathbb{P}$ (cut)=%3.2f'% (red_count/n))
    ax.grid()

draw_sim(samples)    

fig,axs = subplots(1,2)
fig.set_size_inches((12,2))
ax=axs[0]
samples=np.array([[needle_gen(L)['cut'] for i in range(300)] for k in range(100)])

ax.plot(np.mean(samples,axis=1),'o-',alpha=0.3)
ax.axis(xmax=samples.shape[0])
ax.set_xlabel(r'$trial$',fontsize=18)
ax.set_ylabel('$\mathbb{\hat{P}}(cut)$',fontsize=18)
ax.set_title('$\overline{\mathbb{\hat{P}}}(cut) = %3.3f,\hat{\sigma}=%3.3f$'%
             (samples.mean(),(samples.mean(1)).std()))
ax.grid()

ax = axs[1]
ax.hist( samples.mean(1));
ax.set_title('Distribution of Mean Estimates');

fig,ax=subplots()
fig.set_size_inches(12,2)
ax.plot(samples.mean(1))
ax.set_title(r'$\hat{\sigma}=%2.3f,R_e=%2.3f$'% (samples.mean(1).std(),samples.mean(1).std()/samples.mean()), fontsize=18)
ax.set_xlabel('sample history index',fontsize=16)
ax.set_ylabel('$\mathbb{\hat{P}}(cut)$',fontsize=18);

L=0.1

samples_d=[needle_gen(L) for i in range(300*100)]
samples = np.array([k['cut'] for k in samples_d]).reshape(100,300)

def draw_sample_history(samples):
    fig,ax=subplots()
    fig.set_size_inches(12,2)
    mn = samples.mean(1)
    ax.plot(mn)
    ax.set_title(r'$\hat{\sigma}=%2.3f,R_e=%2.3f$'% (mn.std(),
                                                     mn.std()/mn.mean()), fontsize=18)
    ax.set_xlabel('sample history index',fontsize=16)
    ax.set_ylabel('$\mathbb{\hat{P}}(cut)$',fontsize=18);

draw_sample_history(samples)    

draw_sim(samples_d[:500])   

def needle_gen_b(L=0.5):
    '''Drops needle on unit square using a biased strategy. Return dictionary describing 
    needle position and orientation. The `cut` key in the dictionary is True when the
    needle intersects an edge of the unit square and is False otherwise.
    
    L := length of needle (default = 0.5)
    '''
    assert 0<L<1
    yc = rand() # uniform on unit interval
    r = L*rand()/2.
    xc = r if rand()<0.5 else (1-r) # pick either edge 50/50
    ang = rand()*pi # uniform random angle in [0,pi]
    x0,y0 = xc-L/2*cos(ang),yc-L/2*sin(ang) # coordinates for one end of needle 
    xe,ye = xc+L/2*cos(ang),yc+L/2*sin(ang) # coordinates for the other end
    not_cut = (0<x0<1) and (0<y0<1) and (0<xe<1) and (0<ye<1)
    return dict(x0=x0,y0=y0,xc=xc,yc=yc,xe=xe,ye=ye,ang=ang,cut=( not not_cut ))

samples_b=[needle_gen_b(L) for i in range(300*100) ] # create biased samples
draw_sim(samples_b[:500])

no_bias_samples = samples # save for later
samples = np.array([k['cut'] for k in samples_b]).reshape(100,300)

draw_sample_history(samples);

fig,ax=subplots()
fig.set_size_inches(10,3)
ax.hist(no_bias_samples.mean(1),alpha=0.3,normed=1,label='unbiased')
ax.hist(samples.mean(1)*2*L,alpha=0.3,normed=1,label='biased')
ax.legend(loc=0)
ax.set_title('Biased samples more efficiently produce useful samples')

prob_cut= lambda l:(1-(1-2*l/pi)**2)