import matplotlib.pyplot as plt
import matplotlib as mp
import numpy as np
import math

# Esta linea hace que los graficos aparezcan en el notebook en vez de una ventana nueva.
%matplotlib inline

"""
Explanation: Catedra 01
Primera tarea sera anunciada hoy.

Entrega el segundo miercoles, despues de haber tenido 2 auxiliares.

La primera parte del curso

Derivadas e Integrales numericas.
Manejo de errores debido a precision.
End of explanation
"""


# Cambia el tamano de los ticks (los puntos de los ejes)
mp.rcParams['xtick.labelsize'] = 13
mp.rcParams['ytick.labelsize'] = 13

"""
Explanation: Nota:

matplotlib.org -> documentacion de pyplot.
End of explanation
"""


# Parametros iniciales.
e = 1.
k_factorial = 1.
N_max = 10
e_vs_n = [e] # Lista que va a contener los elementos de la serie.
for i in range(1,N_max): # Ciclo que calcula los elementos de la serie y los suma.
    k_factorial *= i
    e += 1. / k_factorial
    e_vs_n.append(e)

# Instruccion para imprimir nuestra aproximacion a e en cada iteracion.
for i in range(N_max):
    print(i, e_vs_n[i])

plt.plot(range(N_max), e_vs_n) # Genera un grafico de e_vs_n v
plt.axhline(math.e, color='0.5') # Genera una linea en 'e'
plt.ylim(0,3) # Cambia los limites del eje y del grafico
# Cambia los labels del grafico.
plt.xlabel('N', fontsize=20)
plt.ylabel('$e_{N}$', fontsize=20)

# Se hace un ciclo for dentro de una lista. i.e. Comprehension list
diferencia = [math.fabs(e_i - math.e) for e_i in e_vs_n] # Nota. fabs convierte a float y luego calcular abs
# Se hace un grafico con escala logaritmica en el eje y.
plt.plot(range(N_max), diferencia)
plt.yscale('log')
plt.xlabel('N', fontsize=20)
plt.ylabel('$e_N - e_{real}$', fontsize=20)

"""
NOTA:
Lo anterior es equivalente a hacer
plt.semilogy(range(N_max), diferencia))
plt.xlabel('N', fontsize=20)
plt.ylabel('$e_N - e_{real}$', fontsize=20)
"""

"""
Explanation: Manejo de errores:
Error de truncacion

Calculo del numero e a traves de la expansion de Taylor: $\sum_{k=0}^{\infty}\frac{x^{k}}{k!}$

Forma basica:
End of explanation
"""


import numpy as np

epsilon = np.logspace(-1, -15, 14, base=10.)
print(epsilon)

print(type(e_vs_n))
print(type(epsilon))

dsindx = (np.sin(1.+epsilon) - np.sin(1.)) / epsilon
print(dsindx)

plt.semilogx(epsilon, dsindx - np.cos(1.), 'o')
plt.axhline(0, color='0.8')
plt.xlabel('epsilon', fontsize=20)
plt.ylabel('$\\frac{d}{dx}\\sin(x) - \\cos(x)$', fontsize=20)
_ = plt.xticks(epsilon[::2])

"""
Explanation: Efectos de precision sobre integrales
Error de Redondeo
Calculando la derivada de $\sin(x)$
Ahora usamos numpy
End of explanation
"""


def integral_trap(f, intervalo):
    '''
    Integracion numerica utilizando el metodo o regla trapezoidal.
    Note que el intervalo debe estar equiespaciado para este metodo.
    
    Parameters
    ----------
    f : numpy.ndarray
        Integrando.
    intervalo : numpy.ndarray
        El intervalo de integracion.
    
    Returns
    -------
    res : double
        Resultado de la integracion.
    '''
    res = f[0]+f[len(intervalo-1)] # Se inicializa la variable que contendra el resultado final.
    dx = intervalo[1] - intervalo[0] # Se escribe el delta x.
    for i in range(len(intervalo)-1): # Se escribe el ciclo que calculara la integral numericamente.
        res += 2*(f[i]+f[i+1])
    
    res *= dx/2.
    
    return res

"""
Explanation: Catedra 02 (Continuacion de Catedra 01)
Integrales Numericas
$$g(x) = \int_{a}^{x} f(x')dx'$$
Equivale a resolver una ecuacion diferencial:
$$g'(x)  = f(x)$$ con condiciones de borde: 
$$\int_{a}^{x}g'(x')dx' = \int_{a}^{x}f(x')dx'$$
$$g(x)-g(a) = \int_{a}^{x}f(x')dx'$$
Integracion Numerica Directa
Regla Trapezoidal
Idea: Aproximar el area de la curva como el area bajo un trapecio entre dos puntos de la funcion

$$
\int_{x_{0}}^{x_{0}+\Delta x}f(x')dx'
$$
Haciendo una expansion de taylor sobre la integral se obtiene:
$$
\int_{x_{0}}^{x_{0}+\Delta x}\left[
    f(x_{0}) + f'(x_{o})(x'-x_{0}) + \frac{1}{2}f''(x_{0})(x'-x_{0})^{2} + \ldots
\right]dx'
$$
$$
= f(x_{0}\Delta x) + \frac{f'(x_{0})(x-x_{0})^{2}}{2}\biggr\vert_{x_{0}}^{x} + \frac{1}{2}f''(x_{0})\frac{(x-x_{0})^{3}}{3}\biggr\vert_{x_{0}}^{x} + \ldots
$$
$$
= f(x_{0})\Delta x + f'(x_{0})\frac{\Delta x^{2}}{2} + \frac{1}{6}f''(x_{0})\Delta x^{3} + \ldots
$$
$$
= \frac{\Delta x}{2}\left[
    f(x_{0}) + \left(
        f(x_{0}) + f'(x_{0})\Delta x + f''(x_{0})\frac{\Delta x^{2}}{2} + \ldots
    \right) - f''(x_{0})\frac{\Delta x^{2}}{6}
\right]
$$
Note que el error de truncacion es del mismo orden en $\Delta x$ que el ultimo termino previo a truncar.
$$
\int_{x_0}^{x_{0}+\Delta x}f(x')dx' = \frac{\Delta x}{2}\left[
    f(x_{0}) + f(x_{0}+\Delta x) + O(\Delta x^2)
\right] = \left(
    f(x_{0}) + f(x_{0}+\Delta x)
\right)\frac{\Delta x}{2} + O(\Delta x^3)
$$
Donde $O(\Delta x^3) \sim -f''(x_{0})\frac{\Delta x^3}{12}$
(graficos van aca)
dividir [a,b] en $N = \frac{b-a}{\Delta x}$ tramos.
$$
\int_{a}^{b}f(x)dx \sim \sum_{i=0}^{N-1}\left[
    \frac{f(a+\Delta xi)+f(a+\Delta x(i+1))}{2}\Delta x + O^{*}(\Delta x^3)
\right]
$$
Nos interesa encontrar una cota superior para nuestro error. Como estamos sumando N tramos, estamos incluyendo N errores, y dado que el error de cada tramo va como $\Delta x^3$ y N va como $\Delta x$, el error finalmente queda como $O(\Delta x^2)$
$$
\sim \left[
    \sum_{i=0}^{N-1}\frac{f(a+\Delta xi)+f(a+\Delta x(i+1))}{2}\Delta x
\right] + NO^{*}(\Delta x^3) \rightarrow O(\Delta x^2)
$$
$$
\int_{a}^{b} f(x)dx \sim \frac{f(a)\Delta x}{2} + \sum_{i=1}^{N-1}f(a+i\Delta x)\Delta x + \frac{f(b)\Delta x}{2} + O(\Delta x^2)
$$
La regla compuesta asumiendo $\Delta x$ constante
Una implementacion simple de este algoritmo es la siguiente:
End of explanation
"""


def simpson(f, a, b, n):
    '''
    Integracion numerica utilizando el metodo o regla trapezoidal.
    Note que el intervalo debe estar equiespaciado para este metodo.
    
    Parameters
    ----------
    f : function
        Integrando.
    a : int or double
        Valor inical del intervalo sobre el cual se va a integrar.
    b : int or double
        Valor final del intervalo sobre el cual se va a integrar.
    n : int
        Numero de particiones del intervalo. Debe ser par.
    
    Returns
    -------
    res : double
        Resultado de la integracion.
        
    Raises
    ------
    ValueError
        Cuando n es impar.
    '''
    if n % 2: # El algoritmo funciona solo para n par. Aca nos aseguramos de recibir un n par.
        raise ValueError("n debe ser par (se recibio n=%d)." % n)
        
    dx = float(b - a) / n # Se inicializa el dx.
    i = a # Se inicializa en a un iterador
    res = 0. # Se inicializa res, en donde se guardara el resultado de la integral.
    sumaPar = sumaImp = 0. # Se inicializan las sumas parciales.
    k = 0
    
    while k < n/2-1: # El ciclo que calcula la integral.
        i += dx
        sumaImp += f(i)
        i += dx
        sumaPar += f(i)
        k += 1
        
    sumaImp += f(i+dx) # Se agrega el ultimo termino de la suma impar.
    res += dx/3.*(f(a) + 4 * sumaImp + 2 * sumaPar + f(b))
    
    return res

"""
Explanation: Regla de Simpson
Idea: Mas evaluaciones (terminos de la expansion) a cambio de mayor precision.

Se hace una suerte de "doble Taylor"
$$
\int_{x_{0}}^{x_{0}+2\Delta x}f(x)dx = \int_{x_{0}}^{x_{0}+2\Delta x}\left[
    f(x_{0} + f'(x_{0})(x-x_{0}) + f''(x_{0})\frac{(x-x_{0})^{2}}{2} + f'''(x_{0})\frac{(x-x_{0})^{3}}{6} + f^{\text{iv}}(x_{0})\frac{(x-x_{0})^{4}}{24} + \ldots
\right]
$$
$$
= f(x_{0})2\Delta x + f'(x_{0})\frac{(2\Delta x)^{2}}{2} + \frac{f''(x_{0})}{2}\frac{(2\Delta x)^{3}}{3} + \frac{f'''(x_{0})}{6}\frac{(2\Delta x)^{4}}{4} + \frac{f^{\text{iv}}(x_{0})}{24}\frac{(2\Delta x)^{5}}{5}
$$
$$
= \frac{\Delta x}{3}\left[
    f(x_{0}) + 4f(x_{0}) + 4f'(x_{0})\Delta x + 4 f''(x_{0})\frac{\Delta x^2}{2} + 4f'''(x_{0})\frac{\Delta x^3}{6} \
    + f(x_{0}) + f'(x_{0})(2\Delta x) + f''(x_{0})\frac{(2\Delta x)^{2}}{2} + f'''(x_{0})\frac{(2\Delta x)^{3}}{6} + O(\Delta x^4)
\right]
$$
$$
= \frac{\Delta x}{3}\left[
    f(x_{0}) + 4f(x_{0}+\Delta x) + f(x_{0} + 2\Delta x)
\right] + O(\Delta x^5)
$$
Note que para la regla compuesta el orden del error va como $\Delta x^4$
Regla compuesta
$$
\int_{a}^{b}f(x)dx = \left[
    \sum_{i=1}^{N-1}\frac{\Delta x}{3}\left[
        f_{i-1} + 4f_{i} + f_{i+1}
    \right]
\right] + O(\Delta x^4)
$$
Donde $f_{i} = f(a+i\Delta x)$
Note que este metodo tambien asume $\Delta x$ constante.
(grafico aca)
Note que la regla anterior tambien puede ser escrita como:
$$
\int_{a}^{b}f(x)dx = \frac{\Delta x}{3}\left(
    f(a) + 4\sum_{i=1}^{N/2-1}f_{2i+1}+2\sum_{i=1}^{N/2-1}f_{2i} + f(b)
\right)
$$
Una implementacion simple de este algoritmo es la siguiente:
End of explanation
"""


