<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 INF-495, v0.01, Claudio Torres, ctorres@inf.utfsm.cl. DI-UTFSM
Textbook
Step1:  First algorithm
Step2:  Second algorithm
Step3:  Third algorithm
Step4:  8.2 Stochastic predator-prey model
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import scipy.sparse.linalg as sp
import sympy as sym
from scipy.linalg import toeplitz
import ipywidgets as widgets
from ipywidgets import IntSlider
import matplotlib.pyplot as plt
%matplotlib inline
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
plt.style.use('ggplot')
import matplotlib as mpl
mpl.rcParams['font.size'] = 14
mpl.rcParams['axes.labelsize'] = 20
mpl.rcParams['xtick.labelsize'] = 14
mpl.rcParams['ytick.labelsize'] = 14
sym.init_printing()
from scipy.integrate import odeint
from ipywidgets import interact
def plot_matrices_with_values2(ax,M,flag_values=False):
    N=M.shape[0]
    cmap = plt.get_cmap('GnBu')
    ax.matshow(M, cmap=cmap)
    if flag_values:
        for i in np.arange(0, N):
            for j in np.arange(0, N):
                # I had to change the order of indices to make it right
                ax.text(i, j, '{:.2f}'.format(M[j,i]), va='center', ha='center', color='r')
ids = np.arange(64) # Making it Python compatible
I = np.reshape(ids,(8,8),order='F')
print(I)
# Defining occupied squares ([row,columns])
occupied = np.array([[2, 3], [5, 5], [6, 2], [6, 3]])-1 # Subtracting one to make it Python compatible!
print(occupied)
for oc in occupied:
    print(I[oc[0],oc[1]])
# Defining the peice
piece = 'Knight' # 'King' or 'Knight'
if piece == 'King':
    # Squares relative to the current one reached by the King
    Mask = np.zeros((8,2))
    k=0
    piece_range=[-1,0,1]
    for i in piece_range:
        for j in piece_range:
            if not ((i==0) and (j==0)):
                Mask[k,0]=i
                Mask[k,1]=j
                k=k+1
    print(Mask)
elif piece == 'Knight':
    # Squares relative to the current one reached by the Knight
    Mask = np.zeros((8,2))
    k=0
    piece_range=[-2,-1,0,1,2]
    for i in [-2,-1,0,1,2]:
        for j in [-2,-1,0,1,2]:
            if np.linalg.norm([i,j],ord=1)==3:
                Mask[k,0]=i
                Mask[k,1]=j
                k=k+1
    print(Mask)
else:
    print('Wrong piece!')
M = np.zeros((64,64))
# Scan through all starting squares
# find all possible destinations
for n in np.arange(8):
    for m in np.arange(8):
        k=I[n,m] # Index to current square
        # Find possible destination square using the mask
        P = np.array([n, m])+Mask
        # Discard occupied squares
        P_tmp=[]
        for p in P:
            flag=0
            for oc_sq in occupied:
                if abs(oc_sq[0]-p[0])+abs(oc_sq[1]-p[1])<1e-10:
                    flag=1
            if flag==0:
                P_tmp.append(p)
        P=np.array(P_tmp,dtype=int)
        # Discarding squares outside the board
        P_tmp=[]
        for p in P:
            if p[0]>=0 and p[0]<=7 and p[1]>=0 and p[1]<=7:
                P_tmp.append(p)
        P=np.array(P_tmp,dtype=int)
        nP = P.shape[0] # number of possible destinations
        # Fill in the probabilities
        for p in P:
            Ip = I[p[0],p[1]] # Destination index
            M[Ip,k] = 1/nP # Moving probability
fig,  (ax1)   = plt.subplots(1, 1, figsize=(10,10))
plot_matrices_with_values2(ax1,M,False)
Explanation: INF-495, v0.01, Claudio Torres, ctorres@inf.utfsm.cl. DI-UTFSM
Textbook: Computational Mathematical Modeling, An Integrated Approach Across Scales by Daniela Calvetti and Erkki Somersalo.
Chapter 8
8.1 Markov process and random walk
End of explanation
def plot_random_walk_output(nsteps=1):
    np.random.seed(0)
    # nsteps = 100000
    prev   = I[1,2]
    ncounts = np.zeros(64)
    squares_ids = [prev]
    for n in np.arange(nsteps):
        # Random draw using the column M[:,prev]
        next_step = np.random.choice(np.arange(64),p=M[:,prev])
        ncounts[next_step] = ncounts[next_step] + 1
        prev = next_step
        squares_ids.append(next_step)
    fig,  (ax1)   = plt.subplots(1, 1, figsize=(10,10))
    plot_matrices_with_values2(ax1,np.log10(np.reshape(ncounts/np.sum(ncounts),(8,8),order='F')),True)
    # print(ncounts)
    # print(squares_ids)
    # print(I)
interact(plot_random_walk_output,nsteps=(0,500000,1))
Explanation: First algorithm: Random walk. Now I added the plot of the resulting probabilities, it is interesting to observe the evolution from the beginning. Compare with the other two outcomes.
End of explanation
def plot_P_n(nj, I, M):
    # Defining initial distribution
    (i0,j0)=(0,1)
    pcurr = np.zeros(64)
    pcurr[I[i0,j0]] = 1
    for j in np.arange(nj):
        pcurr = np.dot(M,pcurr)
    pcurrM = np.reshape(pcurr,(8,8),order='F')
    fig,  (ax1)   = plt.subplots(1, 1, figsize=(10,10))
    plot_matrices_with_values2(ax1,np.log10(pcurrM),True)
    # print(pcurr)
interact(lambda nj=0: plot_P_n(nj, I=I, M=M), nj = (0,2000,1))
Explanation: Second algorithm: Power iteration. Warning: For the case of the Knight, this algorithm shows an interesting behavior and that seems to be related there are two eigenvalues with the magnitud equal to 1, one of the is actually equal to 1 but the other is equally to -1. This does not happen for the case of the King.
End of explanation
w, v = np.linalg.eig(M)
w
v0=np.real(v[:,0])/np.sum(np.real(v[:,0]))
print(v0)
v1=np.real(v[:,1])/np.sum(np.real(v[:,1]))
print(v1)
if piece == 'King':
    v_selected = v0
elif piece == 'Knight':
    v_selected = v1
else:
    print('Wrong piece!')
    v_selected=np.full((64),np.nan)
fig,  (ax1)   = plt.subplots(1, 1, figsize=(10,10))
plot_matrices_with_values2(ax1,np.log10(np.reshape(np.real(v_selected),(8,8),order='F')),True)
Explanation: Third algorithm: Eigenvalues and eigenvectors. We need to look for the eigenvector associated to the eigenvalue equal to 1.
End of explanation
def plot_S_Pred_Prey(seed=0):
    Nj       = 200 # number of X adult pair at t = t(j)
    Mj       = 50  # number of Y adult pair at t = t(j)
    p_mean   = 10  # mean number of X offsprings/adult pair
    p_std    = 3   # std of the number of the X offspring/adult pair
    beta     = 0.3 # Survival probability of X in the absence of Y
    Mstar    = 50  # number of Y that halves the survival prob. of X
    q_mean   = 4   # mean number of Y offspring/adult pair
    q_std    = 2   # std of the number of the Y offsrping/adult pair
    gamma    = 0.6 # Survival prob. of Y in adundance of X
    Nstar    = 400 # Number of X pairs at which the survival prob. of Y is halved
    nrounds  = 40  # Number of generations allowed
    N = np.full(nrounds,np.nan)
    N[0] = Nj
    M = np.full(nrounds,np.nan)
    M[0] = Mj
    Surv = np.full((nrounds-1,2),np.nan)
    np.random.seed(seed)
    for j in np.arange(1,nrounds):
        # Draw the number of offsprings of species X
        Pj = Nj*p_mean+np.sqrt(Nj)*p_std*np.random.randn()
        # Survival probability of the offspring of species X
        s = beta/(1+(Mj/Mstar)**2)
        Nj = s*Pj+np.sqrt(s*Pj*(1-s))*np.random.randn()
        N[j] = 0.5*Nj
        # Draw the number of offsprings of species Y
        Qj = Mj*q_mean+np.sqrt(Mj)*q_std*np.random.randn()
        # Survival probability of the offspring of species Y
        r = gamma*(1-1/(1+(Nj/Nstar)**2))
        Mj = r*Qj+np.sqrt(r*Qj*(1-r))*np.random.randn()
        M[j] = 0.5*Mj
        # Saving the survival probabilities for later analysis
        Surv[j-1,:]=[s,r]
    s=Surv[:,0]
    r=Surv[:,1]
    
    smean = np.sum(s)/(nrounds-1)
    rmean = np.sum(r)/(nrounds-1)
    Surv_c = Surv-np.array([smean,rmean])
    covSurv = np.dot(Surv_c.T,Surv_c)/(nrounds-1)
    corrcoeff = covSurv[0,1]/np.sqrt(covSurv[0,0]*covSurv[1,1])
        
    fig,  ((ax1,ax2),(ax3,ax4))   = plt.subplots(2, 2, figsize=(10,10))
    
    ax4.plot(N,'r.-')
    ax4.plot(M,'b.-')
    ax4.grid(True)
    
    s=Surv[:,0]
    r=Surv[:,1]
    ax1.plot(s,r,'.')
    ax1.set_xlim(0,1)
    ax1.set_ylim(0,1)
    ax1.set_xlabel(r'$s_i$')
    ax1.set_ylabel(r'$r_i$')
    ax1.text(0.5, 0.2, 'c={:.5f}'.format(corrcoeff), va='center', ha='center', color='k')
    
    ax2.hist(r,orientation='horizontal')
    ax2.set_ylim(0,1)
    
    ax3.hist(s)
    ax3.set_xlim(0,1)
    
    plt.show()
interact(plot_S_Pred_Prey,seed=(0,100,1))
Explanation: 8.2 Stochastic predator-prey model
End of explanation
<END_TASK>