import numpy as np
import PIL.Image

im = PIL.Image.open("/Users/valeriaalvarez/Documents/rh.jpeg")
col,row =  im.size
A = np.zeros((row*col, 5))
pixels = im.load()
print(pixels[187,250])
for i in range(col):
    for j in range(row):
        #print("i=%d, j=%d" % (i,j))        
        r,g,b =  pixels[i,j]
        radiohead[i*col + j,:] = r,g,b,i,j

im

im.size

print(radiohead)

radiohead.shape

"""
Explanation: Tarea 2
Valeria Alvarez Navarro
Preguntas Teóricas
1.- ¿Por qué una matriz equivale a una transformación lineal entre espacios vectoriales?
Una matriz vive en un espacio R mxn (pudiéndose descomponer en n vectores de dimensión Rm), al multiplicar una matriz por un vector de dimensión n (x elemento de Rn), se le está aplicando una transformación lineal que va a llevar a este vector al espacio Rm con lo cual se está llevando una transformación lineal entre espacios vectoriales (Rn a Rm).
2.- ¿ Cuál es el efecto de una transformación lineal de una matriz diagonal y el de una matriz ortogonal?
En el caso de la matriz diagonal, el efecto de una transformación lineal es que los vectores transformados pueden expandirse o contraerse.
En el caso de la matriz ortogonal, el efecto de una transformación lineal es que todos los ángulos y longitudes de los vectores transformados permanece sin cambio. Esto es, al multiplicar por una escalar c, los vectores van a conservar el mismo ángulo entre ellos y van a preservar la longitud original cuando la transformación los haya rotado.
3.- ¿Qué es la descomposición en valores singulares de una matriz?
Es una manera de transformar los datos contenidos en una matriz (sobre todo cuando estos tienen una dimensión grande). Esta descomposición existe para todas las matrices sin importar su dimensión (al contrario de la descomposición en valores propios).Para llevar a cabo la descomposición en valores singulares (SVD), se utliza la factorización de la matriz A de mxn en tres matrices: dos matrices ortogonales que contienen los vectores singulares (U y V) y una matriz cuasi diagonal que contiene los valores singulares de A, denominada Sigma (S). Para lograr que la matriz S sea diagonal, se le añaden ceros. Los valores de la diagonial de S son siempre no negativos y las entradas se encuentran ordenadas de mayor a menor.
La SVD es un método para obtener la aproximación más cercana de oreden k al rango de la matriz A. La elección de k es un trade off entre precisión y velocidad de procesamiento de los datos.
A se reconstruye multiplicando las tres matrices: U, S, VT (V transpuesta)
4.- ¿Qué es diagonalizar una matriz y qué representan los eigenvalores?
Diagonalizar una matriz cuadrada (A), se refiere al proceso por el cual se encuentra una matriz (P) tal que al multiplicar P inversa por A por P, se obtiene una matriz que únicamente tiene elementos en la diagonal y en todas sus demás entradas tiene ceros, la cual se llama matriz Diagonal (D). 
Las entradas de la matriz diagonal D, son los eigenvalores de A, los cuáles representan las raíces de su polinomio característico det|(A-LamdaI)|=0.
Para que la matriz A sea pueda ser diagonalizable sus n vectores deben ser linealmente independientes. Más aun, P debe ser una matriz invertible en la cual sus columnas son los n eigenvectores linealmente independientes de A.
Los eigenvalores representan el factor escalar por el cual un eigenvector ha sido multiplicado al llevar a cabo una transformación lineal; los cuales a su vez no son afectados por transformaciones lineales en cuanto a su dirección (permanecen como los originales).
5.- Intuitivamente, ¿Qué son los eigenvectores?
Son los vectores que al aplicar una transformación lineal, no son alterados en su dirección o la transformación sólo hace que se multipliquen por un escalar, lo cual tampoco cambia su dirección. 
Son los vectores que caracterizan las transformación lineal.
6.- ¿Cómo interpretas la descomposición en valores singulares como una composición de tres tipos de transformaciones lineales simples?
La primera matriz al ser ortogonal, debe llevar a cabo una rotación. Posteriormente al aplicar la matriz diagonal se lleva a cabo una expansión o contracción de ejes. Por último se vuelve a llevar a cabo una rotación ya que la tercera matriz también es ortogonal.
Preg. 2 y 3.
7.- ¿Qué relación hay entre la descomposición en valores singulares y la diagonalización?
La diagonalización sólo se puede llevar a cabo en matrices cuadradas y debe existir la matriz inversa (P, pregunta 4), ese proceso se llama también eigen descomposición. En cuando a la SVD, no importa el tamaño de las matrices y no se requiere encontrar P inversa, ya que se puede construir una matriz pseudo.inversa.
Adicionalmente, la SVD es más estable y con menos limitaciones para manejar grandes volúmenes de datos.
8.- ¿Cómo se usa la descomposición en valores singulares para dar una aproximación de rango menor a una matriz?
La descomposiciòn SVD de A, proporciona una representación  de una matriz en la forma de combinaciones lineales de vectores los cuales están ordenados por importancia (matriz S), quedarse con los "k" más importantes o los primeros es equivalente a dar una aproximaciòn de rango menor (orden k) de la matriz A; donde k<=R. Para llevar a cabo este proceso, se deben quedar las primeras k columnas de la matriz U, los primeros k valores de la matriz S, y los primeros k renglones de la matriz Vt).
9.- Describe el método de minimización por descenso gradiente
El método de minimización por descenso gradiente, es algoritmo de aproximación de mínimos locales de una función. El descenso gradiente inicia con un cojunto de parámetros iniciales y de manera iterativa se mueve hacia aquellos que minimizan la función. La manera en que lleva a cabo las iteraciones es a través de cálculo diferencial, siguiendo la dirección de la pendientes de la superficie creada por la función objetivo en la dirección opuesta a la del gradiente.
10.- Menciona 4 ejemplos de problemas de optimización (dos con restricciones y dos sin restricciones) que te parezcan interesantes como Científico de Datos


El algoritmo de asignación de registros que pretenecen a una misma clase (Data Quality)


La tasa de interés vs el riesgo de crédito


La partición (course classing) de variables tal que la pendiente de la variable dependiente sea máxima


El punto de corte de un score de originación de tarjetas de crédito sujeto a mantener una tasa de rechazo menor al 10% y una tasa de incumplimiento menor al 15%.


Parte 2: Aplicaciones de Python
Ejercicio 1
- Importar imagen
End of explanation
"""


import numpy as np
RH = np.linalg

U, s, Vh = RH.svd(radiohead, full_matrices=False)
assert np.allclose(radiohead, np.dot(U, np.dot(np.diag(s), Vh)))


R_H = np.dot(np.dot(U, np.diag(s)), Vh)
print(np.std(radiohead), np.std(R_H), np.std(radiohead - R_H))

U

s

s.shape

Vh

import scipy.linalg as sc

S = sc.diagsvd(s, 5, 5)
S

"""
Explanation: - Realizar la descomposición SVD
End of explanation
"""


U @ S @ Vh

# otra forma de hacer la verificación
rd_head=dot(U, np.dot(np.diag(s), Vh))
rd_head

"""
Explanation: - Verificar la descomposición SVD
End of explanation
"""


from scipy import *
from pylab import *
img = imread("/Users/valeriaalvarez/Documents/rh.jpeg")[:,:,0]
gray()
figure(1)
imshow(img)

m,n = img.shape
m,n

U,s,Vt = svd(img)

Vt.shape

s.shape

S = resize(s,[m,1])*eye(m,n)
S

S.shape

m,n = img.shape
U,s,Vt = svd(img)
S = resize(s,[m,1])*eye(m,n)
k = 5
imshow(dot(U[:,1:k],dot(S[1:k,1:k],Vt[1:k,:])))
show()

k = 10
imshow(dot(U[:,1:k],dot(S[1:k,1:k],Vt[1:k,:])))
show()

k = 25
imshow(dot(U[:,1:k],dot(S[1:k,1:k],Vt[1:k,:])))
show()

"""
Explanation: con el método anterior, no logré imprimir la imagen con el imshow
- Hacer la aproximación de orden k
End of explanation
"""


from numpy import *
import numpy as np

def gen_mat(i,j):
    A = floor(random.rand(i,j)*20-0) # se está haciendo una matriz aleatoria de 4x4
    b = floor(random.rand(j,1)*20-0) # este es el vector de resultados b
    return A,b

A,b= gen_mat(4,4)

A

b

#Esto sólo sirve para matrices cuadradas.
def Inversa(A):
    if((A.shape[0] == A.shape[1])):
        U,s,V=np.linalg.svd(A)
        Inversa = np.dot(np.dot(V.T,linalg.inv(diag(s))),U.T)
        return Inversa
    else:
        return "La Matriz no es cuadrada, calcula la pseudoinversa"

w=Inversa(A)
w

m,n=A.shape
m,n

U, s, V = np.linalg.svd(A)

s

U.shape

V.shape

S=eye(m,n)
S

S = resize(s,[m,1])*eye(m,n)
S

Sigma = np.zeros([U.shape[1],V.shape[0]])
Sigma

def P_Inversa(A):
    import numpy as np
    if isinstance(A, np.ndarray): 
        U, s, V = np.linalg.svd(A)
        m,n=A.shape
        S = resize(s,[m,1])*eye(m,n) 
        for i in range(S.shape[0]):
            for j in range(S.shape[1]):
                if (i == j):
                    if (s[i] == 0):
                        S[i,j] == 0
                    else:
                        S[i,j] = 1/s[i]
                        
        P_Inversa = np.dot(np.dot(V.T,S.T),U.T)
        return(P_Inversa)
    else:
        return "Errores en la especificación"

A_I=P_Inversa(A)
A_I

"""
Explanation: ¿Qué tiene que ver este proyecto con compresión de imágenes?
La descomposición de valores singulares tiene una de sus mayores aplicaciones en la compresión de imágenes. Lo que entiendo es que esta aplicación es importante conocerla pues la SVD es el mecanismo para que estas imágenes puedan ser trasmitidas de manera eficiente a través de internet con la cantidad mínima de información necesaria para que la imagen no pierda sus elemetos escenciales y con utilizando menor capacidad de almacenamiento.
En la compresión de imágenes es probablemente, donde se puede ver con mayor claridad el potencial del SVD para el análisis de datos.
Ejercicio 2
Aplicación de pseudoinversa y sistema de ecuaciones

Programar una función que dada cualquier matriz devuelva la pseudainversa usando la descomposición SVD. 
Hacer otra función que resuelva cualquier sistema de ecuaciones de la forma Ax=b usando esta pseudoinversa.
End of explanation
"""


B = np.linalg.pinv(A)
np.allclose(A, np.dot(A, np.dot(B, A)))
B

"""
Explanation: Comprobación de que el proceso anterior sale correctamente.
End of explanation
"""


def Solucion(A,b):
    import numpy as np
    if isinstance(A, np.ndarray):
        if isinstance(b, np.ndarray):
            if((A.shape[1] == b.shape[0])):#la matriz y el vector son compatibles
                if((A.shape[0] == A.shape[1])):
                    A_inv=Inversa(A)
                    x_sol = np.dot(A_inv,b)
                else:
                    A_inv=P_Inversa(A)
                    x_sol = np.dot(A_inv.T,b)
                return(x_sol)   
            else:
                return "A y b son incompatibles"
        else:
            return "Problemas con b"
    else:
        return "Problemas con A"

x = Solucion(A,b)
x

"""
Explanation: Para obtener la solución para cualquier matriz, hay que juntar las 2 funciones anteriores
End of explanation
"""


import numpy as np
A1 =np.array([[1,1],[0,0]])
A1

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

x1 = Solucion(A1,b1)
x1

"""
Explanation: Jugar con el sistema Ax=b donde A = [[1,1],[0,0]] y b puede tomar distintos valores.
a) Observar que pasa si b está en la imagen de A (contestar cuál es la imagen) y si no está (ej. b = [1, 1]).
b) Contestar, ¿la solución resultante es única? Si hay más de una solución, investigar que caracteriza a la solución directa.
c) Repetir cambiando A = [[1,1],[0, 1e-32]], ¿en este caso, la solución es única? ¿Cambia el valor devuelto de x en cada posible valor de b en el punto anterior?
End of explanation
"""


d=np.linalg.det(A1) 
d

def Solucion(A,b):
    import numpy as np
    import warnings
    if isinstance(A, np.ndarray):
        if isinstance(b, np.ndarray):
            if((A.shape[1] == b.shape[0])):#la matriz y el vector son compatibles
                if((A.shape[0] == A.shape[1]) and np.linalg.det(A) != 0):
                    A_inv=Inversa(A)
                    x_sol = np.dot(A_inv,b)
                    return(x_sol)
                elif (np.linalg.det(A)==0):
                    A_inv=P_Inversa(A)
                    x_sol = np.dot(A_inv.T,b)
                    print(x_sol) 
                    print("La Matriz A es singular")
                else:
                    A_inv=P_Inversa(A)
                    x_sol = np.dot(A_inv.T,b)
                    return(x_sol)   
            else:
                return "A y b son incompatibles"
        else:
            return "Problemas con b"
    else:
        return "Problemas con A"

"""
Explanation: Al hacer el ejercicio con A1 = [[1,1],[0,0]] y b=([[1],[1]]), se obtiene un error de que el sistema no se puede resolver porque A es una matriz singular, las columnas de A son iguales por lo que su determinante es igual a 0. Para evitar este error y que se resuelva por la pseudoinversa, se mete la condición del determinante en la función.
End of explanation
"""


x1 = Solucion(A1,b1)
x1

"""
Explanation: Se vuelve a correr el ejercicio con las modificaciones al código
Se vuelve a hacer la prueba con el vector b=[1,1]
End of explanation
"""


b2=np.array([[1],[2]])
x2 = Solucion(A1,b2)
x2

b3=np.array([[1],[0]])
x3 = Solucion(A1,b3)
x3

"""
Explanation: Se prueba con vectores b diferentes:
End of explanation
"""


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

x4 = Solucion(A2,b1)
x4

x5 = Solucion(A2,b2)
x5

x6 = Solucion(A2,b3)
x6

b4=np.array([[5],[7]])
x7 = Solucion(A1,b4)
x7

x8 = Solucion(A2,b4)
x8

b5=np.array([[5],[1e-13]])
x9 = Solucion(A1,b5)
x9

x10 = Solucion(A2,b5)
x10

"""
Explanation: Las soluciones parecen ser de la forma [x,0]
La Imagen de A, son todas las posibles combinaciones lineales de las columnas de A. 
Cambiar el valor de la posición [2,2] de la Matriz A1 por el valor 1e-32.
End of explanation
"""


import numpy as np
from pandas import DataFrame
import pandas as pd
import os

"""
Explanation: Al cambiar la entrada [2,2] de la matriz A original de 0 a 1-e32 y correr para diferentes valores de b, se observa que el resultado esta determinado por la elección de b[2,1] sin importar el valor de b[1,1]. En este caso las soluciones parecen ser de la forma [-x,x] donde x=b[2,1].
En comparación con la matriz original, el valor de b[2,1] ahora no es cero, sino que toma un valor muy grande. El pequeño cambio en la Matriz A hace que las soluciones sean muy diferenetes.
Ejercicio 3
Utilizar la paquetería pandas para trabajar el ajuste por mínimos cuadrados a un conjunto de datos.
End of explanation
"""


tabla = pd.read_csv("/Users/valeriaalvarez/Documents/tarea2.csv")
tabla

"""
Explanation: Programar un script que lea el archivo study_vs_sat.csv y lo almacene como un dataframe en pandas.
End of explanation
"""


df=pd.DataFrame(tabla)
df

"""
Explanation: Plantear un problema de optimización que intente hacer una aproximación de la forma: sat_score ~ alpha + beta*study_hoursi, con un valor para cada individuo
End of explanation
"""


import statsmodels.formula.api as sm

result = sm.ols(formula="sat_score ~ study_hours", data=df).fit()
print (result.params)

study_hours=tabla["study_hours"]
sat_score=tabla["sat_score"]

Suma_study_hours = sum(study_hours)
Suma_sat_score = sum(sat_score)
Suma_hours_score = sum(study_hours*sat_score) 
Suma_hours_2 = sum(study_hours**2)
Suma_score_2 = sum(sat_score**2)
obs = len(study_hours)


Alpha = Suma_sat_score/obs - (Beta*Suma_study_hours)/obs
Alpha

Beta = (obs*Suma_hours_score - Suma_study_hours*Suma_sat_score)/(obs*Suma_hours_2 - Suma_study_hours**2)
Beta

"""
Explanation: Plantear como un problema de optimización que intente hacer una aproximación de la forma <code>sat_score ~ alpha + beta*study_hours minimizando la suma de los errores de predicción al cuadrado. Pueden consultar este link https://en.wikipedia.org/wiki/Simple_linear_regression ¿Cuál es el gradiente de la función que se quiere optimizar (hint: las variables que queremos optimizar son alpha y beta)?

Antes de hacer las funciones, se prueba con el resultado de las librerías de python.
End of explanation
"""


def mco(b0, b1, X):
    rows = len(X)
    sat_score_e1 = np.array([b0 + b1*X[i] for i in range(rows)])
    return(sat_score_e1)

sat_score_e1=mco(Alpha, Beta, study_hours)
sat_score_e1

"""
Explanation: Programar una función que reciba valores de alpha, beta y el vector sat_score y devuelva un vector array de numpy de predicciones alpha + beta*study_hours_i, con un valor por cada individuo
End of explanation
"""


len(df["study_hours"])
col_b0 = [1 for i in range(len(x))]
X_m = np.array([col_b0,df["study_hours"]])
X_m

parametros = np.array([Alpha,Beta])
sat_score_e2=np.dot(X_m.T,parametros)
sat_score_e2

"""
Explanation: Definir 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] devuelve alpha + betastudy_hours_i en cada entrada y que entonces el problema se vuelve sat_score ~ X*[alpha,beta].
End of explanation
"""


np.linalg.pinv((X_m).dot(X_m.T)).dot(X_m).dot(sat_score)

"""
Explanation: Calcular la pseudoinversa B=(X^TX)^(-1)X^T*sat_score para obtener alpha y beta soluciones
End of explanation
"""


import matplotlib.pyplot as plt

plt.subplot(223)
plt.scatter(study_hours, sat_score_e1, label="pronosticado")
plt.scatter(study_hours, sat_score, label="observado")
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.ylabel('sat score')
plt.xlabel('study hours')
plt.title('Ajuste Modelo')




print(plt.show())

"""
Explanation: Comparar la solución para obtener los parámetros anteriormente planteada y la que se obtienen al hacel la multiplicación de matrices directa.

La solución es la misma excepto por los decimales que se redondean en la de la pseudoinversa.

Avanzado: Usen la libreria matplotlib para visualizar las predicciones  con alpha y beta solución contra los valores reales de sat_score.
End of explanation
"""


