# -*- coding: utf-8 -*-
"""Untitled13.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1YgzGZJkGYz-1W9-dJG2N-AVO5RWqW_JZ
"""

#Tarea 4 Metodos Copmutacionales 
#Samuel David Cañas Molina - 201922001

#A Integrales

"""$f(x$) = $\int_{0}^{1}x^{-1/2}\,\text{d}x$

$g(x)$ = $\int_{0}^{\infty}e^{-x}\ln{x}\,\text{d}x$

$h(x)$ = $\int_{0}^{\infty}\frac{\sin{x}}{x}\,\text{d}x$
"""

import numpy as np
import pandas as pd
import scipy as sp
import sklearn as sl
import time

#Primer metodo Cuadraturas Gaussianas

xi = 1/2
xf = 1
xf2 = 10
n = 11

X = np.linspace(xi,xf, n)
X2 = np.linspace(xi,xf2, n)

V = np.matrix([X**i for i in range(len(X))]  )
V_inv = np.linalg.inv(V)

V2 = np.matrix([X2**i for i in range(len(X2))]  )
V2_inv = np.linalg.inv(V2)

I = np.matrix([ (   ( xf**(n+1) - xi**(n+1) ) /(n+1)   ) for n in range(len(X)) ]).T
I2 = np.matrix([ (   ( xf2**(n+1) - xi**(n+1) ) /(n+1)   ) for n in range(len(X2)) ]).T

def f(x):
    return x**(-1/2)

def g(x):
    return np.exp(-x)*np.log(x)

def h(x):
    return np.sin(x)/x


vec_f = np.matrix([ f(X[i]) for i in range(len(X))]  ).T
vec_g = np.matrix([ g(X2[i]) for i in range(len(X2))]  ).T
vec_h = np.matrix([ h(X2[i]) for i in range(len(X2))]  ).T

print("El calculo de las integrales por medio del metodo de  Cuadraturas Gaussianas es: ")
print("El calculo de la integral f(x) es: ", (vec_f.T*V_inv*I)[0,0])
print("El calculo de la integral g(x) es: ", (vec_g.T*V2_inv*I)[0,0])
print("El calculo de la integral h(x) es: ", (vec_h.T*V_inv*I)[0,0])

#Segundo Metodo Trapecio 

inicio = time.time()

n=100
x_f = np.linspace(0.000001,1,1000001)
x = np.linspace(0.000001,n,1000001)


def f(x):
    return x**(-1/2)

def g(x):
    return np.exp(-x)*np.log(x)

def h(x):
    return np.sin(x)/x

def integrate(x,f):
    calculo=0
    for i in range(len(x)-1):
        calculo = calculo + (f(x[i+1])+(f(x[i])))*abs(x[i+1]-x[i])/2
    return(calculo)

    
integral_f = integrate(x_f,f)    
integral_g = integrate(x,g)
integral_h = integrate(x,h)


print("El calculo de las integrales por medio del metodo del Trapecio es: ")
print(f"El calculo de la integral f(x) es: {integral_f}")
print(f"El calculo de la integral g(x) es: {integral_g}")
print(f"El calculo de la integral h(x) es {integral_h}")


#Tercer Metodo Simspon 

inicio = time.time()

x_f_ini = 0.000001
x_f_fin = 1
x_ini =0.000001
x_fin =100
nodos=1000001


def f(x):
    return x**(-1/2)
def g(x):
    return np.exp(-x)*np.log(x)
def h(x):
    return np.sin(x)/x
def simpson_integrate(desde, hasta, nodos, f):
    x, delta_x = np.linspace( desde, hasta, num=nodos-1 , retstep=True  )
    return (delta_x/3)*(   f(x[0]) + 2*np.sum(f(x[2:len(x)-1:2])) + 4*np.sum(f(x[1::2])) + f(x[-1]) )
  
integral_f = simpson_integrate(x_f_ini, x_f_fin, nodos, f)    
integral_g = simpson_integrate(x_ini, x_fin, nodos, g)
integral_h = simpson_integrate(x_ini, x_fin, nodos, h)

print("El calculo de las integrales por medio del metodo de Simpson es: ")
print(f"El calculo de la integral f(x) es: {integral_f}")
print(f"El calculo de la integral g(x) es: {integral_g}")
print(f"El calculo de la integral h(x) es: {integral_h}")

#B Fourier 

#Datos Tarea 3 
#Tarea 3 Metodos Computacionales 
#Samuel David Cañas Molina - 201922001 

import numpy as np
import matplotlib as mpl
import scipy as sp
import pandas as pd
import sklearn as sl
import seaborn as sns; sns.set()

from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm

# get_ipython().run_line_magic('matplotlib', 'inline')


#Tarea 3: Encuentre la regresión
#Usted recibe unos datos x y y cómo se muestran a continuación. Usted debe responder cuatro preguntas a partir de estos datos. 
#Suponga que usted tiene un modelo tal que y=f(x) más aún desconoce f.

df = pd.read_pickle("ex1.gz")
sns.scatterplot(x ='x',y ='y',data=df, label = "Datos")
plt.legend(loc="upper right")
plt.show()
df

# Punto (A) Pendiente e intercepto
# Determine la pendiente de los datos en el intervalo [0,1.5] y el valor del intercepto con el eje y. 
#Es decir, f(0)=?. ¿Cuál es el valor de r^2?

lx = df["x"]
ly = df["y"]


xt = []
yt = []
for i in range(len(lx)):
    if(0 <= lx[i] <= 1.5):
        xt.append(lx[i])
        yt.append(ly[i])

x = np.array(xt)
y = np.array(yt)


n = len(x)
sumx = sum(x)
sumy = sum(y)
sumxy = sum(x*y)
sumx2 = sum(x**2)
sumy2 = sum(y**2)
promx = sumx/n
promy = sumy/n

m = (sumx*sumy-(n*sumxy))/(sumx**2-(n*sumx2))
b = promy - m*promx

m, b

sns.scatterplot(x='x',y='y',data=df, label = "Datos")
plt.plot(x, m*x + b, color = "k", label = "Ajuste")
plt.legend(loc="upper right")
plt.title("Regresión polinomial en intervalo de 0 a 1.5")
plt.show()

sigmax = np.sqrt((sumx2/n) - promx**2)
sigmay = np.sqrt(sumy2/n - promy**2)
covar = (sumxy/n) - (promx*promy)
r2 = (covar/(sigmax*sigmay))**2


#Punto (B) Regresión polinomial
# Suponga que quiere realizar la siguiente regresión polinomial:
#y=\beta+\beta_1x+\beta_2x^2+\beta_2x^3+\beta_2x^4+\beta_2x^5.
# Plantee la función de costo que le permita calcular los coeficientes y calcule: 
#$\beta_1, \beta_2, \beta_3, \beta_4, y \beta_5. ¿Cuál es el r^2?

# Calcule f(0) y compare con los resultados anteriores

A = []
for i in range(len(lx)):
    A.append([1, lx[i], lx[i]**2, lx[i]**3, lx[i]**4, lx[i]**5])

Matriz = np.matrix(A)

b = np.array([ly]).reshape(-1,1)

Mb = np.matrix(b)


def L(k,Matriz,Mb):
    
    m,n = Matriz.shape
    X = np.matrix(k).T
    DeltaB=(Matriz*X - Mb) 
    return (DeltaB.T*DeltaB)[0,0]/m 

def dLdx(k,Matriz,Mb):
    
    m,n = Matriz.shape
    X = np.matrix(k).T
    DeltaB=(Matriz*X - Mb) 
    return (2/m)*np.array(Matriz.T*DeltaB).flatten() 

e = 1e-10

res1 = sp.optimize.minimize(fun=L,x0=np.zeros(Matriz.shape[1]), args = (Matriz,Mb), tol=1e-10)
res1['x']

#Calculo de f(0) y de cada coeficiente 
# $f(0)$ =  0.35735768,
# $\beta_1$ = -0.84260174,
# $\beta_2$ = 3.78478429,
# $\beta_3$ = -3.00338516,
# $\beta_4$ =  0.85844412,
# $\beta_5$ = -0.08305811


L(res1['x'],Matriz,Mb)


#Valor de la función de costo
# L = 0.012550626817337412

k =  np.linspace(0,4,100)

def Poli(c5,c4,c3,c2,c1,b,x):
    return  b +  c1*x +c2*(x**2) + c3*(x**3)+ c4*(x**4) +c5*(x**5) 

pre = Poli(res1['x'][5],res1['x'][4],res1['x'][3],res1['x'][2],res1['x'][1],res1['x'][0],k)

r2 = 1-np.sum((pre - ly)**2)/np.sum((ly-ly.mean())**2)
r2, pre[0]


#r^2 = 0.9111078914676154
#f(0)= 0.3573580251551316

sns.scatterplot(x ='x',y ='y',data=df, label = "Datos")
plt.plot(lx,pre,"--", color='k', label = "Ajuste")
plt.title("Regresión polinomial")
plt.legend(loc="upper right")



# Punto (C) Regresión polinomial exacta
# Resulta, que cuando se quiere hacer alguna regresión polinomial esta se puede hacer de forma exacta. 
#¿Cómo? Suponga que ud va a considerar que su problema en lugar de tener $1$ variable ($x$) tiene $n+1$, 
#siendo $n$ el orden del polinomio a ajustar. Es decir, sus nuevas variables van a ser:
#$\{x_0,\,x_1,\,x_2,\,x_3,\dots,\,x_n\}$ definiendo $x_j=x^j$. 
#Así pues, siguiendo el mismo procedimiento para la regresión lineal multidimensional que realizamos para el ejercicio de datos inmobiliarios, puede encontrar los valores de los coeficientes $\beta_1$, $\beta_2$, $\beta_3$, $\beta_4$, y $\beta_5$. Encuentre estos valores y compare con los resultados en la sección (B).
# Calcule f(0) y compare con los resultados anteriores.
# Si ud se pregunta si esto es posible la respuesta es sí. Inclusive, esto se puede extender a cualquier a cualquier conjunto de funciones, tal que $x_j=f_j(x)$, que represente un conjunto "linealmente independiente" (¡Me estoy adelantando a Fourier!). Para quienes quieran explorar algunas curiosidades matemáticas, cuando $n+1$ es igual al número de puntos o valores de $x$ (y todos diferentes) la matriz es siempre invertible y resulta ser la inversa de una matriz de Vandermonde.

lx = df["x"]
ly = df["y"]

X = np.array(lx).reshape(-1, 1)
Y = np.array(ly).reshape(-1, 1)

P = np.array([np.ones([len(lx), 1]), X, X**2, X**3, X**4, X**5]).reshape(6, len(lx)).T
coeffs = np.linalg.inv(P.T @ P) @ P.T @ Y

b, c1, c2, c3, c4, c5 = coeffs

coeffs

Ajuste = b + (c1*X) + (c2*X**2) + (c3*X**3) +  (c4*X**4) + (c5*X**5)

plt.figure()
sns.scatterplot(x='x',y='y',data = df, label = "datos")
plt.plot(X, Ajuste, c ='k', label = "Ajuste")
plt.legend(loc="upper left")
plt.title("Regresión polinomial exacta")
plt.ylabel('y')
plt.xlabel('x')
plt.show()

p = np.array(P@coeffs).flatten()
o = np.array(Y).flatten()
k = o.mean()

r2 = 1 - (np.linalg.norm(p - o)**2)/(np.linalg.norm(o - k)**2)

# r^2 = 0.9199225041318263
#f(0)=0.35644669

# Punto (D) Regresión a un modelo teórico
# Suponga que su modelo teórico es el siguiente:y=\frac{a}{\left[(x-b)^2+c\right]^\gamma}.Halle a, b, c y \gamma.
# Calcule f(0) y compare con los resultados anteriores

def f(parametro,x): 
    return (parametro[0])/((x-parametro[1])**2 + parametro[2])**parametro[3]

def Lfit(parametro,x,y): 
    deltaY=f(parametro,x) - y
    return np.dot(deltaY,deltaY)/len(y)

e = 1e-8

res2 = sp.optimize.minimize(fun=Lfit, x0=np.array([0,1,0,0]), args = (x,y), method='L-BFGS-B', tol=e)

res2["x"]


# a =0.94519587,
#b = 1.43859817
#c = 0.7390972 
#\gamma = 1.12724243

pre2 = f(res2.x,lx)

b = f(res2.x,0)

# f(0)=0.2950897972713322

plt.figure()
plt.scatter(lx,ly, label = "Datos")
plt.plot(lx,pre2, c = 'k', label = "Ajuste")
plt.legend(loc="upper left")
plt.title("Regresión a un modelo teórico")
plt.ylabel('y')
plt.xlabel('x')
plt.show()

"""Calcule la transformada rápida de Fourier para la función de la Tarea 3 (D) en el intervalo $[0,4]$ ($k$ máximo $2\pi n/L$ para $n=25$). Ajuste la transformada de Fourier para los datos de la Tarea 3 usando el método de regresión exacto de la Tarea 3 (C) y compare con el anterior resultado. Para ambos ejercicios haga una interpolación y grafique para comparar."""

# Fourier 
import numpy as np
import matplotlib as mpl
import scipy as sp
import pandas as pd
import sklearn as sl
import seaborn as sns; sns.set()

from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm

df = pd.read_pickle("ex1.gz")
sns.scatterplot(x ='x',y ='y',data=df, label = "Datos")
plt.legend(loc="upper right")
plt.show()
df

x = df["x"]
y = df["y"]

listax = []
listay = []
for i in range(len(x)):
    if x[i]<=1.5 :
        listax.append(x[i])
        listay.append(y[i])
        
x_a = np.array(listax)
y_a = np.array(listay)

def f(vec,x_a): 
    return (vec[0])/((x_a-vec[1])**2 + vec[2])**vec[3]

def L_fit(vec,x_a,y_a): 
    deltaY=f(vec,x_a) - y_a
    return np.dot(deltaY,deltaY)/len(y_a)


res2 = sp.optimize.minimize(fun=L_fit, x0=np.array([0,1,0,0]), args = (x_a,y_a), method='L-BFGS-B', tol=1e-8)
b = f(res2.x,0)

plt.figure()
plt.scatter(x,y, label = "Datos")
plt.plot(x,f(res2.x,x), c = 'r', label = "Ajuste")
plt.xlabel('x')
plt.ylabel('y')
plt.show()