egImage = "https://upload.wikimedia.org/wikipedia/commons/7/7d/Dog_face.png"
img = io.imread(egImage, as_grey = True)
plt.imshow(img, cmap = cm.gray)
plt.show()

print (type(img))
print (str(img.shape))

"""
Explanation: I) Aplicación de la SVD a la compresión de imágenes
Recibir el path de un archivo de imagen png (i) y convertirlo en una matriz numerica con coeficientes en (0,1) (ii)
End of explanation
"""


img2 = img[5:70,0:70]
plt.imshow(img2, cmap = cm.gray)
plt.show()

?imread

img.shape

img2.shape

"""
Explanation: Le cortamos unos pedazos a la matriz y generamos una nueva
End of explanation
"""


?linalg.svd #Full SVD reconstruction means diagonlization? Reduced SVD reconstruction means SVD?

"""
Explanation: Realizar y verificar la descomposición SVD (iii)
End of explanation
"""


u, s, vt = linalg.svd(img, full_matrices=True) #img es simetrica
print(u.shape)
S = np.diag(s)
print(S.shape)
print(vt.shape)


"""
Explanation: Caso A simetrica
End of explanation
"""


S=np.diag(s)
img_rec = np.matmul(np.matmul(u,S),vt)
print(img_rec.shape)

"""
Explanation: Comprobamos que el rango de la aproximación es el mismo
End of explanation
"""


u, s, vt = linalg.svd(img2, full_matrices=False) #img2 NO es simetrica
print(u.shape)
S = np.diag(s)
print(S.shape)
print(vt.shape)

"""
Explanation: Caso A NO simetrica
End of explanation
"""


img2_rec = np.matmul(np.matmul(u,S),vt)
print(img2_rec.shape)

"""
Explanation: Comprobamos que el rango de la aproximacion es el mismo
End of explanation
"""


plt.imshow(img2_rec, cmap = cm.gray)
plt.show()

"""
Explanation: Y, graficamos la aproximación SVD...
End of explanation
"""


k=1000


u_k = u[:,:k]
S_k = S[:k,:k]
vt_k = vt[:k,:]
img_k = np.matmul(u_k,np.matmul(S_k,vt_k))
plt.imshow(img_k, cmap = cm.gray)
plt.show()

"""
Explanation: Aproximaciones de la SVD con distintos valores de k (iv)...
End of explanation
"""


#Primero hagamosla para una matriz particular y despues veamos el tema de la generalización. Sea A una matriz cualquiera:
A = np.array([[1,2,3],[4,5,6],[7,8,9]])
Au, As, Avt = linalg.svd(A, full_matrices=True)
#As = np.diag(As)
A_inv = np.linalg.inv(A)
Au_inv = np.linalg.inv(Au)
As_inv = np.diag(np.array([1/x if x != 0 else 0 for x in As]))
Avt_inv = np.linalg.inv(Avt)
print('                 A_inv')
print(A_inv)
print('                 Au_inv') 
print(Au_inv)
print('                 As_inv')
print(As_inv)
print('                 Avt_inv')
print(Avt_inv)

#Verificacion del calculo de las inversas
print('         Au y Au inversa')
print(np.dot(Au,Au_inv))
print('         As y As inversa')
print(np.dot(As,As_inv))
print('         Avt y Avt inversa')
print(np.dot(Avt,Avt_inv))

"""
Explanation: La SVD sirve para almacenar de forma mas compacta imagenes ya que la matriz de aproximación de grado k contiene la información necesaria para replicar la imagen original
II) Aplicaciones de la pseudoinversa en resolucion de sistemas de ecuaciones
Programar una funcion que dada cualquier matriz devuelva la pseudoinversa usando la descomposición SVD (i). Hacer otra funcion que resuleva sistemas de ecuaciones de la forma $\mathbf b=A\overrightarrow{x} $ usando la pseudoinversa (ii)
End of explanation
"""


#NO ME DA PREGUNTAR MAÑANA
A_pseudo_inv = np.dot(Avt_inv,np.dot(As_inv,Au_inv))
print("           Metodo propio")
print (A_pseudo_inv)
print("           Metodo numpy ")
print(A_inv)


"""
Explanation: La definición de la psuedo inversa es la siguiente. Partiendo de la SVD $ \mathbf A = U \Sigma V^{T} $ por lo que la pseudoinversa de A será $ A^{+}=V \Sigma^{-1} U^{T} $
End of explanation
"""


def get_pseudoinverse(matriz):
    "Regresa la pseudo inversa de la matriz dada"
    Au, As, Avt = linalg.svd(matriz, full_matrices=True)
    #As = np.diag(As)
    Au_inv = Au.transpose()
    As_inv = np.diag(np.array([1/x if x!=0 else 0 for x in As]))
    Avt_inv = Avt.transpose()
    
    pseudo_matriz = np.dot(Avt_inv,np.dot(As_inv, Au_inv))
    
    return pseudo_matriz

"""
Explanation: La función:
End of explanation
"""


#Verifiquemos que sirve
A = np.array([[1,2,4],[1,4,6],[2,1,6]])
Ainv = get_pseudoinverse(A)
iden = np.dot(A,Ainv)
print("La matriz A por su pseudoinversa debe ser la identidad")
iden

"""
Explanation: Verifiquemos que sirve...
End of explanation
"""


def resuelve_ec(A,b):
    #Primero obtengamos la pseudoinversa con el metodo get_pseudoinverse recien creado
    A_inv = get_pseudoinverse(A)
    x = np.dot(A_inv, b)
    return x

"""
Explanation: Función para solucionar sistemas de ecuaciones lineales de la forma: $\mathbf b=A\overrightarrow{x} $ usando la pseudoinversa de A talque $ \mathbf A^{+}b=\overrightarrow{x}$
End of explanation
"""


A = np.array([[1,2],[2,1]])
b = np.array([3,3])
print('\n  matriz A =\n',A,'\n  b=',b)
x = resuelve_ec(A,b)
print('\n  x = ',x)
y = np.dot(A,x)
print('\n  Ax = ',y)

"""
Explanation: Verifiquemos que sirve...
End of explanation
"""


# Usemos nuestras funciones para jugar con distintos valores de b
A = np.array([[1,1],[0,0]])
b = np.array([1,1])
x = resuelve_ec(A,b)
print(x)

"""
Explanation: Jugar con el sistema $ \mathbf Ax=b $ donde $ \mathbf A = [[1,1],[0,0]] $ y $ \mathbf b $ puede tomar distintos valores. ¿Cúal es la imagen de A?¿Qué pasa con la solución a la ecuación si $ \mathbf b $ no esta en la imagen de A, la solución resultante es única?

$ \mathbf Ax=b \rightarrow Img(A) = \left { {y \in \mathbb{R}^{2} : \exists x : Ax=y} \right }  $ <p>Por lo que la solución al sistema consiste en encontrar los valores de $\mathbf  x_{1}, x_{2} $ para los cuales se cumple que:

</p>
$$ \begin{pmatrix}
 a, b \ 
 c, d 
\end{pmatrix} \begin{pmatrix}
 x_{1} \ 
 x_{2} 
\end{pmatrix} =  \begin{pmatrix}
 b_{1} \ 
 b_{2} 
\end{pmatrix} $$
$$ \rightarrow 
ax_{1}+bx_{2} = b_{1} \cap  cx_{1}+dx_{2} = b_{2} $$
<p>Por lo que si $ \mathbf c=0,d=0,b_{2}=1 $ $ \mathbf b \notin Img(A) $ </p>
<p>Pues tendría suceder:</p>
$$ \mathbf x_{1}+x_{2}=1 $$ 
<p> y </p>
$$ 0=1 $$
 $$ \rightarrow \mathbf Img(A) =\left {   {y \in \mathbb{R}^{2} : y_{1} \in \mathbb{R} \cap 1 = 0} \right } $$ 
<p> Si $\mathbf b $ no se encuentra en la $\mathbf Img(A) $  es necesario encontrar una proyección de este vector en el espacio generado por las columnas de la matriz $\mathbf A$, generalmente se usa a la matriz de proyección:</p>
$$\mathbf P_{x} = x \left[  xx^{T}  \right]^{-1} x^{T} $$
<p>Lo que implica que: </p>
$$  \widehat{b} \mathbf = AP_{x} = A \left( x \left[  xx^{T}  \right]^{-1} x^{T} \right) $$
<p>Para resolver este sistema se utiliza el método de mínimos cuadrados y la solución del método es única</p>

Matriz A con determinante igual a 0:
<p>En este caso la solución es {$  \mathbf x \in \mathbb{R}^{2} : x_{1}+x_{2}=1 \cap $ ¡$ 1=0 $! $ $} </p>
<p>Que significa que el vector $ \mathbf b \notin Img(A) $
End of explanation
"""


A = np.array([[1,1],[0,1e-32]])
b = np.array([1,1])
y = resuelve_ec(A,b)
x = np.linalg.solve(A,b)
print(y)
print(x)


"""
Explanation: Matriz A con determinante distinto de 0:
<p>En este caso la solución es {$  \mathbf x \in \mathbb{R}^{2} : x_{1}+x_{2}=1 \cap  1=e^{-32}x_{2}  $} </p>
<p>Que significa que el vector $ \mathbf b \in Img(A) $
End of explanation
"""


import pandas as pd
from sklearn.linear_model import LinearRegression

base = pd.read_csv("study_vs_sat.csv", header=0)
print(base.head())
print(base.dtypes)

"""
Explanation: Importar un archivo csv usando la libreria pandas y resolver un modelo de minimos cuadrados
End of explanation
"""


base_reg = np.matrix(base)
X, Y = base_reg[:,0], base_reg[:,1]
mdl = LinearRegression().fit(X, Y)
m = mdl.coef_[0]
b = mdl.intercept_
print("E(Y|X): E(y|x) = {0}E(x) + {1}".format(m,b))

"""
Explanation: <p>Se comienza por postular a la forma del valor esperado del sat_score $ \mathbf y $ condicional a las study_hours $ \mathbf x $ como una relación lineal (MPL) en el cual si se puede asumir que $ \mathbf E(\vec{u} \mid \vec{x})=0 $ , se tiene: </p>
$$ \vec{y} = \vec{ \beta } \vec{x} + \vec{u} \Rightarrow E( \vec{y} \mid \vec{x} ) = \vec{ \beta} \vec{x} $$
<p>Con: </p>
$$ \vec{ \beta } = (1, \beta_{1}, ... , \beta_{n}) $$
<p>Por lo que, el problema que se desea resolver es el siguiente: </p>
$$ \mathbf min_{ \vec{ \beta }} E(UU_{T}) = min_{ \vec{ \beta }} E(( \vec{y} - E( \vec{y} \mid \vec{x}))^2) $$
<p>Con función objetivo: $ \mathbf E(( \vec{y}- \vec{ \beta } \vec{x})( \vec{y}- \vec{ \beta } \vec{x})^{T}) $ </p>
<p>El gradiente de este problema se encuentra al derivar con respecto a $ \vec{ \beta } $ la funcion objetivo: </p>
$$ \bigtriangledown_{ \vec{ \beta } } = \frac{\partial  E(( \vec{y}- \vec{ \beta } \vec{x})( \vec{y}- \vec{ \beta } \vec{x})^{T}) }{\partial \vec{ \beta }} $$
<p>La C.P.O de este problema obligan a que en el optimo se cumpla : </p>
$$ \bigtriangledown_{ \vec{ \beta } } = \vec{0} $$
<p>La libreria sklearn ya contiene el modulo de regresion lineal que permite encontrar los coeficientes buscados a travéz del método fit </p>
End of explanation
"""


def MejorPredictorLineal(alpha, beta, x):
    "Regresa el resultado de la ecuacion postulada: Mejor Predictor Lineal"
    X = np.array(x)
    y_arr = b + m*X 
    return y_arr


"""
Explanation: Construye una función que reciba valores de alpha, beta y horas de estudio para devolver en vector de numpy con la predección del modelo de mínimos cuadrados correspondiente
End of explanation
"""


ListaDeNuevasCalificaciones = np.array([[6, 7, 8, 9, 7]])
ResultadosPredecidos = MejorPredictorLineal(353.16, 25.33, ListaDeNuevasCalificaciones)

print('\n Resultados Predecidos  E(Y|X) =', ResultadosPredecidos)


import itertools
import operator
print(base_reg[:,0].size)
horas_est = np.array([base_reg[:5,0]])
print(horas_est)
item = print(horas_est.size)
aux = horas_est.tolist()
merged = list(itertools.chain(*aux))
print(merged)
b = np.array([[1, 1, 1, 1, 1],aux])
print(b)
b

sum(aux,[])

base_reg_1 = np.matrix(base)+ b +1

alpha, beta = base_reg[:,0], base_reg[:,1]


"""
Explanation: Veamos que sucede si hay obtenemos una muestra distinta de individuos y les aplicamos el Mejor Predictor Lineal que postulamos:
End of explanation
"""


