In [1]:

    

import numpy as np
import numpy.random as npr
import pylab as pl
%matplotlib inline


    

In [43]:

    

def simulate(length=100000):
    epsilon = 0.05 # variance of fluctuation
    alpha = 1.0
    mu = 1.0
    
    steps_per_unit_time=512
    tstep = 1.0 / steps_per_unit_time
    
    samples = 1
    
    locs = np.zeros((length,2))
    
    for i in range(samples):
        x=1.0
        y=0.0
        t=0
        
        while x>0.0 and t < length-1:
            t+=1
            grx,gry = np.sqrt(epsilon*tstep)*npr.randn(2)
            
            xa = x+tstep*x*(1.0-x**2 - alpha*(y**2))+grx
            ya = y-tstep*mu*(1.0+x**2)*y+gry
            
            xf = x+tstep*xa*(1.0-xa**2-alpha*ya**2)+grx
            yf = y-tstep*mu*(1.0+xa**2)*ya+gry
            
            x=(xf+xa)/2.0
            y=(yf+ya)/2.0
            
            locs[t] = (x,y)
            
    return locs


    

In [44]:

    

from numba import jit, autojit, double, float64, float32, void, int32


    

In [45]:

    

simulate_jit = jit(double[:,:](int32))(simulate)


    

In [47]:

    

%timeit simulate_jit(100000)


    

    




1 loops, best of 3: 2.15 s per loop

In [48]:

    

%timeit simulate(100000)


    

    




1 loops, best of 3: 8.2 s per loop

In [7]:

    

locs = np.array(locs)


    

In [17]:

    

thin = 10
pl.scatter(locs[::thin,0],locs[::thin,1],
           linewidths=0,alpha=0.5,
           cmap='Blues',
           c=range(len(locs[::thin])))
pl.colorbar()


    

    
Out[17]:





<matplotlib.colorbar.Colorbar instance at 0x1171330e0>






    

In [18]:

    

# path sampling
time = 16 # total time
tbins = time*steps_per_unit_time # total time steps


    

In [19]:

    

relax = tbins #MC moves per sample


    

In [22]:

    

x,y = np.zeros(tbins),np.zeros(tbins)
grx,gry = npr.randn(tbins),npr.randn(tbins)


    

In [23]:

    

transition=tbins


    

In [ ]:

    

def propagate(tbins,grx,gry):
    for t in range(tbins-1):
        xa = 
        
        xa = x+tstep*x*(1-x**2 - alpha*(y**2))+grx
        ya = y-tstep*mu*(1+x**2)*y+gry
        
        xf = x+tstep*xa*(1-xa**2-alpha*ya**2)+grx
        yf = y-tstep*mu*(1+xa**2)*ya+gry
        
        x=(xf+xa)/2
        y=(yf+ya)/2
    return transition


    

In [ ]:

    

#generate initial transition path
while transition==tbins:
    oldx=x[tbins-1]
    # generate trial trajectory
    t = npr.randint(tbins)
    oldgrx = grx[t]
    oldgry = gry[t]
    grx[t],gry[t] = np.sqrt(epsilon*tstep)*npr.randn(2)
    transition=propagate(tbins,tstep,)