import matplotlib.pyplot as plt
%matplotlib inline
from skimage import io
from skimage.viewer import ImageViewer
from numpy import linalg
import numpy as np


matrix = io.imread('face.png', as_grey=True)
print(matrix.shape)

#### La clase image viewer te permite ver imagenes
#viewer = ImageViewer(img)
#viewer.show()

plt.imshow(matrix, cmap='gray')
plt.show()

"""
Explanation: SVD  imagenes

Svd de una imagen en grayscale.
End of explanation
"""


matrix.shape

u, s, vt = linalg.svd(matrix, full_matrices=False)
print(u.shape)
S = np.diag(s)
print(S.shape)
print(vt.shape)


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

"""
Explanation: Ahora la descomponemos en SVD.
End of explanation
"""


plt.imshow(m_check, cmap='gray')
plt.show()

"""
Explanation: Sea m la matrix en blanco y negro correspondiente a una imagen

Las filas de u son los eigenvectores asociados a $m$ $x$ $mT$
Las columnas de vt son los eigenvectores asociados a $mT$ $x$ $m$
$S$ es la matriz diagonal cuya diagonal contiene los valores singulares asociados a $m$

A continuación se grafica la comprobación
End of explanation
"""


k=10


u_r = u[:,:k]
S_r = S[:k,:k]
vt_r = vt[:k,:]

matrix_k = np.matmul(u_r, np.matmul(S_r,vt_r))
plt.imshow(matrix_k, cmap='gray')
plt.show()

"""
Explanation: Ahora usaremos la SVD para aproximar un rango menor de m. K será 100.
End of explanation
"""


# A es una matriz cualquiera.
A = np.array([[1,2,4],[1,4,6],[2,1,6]])
Au, As, Avt = linalg.svd(A, full_matrices=True)
AS = np.diag(As)


Ainv = np.linalg.inv(A)
Au_inv = np.linalg.inv(Au)
As_inv = np.linalg.inv(AS)
Avt_inv = np.linalg.inv(Avt)
print(Au_inv,'\n', As_inv,'\n', Avt_inv)

"""
Explanation: Esta descomposición es fundamental cuando se trata de comprimir imágenes. La matriz aproximación K mostrada previamente es una compresión de la matriz original. 

Ahora veremos la aplicación a pseudoinversa y sistemas de ecuaciones
End of explanation
"""


AS = np.diag(np.array([1/x if x != 0 else 0 for x in As]))
AS

# Verificación
print(np.dot(Au, Au_inv))
print(np.dot(As, As_inv))
print(np.dot(Avt, Avt_inv))

"""
Explanation: Para Au y Avt, las transpestas son las inversas, en cuanto a S, por ser diagonal, los reciprocos de su diagonal, son la diagonal principal de la inversa de S, siempre y cuando distitnos de cero.
End of explanation
"""


print(np.dot(Avt_inv,np.dot( As_inv,Au_inv)) )
# Ahora todo va a coincidir
print(np.linalg.pinv(A))
print(Ainv)

"""
Explanation: Por que la definición es muy clara. Si la SVD de $A = U\sum V^{T}$
La pseudoinversa de A, llamemosla $A^{+} = V\sum^{-1}U^{T} $
Multiplicar de esta forma sería un error   $U^{T}\sum^{-1}V$
End of explanation
"""


def get_pseudoinverse(matrix):
    
    """
    Get pseudo inverse of a matrix
    """
    
    Au, As, Avt = linalg.svd(matrix, full_matrices=False)
    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_matrix = np.dot(Avt_inv,np.dot( As_inv, Au_inv))
    
    return pseudo_matrix

#Verificamos
A = np.array([[1,2,4],[1,4,6],[2,1,6]])
Ainv = get_pseudoinverse(A)
iden = np.dot(A,Ainv)
iden

def solve_system(coefficient_matrix, image_vector):
    
    pseudoinverse = get_pseudoinverse(coefficient_matrix)
    domain_vector = np.dot(pseudoinverse,image_vector)
    
    return domain_vector
    

### Prueba de calidad 
B = np.array([[-10,9],[10,5]])
b = np.array([-9,-5])
print(B,'\n',b)
x = solve_system(B,b)
print(x)
np.dot(B,x)

###### Funciona!!!!!!!!

"""
Explanation: Programar una función que dada cualquier matriz devuelva la pseudainversa usando la descomposición SVD.
End of explanation
"""


A = np.array([[1,1],[0,0]])
A

"""
Explanation: Ejercicio 3
End of explanation
"""


b = np.array([1,1])

"""
Explanation: La imagen de A es el siguiente conjunto $  b \hspace{3mm} x \hspace{3mm} \left(\begin{array}{ccc}  1 \ 0 \end{array} \right)$ con $b \in \Bbb R$
End of explanation
"""


solve_system(A,b)

"""
Explanation: Evidentemente b no esta en la imagen
La solución no trivial al sistema homogeneo es $ x \hspace{3mm} = \hspace{3mm} \left(\begin{array}{ccc}  a \ -a \end{array} \right)$ con $a \in \Bbb R$
End of explanation
"""


np.linalg.solve(A,b)

"""
Explanation: Por la pseudoinversa se encuentra una solución que no es otra cosa más que el vector $b$ proyectado en el espacio generado por la matriz A.
End of explanation
"""


A = np.array([[1,1], [0,1e-32]])
A

np.linalg.det(A)

"""
Explanation: La matriz A es singular por definición no tiene inversa. Sin embargo, por
el método de la pseudo inversa podemos aproximarla solución aún cuando no esté en el espacio de la imagen de x o, equivalentemente, el espacio generado por $A$ el método que diseñamos, obtiene la pseudoinversa y aproxima una solución proyectando el vector $b$ sobre el espacio generado por la matriz A. 


Si usaramos el tradicional de python o R, de obten la inversa, alzaría una excepción del tipo
Error, la matriz es singular.


Ahora $ A =  \left( \begin{array}{ccc}
 1 & 1 \
 0 & 1e-32 \end{array} \right)$
End of explanation
"""


b = np.array([2,0])
c = np.array([1,1])

solution_pseudo = solve_system(A,b)
solution_inv = np.linalg.solve(A,b)
print(np.dot(A,solution_pseudo))
print(np.dot(A,solution_inv))


solution_pseudo_c = solve_system(A,c)
solution_inv_c = np.linalg.solve(A,c)
print(solution_pseudo_c)
print(solution_inv_c)
#print(np.dot(A,solution_pseudo_c))
print(np.dot(A,solution_inv_c))

"""
Explanation: El determinante es distinto de cero, la inversa existe.
End of explanation
"""


print(solution_pseudo_c)
print(np.dot(A,solution_pseudo_c))

"""
Explanation: El método de la inversa sí encontró una solución, sin embargo, no es el vector (1,1)
End of explanation
"""


import pandas as pd

base =pd.read_csv('study_vs_sat.csv')
base.head(5)

"""
Explanation: El método gradiente descendiente
End of explanation
"""


def predict_score( alpha, beta, study_hours):
    
    predictions = np.array([alpha + beta * study_time for study_time in study_hours])
    
    return predictions
    

"""
Explanation: Planteamiento del problema de minimización
Sea $ y $ sat_score y la queremos explicar mediante una relación lineal con $ x $ study_hours, 
entonces queremos encontrar $ \alpha , \beta $ tal que para cada observación  aproximamos $ y_i$ con   $ \hat{y}_i = \alpha_i + \beta x_i $
La función de error que deseamos minimizar:  $ E  = \sum\limits_{i=1}^n (y_i - \hat{y}i )^2 = \sum\limits{i=1}^n (y_i - \alpha_i + \beta x_i )^2 $
Encontrar $ \alpha $ y $ \beta $ tal que minimizen E. 
Las condiciones de primer orden implican obtener las derivadas parciales con el vector gradiente e igualarlas 0.
Vector Gradiente $$ \nabla F(\alpha,\beta) = (\frac{\partial E}{\partial \alpha},\frac{\partial E}{\partial \beta}) $$
Condiciones de Primer Orden
$$ \frac{\partial E}{\partial \alpha} = -2 * \sum\limits_{i=1}^n (y_i - \alpha_i + \beta x_i )  = 0 $$
$$ \frac{\partial E}{\partial \beta} =  -2 * \sum\limits_{i=1}^n (y_i - \alpha_i + \beta x_i )*x_i = 0 $$
La siguiente función recibe los parámetros alfa, beta y un numpy array de una dimensión con las horas de estudio y regresa un numpy array con predicciones de la forma $ y_i = \alpha + \beta x_i $
End of explanation
"""


base =pd.read_csv('study_vs_sat.csv')
base.head(5)
study_hours = base.iloc[:,0]
study_hours = study_hours.values

alpha, beta = 2,3

X_array = np.empty( shape=[20,2], dtype=int)

for i in range(20):
    
    X_array[i] = [1, study_hours[i]]

X_array * [alpha, beta]

"""
Explanation: Definan un numpy array X de dos columnas, la primera con unos en todas sus entradas y la segunda con la variable study_hours. Observen que X[alpha,beta] nos devuelve alpha + betastudy_hours_i en cada entrada y que entonces el problema se vuelve sat_score ~ X*[alpha,beta]
End of explanation
"""


#La pseudoinversa de X es la siguiente 
X_pseudoinv = get_pseudoinverse(X_array)
X_pseudoinv

sat_score = base['sat_score'].values
solution_vector = np.dot(X_pseudoinv, sat_score)

print('Alfa óptima con pseudoinversa es :',solution_vector[0])
print('\n'*1)
print('Beta óptima con pseudoinversa es :', solution_vector[1])

"""
Explanation: Calculen la pseudoinversa X^+ de X y computen (X^+)*sat_score para obtener alpha y beta soluciones
End of explanation
"""



J = np.linalg.inv(np.dot(X_array.transpose(), X_array))
K = np.dot(J, X_array.transpose())
sol_vector = np.dot(K, sat_score)
print('Alfa óptima con inversa tradicional es :',sol_vector[0])
print('\n'*1)
print('Beta óptima con inversa tradicional es :', sol_vector[1])
print('La solución es la misma que con la pseudoinversa')

"""
Explanation: Comparen la solución anterior con la de la fórmula directa de solución exacta (alpha,beta)=(X^tX)^(-1)X^t*sat_score
Resolvamos por partes, queremos 
$((X^{t}X))^{-1}X^{t}$ * study_hours
Sea $ J = ((X^{t}X))^{-1}$
End of explanation
"""


