%pylab inline
import numpy as np
import scipy.optimize
import scipy.sparse
import matplotlib.pyplot as plt
from IPython.display import clear_output
import time

delta=0.0045; tau1=0.02; tau2=0.2; alpha=0.899; beta=-0.91; gamma=-alpha;

def f(u,v):
    return alpha*u*(1-tau1*v**2) + v*(1-tau2*u);

def g(u,v):
    return beta*v*(1+alpha*tau1/beta*u*v) + u*(gamma+tau2*v);

def five_pt_laplacian_sparse_periodic(m,a,b):
    """Construct a sparse matrix that applies the 5-point laplacian discretization
       with periodic BCs on all sides."""
    e=np.ones(m**2)
    e2=([1]*(m-1)+[0])*m
    e3=([0]+[1]*(m-1))*m
    h=(b-a)/(m+1)
    A=scipy.sparse.spdiags([-4*e,e2,e3,e,e],[0,-1,1,-m,m],m**2,m**2)
    # Top & bottom BCs:
    A_periodic = scipy.sparse.spdiags([e,e],[m-m**2,m**2-m],m**2,m**2).tolil()
    # Left & right BCs:
    for i in range(m):
        A_periodic[i*m,(i+1)*m-1] = 1.
        A_periodic[(i+1)*m-1,i*m] = 1.
    A = A + A_periodic
    A/=h**2
    A = A.todia()
    return A

A = five_pt_laplacian_sparse_periodic(4,-1.,1.)
plt.spy(A)

def one_step(u,v,k,A,delta,D1=0.5,D2=1.0):
    u_new = u + k * (delta*D1*A*u + f(u,v))
    v_new = v + k * (delta*D2*A*v + g(u,v))
    
    return u_new, v_new

delta=0.0021; tau1=3.5; tau2=0; alpha=0.899; beta=-0.91; gamma=-alpha;

def step_size(h,delta):
    return h**2/(5.*delta)

def pattern_formation(m=10,T=1000):
    r"""Model pattern formation by solving a reaction-diffusion PDE on a periodic
        square domain with an m x m grid."""
    D1 = 0.5
    D2 = 1.0
    
    # Set up the grid
    a=-1.; b=1.
    h=(b-a)/m;                 # Grid spacing
    x = np.linspace(a,b,m)     # Coordinates
    y = np.linspace(a,b,m)
    Y,X = np.meshgrid(y,x)

    # Initial data
    u=np.random.randn(m,m)/2.;
    v=np.random.randn(m,m)/2.;

    #plt.clf(); plt.hold(False)
    #plt.pcolormesh(x,y,u); plt.colorbar(); plt.axis('image');
    #plt.draw()

    frames = [u]
    
    u=u.reshape(-1)
    v=v.reshape(-1)

    A=five_pt_laplacian_sparse_periodic(m,-1.,1.)

    t=0.                      # Initial time
    k = step_size(h,delta)    # Time step size
    N = int(round(T/k))      # Number of steps to take

    
    #Now step forward in time
    next_plot = 0
    for j in range(N):

        u,v = one_step(u,v,k,A,delta)
        
        t = t+k;

        #Plot every t=5 units
        if t>next_plot:
            next_plot = next_plot + 5
            U=u.reshape((m,m))
            #plt.clf()
            #plt.pcolormesh(x,y,U)
            #plt.colorbar()
            #plt.axis('image')
            #plt.title(str(t))
            #time.sleep(0.2)
            #clear_output()
            #fi=plt.gcf()
            #display(fi)  
            #plt.clf()
            frames.append(U)
            
    return X,Y,frames

from matplotlib import animation
import matplotlib.pyplot as plt
from clawpack.visclaw.JSAnimation import IPython_display
import numpy as np

fig = plt.figure(figsize=[4,4])

U = frames[0]

# This essentially does a pcolor plot, but it returns the appropriate object
# for use in animation.  See http://matplotlib.org/examples/pylab_examples/pcolor_demo.html.
# Note that it's necessary to transpose the data array because of the way imshow works.
im = plt.imshow(U.T, vmin=U.min(), vmax=U.max(),
           extent=[x.min(), x.max(), y.min(), y.max()],
           interpolation='nearest', origin='lower')

def fplot(frame_number):
    U = frames[frame_number]
    im.set_data(U.T)
    return im,

x,y,frames = pattern_formation(m=200,T=300)

animation.FuncAnimation(fig, fplot, frames=len(frames), interval=20)

def one_step(u,v,k,A,delta,D1=0.5,D2=1.0):
    # Your code here
    return u_new, v_new

def step_size(h,delta):
    return h/(10.*delta)

pattern_formation(m=150,T=300)

epsilon = 0.05

# Grid
m = 64
x = np.arange(-m/2,m/2)*(2*np.pi/m)
k = 1./m**2
tmax = 30.

# Initial data
u = np.sin(x)**2 * (x<0.)
v = np.fft.fft(u)

xi = np.array([range(0,m/2) + [0] + range(-m/2+1,0)])
eps_xi2 = epsilon * xi**2.
E = np.exp(-k * epsilon * xi**2.)

nplt = np.floor((tmax/25)/k)
nmax = int(round(tmax/k))

for n in range(1,nmax+1):
    v = E*v
    t = n*k
    
    # Plotting
    if np.mod(n,nplt) == 0:
        u = np.squeeze(np.real(np.fft.ifft(v)))
        plt.clf()
        plt.plot(x,u,linewidth=3)
        plt.title('t='+str(t))
        plt.xlim((-np.pi,np.pi))
        plt.ylim((0.,1.))
        time.sleep(0.2)
        clear_output()
        fi=plt.gcf()
        display(fi)  
        plt.clf()



delta=0.0045; tau1=2.02; tau2=0.; alpha=2.0; beta=-0.91; gamma=-alpha;
delta=0.0005; tau1=2.02; tau2=0.; alpha=2.0; beta=-0.91; gamma=-alpha;
delta=0.0021; tau1=3.5; tau2=0; alpha=0.899; beta=-0.91; gamma=-alpha;
delta=0.0045; tau1=0.02; tau2=0.2; alpha=1.9; beta=-0.85; gamma=-alpha;
delta=0.0001; tau1=0.02; tau2=0.2; alpha=0.899; beta=-0.91; gamma=-alpha;
delta=0.0045; tau1=0.02; tau2=0.2; alpha=1.9; beta=-0.91; gamma=-alpha;