import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

%matplotlib inline

def df_dt(x, t, a, b, c, d):
    
    dx = a * x[0] - b * x[0] * x[1]
    dy = - c * x[1] + d * x[0] * x[1]
    
    return np.array([dx, dy])

# Parámetros
a = 0.1
b = 0.02
c = 0.3
d = 0.01

# Condiciones iniciales
x0 = 40
y0 = 9
conds_iniciales = np.array([x0, y0])

# Condiciones para integración
tf = 200
N = 800
t = np.linspace(0, tf, N)

solucion = odeint(df_dt, conds_iniciales, t, args=(a, b, c, d))

plt.style.use('pybonacci')

plt.figure("Evolución temporal", figsize=(8,5))
plt.title("Evolución temporal")
plt.plot(t, solucion[:, 0], label='presa')
plt.plot(t, solucion[:, 1], label='depredador')
plt.xlabel('tiempo')
plt.ylabel('población')
plt.legend()
plt.savefig('evolucion_temporal.png')

plt.figure("Presas vs depredadores", figsize=(8,5))
plt.plot(solucion[:, 0], solucion[:, 1])
plt.xlabel('presas')
plt.ylabel('depredadores')
plt.savefig('presas_vs_depredadores.png')

x_max = np.max(solucion[:,0]) * 1.05
y_max = np.max(solucion[:,1]) * 1.05

x = np.linspace(0, x_max, 25)
y = np.linspace(0, y_max, 25)

xx, yy = np.meshgrid(x, y)
uu, vv = df_dt((xx, yy), 0, a, b, c, d)
norm = np.sqrt(uu**2 + vv**2)
uu = uu / norm
vv = vv / norm

plt.figure("Campo de direcciones", figsize=(8,5))
plt.quiver(xx, yy, uu, vv, norm, cmap=plt.cm.gray)
plt.plot(solucion[:, 0], solucion[:, 1])
plt.xlim(0, x_max)
plt.ylim(0, y_max)
plt.xlabel('presas')
plt.ylabel('depredadores')
plt.savefig('campo_direcciones.png')

n_max = np.max(solucion) * 1.10

fig, ax = plt.subplots(1,2)

fig.set_size_inches(12,5)

ax[0].quiver(xx, yy, uu, vv, norm, cmap=plt.cm.gray)
ax[0].plot(solucion[:, 0], solucion[:, 1], lw=2, alpha=0.8)
ax[0].set_xlim(0, x_max)
ax[0].set_ylim(0, y_max)
ax[0].set_xlabel('presas')
ax[0].set_ylabel('depredadores')

ax[1].plot(t, solucion[:, 0], label='presa')
ax[1].plot(t, solucion[:, 1], label='depredador')
ax[1].legend()
ax[1].set_xlabel('tiempo')
ax[1].set_ylabel('población')
plt.savefig('campo_direcciones_ev_temporal.png')

def C(x, y, a, b, c, d):
    return a * np.log(y) - b * y + c * np.log(x) - d * x

x = np.linspace(0, x_max, 100)
y = np.linspace(0, y_max, 100)
xx, yy = np.meshgrid(x, y)
constant = C(xx, yy, a, b, c, d)

plt.figure('distintas_soluciones', figsize=(8,5))
plt.contour(xx, yy, constant, 50, cmap=plt.cm.Blues)
plt.xlabel('presas')
plt.ylabel('depredadores')
plt.savefig('distintas_soluciones.png')

#n_max = np.max(solucion) * 1.10

fig, ax = plt.subplots(1,2)

fig.set_size_inches(12,5)

ax[0].plot(solucion[:, 0], solucion[:, 1], lw=2, alpha=0.8)
ax[0].scatter(c/d, a/b)
levels = (0.5, 0.6, 0.7, 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.775, 0.78, 0.781)
ax[0].contour(xx, yy, constant, levels, colors='blue', alpha=0.3)
ax[0].set_xlim(0, x_max)
ax[0].set_ylim(0, y_max)
ax[0].set_xlabel('presas')
ax[0].set_ylabel('depredadores')

ax[1].plot(t, solucion[:, 0], label='presa')
ax[1].plot(t, solucion[:, 1], label='depredador')
ax[1].legend()
ax[1].set_xlabel('tiempo')
ax[1].set_ylabel('población')
plt.savefig('distintas_soluciones_ev_temporal.png')

def logistic_curve(t, a=1, m=0, n=1, tau=1):
    e = np.exp(-t / tau)
    return a * (1 + m * e) / (1 + n * e) 

x_ = np.linspace(0,10)
plt.figure('función logística', figsize=(8,5))
plt.plot(x_, logistic_curve(x_, 1, m=10, n=100, tau=1))
plt.savefig('funcion_logistica.png')

def df_dt_logistic(x, t, a, b, c, d, r):
    
    dx = a * x[0] - r * x[0]**2 - b * x[0] * x[1]
    dy = - c * x[1] + d * x[0] * x[1]
    
    return np.array([dx, dy])

# Parámetros
a = 0.1
b = 0.02
c = 0.3
d = 0.01
r = 0.001

# Condiciones iniciales
x0 = 40
y0 = 9
conds_iniciales = np.array([x0, y0])

# Condiciones para integración
tf = 200
N = 800
t = np.linspace(0, tf, N)

solucion_logistic = odeint(df_dt_logistic, conds_iniciales, t, args=(a, b, c, d, r))

n_max = np.max(solucion) * 1.10

fig, ax = plt.subplots(1,2)

fig.set_size_inches(12,5)

x_max = np.max(solucion_logistic[:,0]) * 1.05
y_max = np.max(solucion_logistic[:,1]) * 1.05

x = np.linspace(0, x_max, 25)
y = np.linspace(0, y_max, 25)

xx, yy = np.meshgrid(x, y)
uu, vv = df_dt_logistic((xx, yy), 0, a, b, c, d, r)
norm = np.sqrt(uu**2 + vv**2)
uu = uu / norm
vv = vv / norm

ax[0].quiver(xx, yy, uu, vv, norm, cmap=plt.cm.gray)
ax[0].plot(solucion_logistic[:, 0], solucion_logistic[:, 1], lw=2, alpha=0.8)
ax[0].set_xlim(0, x_max)
ax[0].set_ylim(0, y_max)
ax[0].set_xlabel('presas')
ax[0].set_ylabel('depredadores')

ax[1].plot(t, solucion_logistic[:, 0], label='presa')
ax[1].plot(t, solucion_logistic[:, 1], label='depredador')
ax[1].legend()
ax[1].set_xlabel('tiempo')
ax[1].set_ylabel('población')
plt.savefig('campo_direcciones_ev_temporal_caso2.png')

from IPython.html.widgets import interact

def solucion_temporal_interact(a, b, c, d, x0, y0, tf):
    
    conds_iniciales = np.array([x0, y0])

    # Condiciones para integración
    N = 800
    t = np.linspace(0, tf, N)
    
    solucion = odeint(df_dt, conds_iniciales, t, args=(a, b, c, d))
    
    plt.figure("Evolución temporal", figsize=(8,5))
    plt.title("Evolución temporal")
    plt.plot(t, solucion[:, 0], label='presa')
    plt.plot(t, solucion[:, 1], label='depredador')
    plt.xlabel('tiempo')
    plt.ylabel('población')
    plt.legend()

interact(solucion_temporal_interact,
         a=(0.01,0.5), b=(0.01,0.5),
         c=(0.01,0.5), d=(0.01,0.5),
         x0=(1,80), y0=(1,50),
         tf=(50,300));

def mapa_fases_interact(a, b, c, d, x0, y0, tf):
    
    conds_iniciales = np.array([x0, y0])

    # Condiciones para integración
    N = 800
    t = np.linspace(0, tf, N)
    
    solucion = odeint(df_dt, conds_iniciales, t, args=(a, b, c, d))
    
    x_max = np.max(solucion[:,0]) * 1.05
    y_max = np.max(solucion[:,1]) * 1.05

    x = np.linspace(0, x_max, 25)
    y = np.linspace(0, y_max, 25)

    xx, yy = np.meshgrid(x, y)
    uu, vv = df_dt((xx, yy), 0, a, b, c, d)
    norm = np.sqrt(uu**2 + vv**2)
    uu = uu / norm
    vv = vv / norm

    plt.figure("Campo de direcciones", figsize=(8,5))
    plt.quiver(xx, yy, uu, vv, norm, cmap=plt.cm.gray)
    plt.plot(solucion[:, 0], solucion[:, 1])
    plt.xlim(0, x_max)
    plt.ylim(0, y_max)
    plt.xlabel('presas')
    plt.ylabel('depredadores')
    plt.savefig('campo_direcciones.png')

interact(mapa_fases_interact,
         a=(0.01,0.5), b=(0.01,0.5),
         c=(0.01,0.5), d=(0.01,0.5),
         x0=(1,80), y0=(1,50),
         tf=(50,300));
# Install IPython magic extension in a notebook to display information about
# which versions of dependency package that was used to run the notebook.
%install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py
%load_ext version_information

%version_information, numpy, matplotlib, scipy