#!/usr/bin/env python
# -*- coding: utf-8 -*-

'''
Introducción
============

El práctico consta de dos partes:

1) Clustering, en donde se deberá implementar el algoritmo k-means

2) Clasificación, en donde se deberá entrenar un clasificador lineal

Se trabajará sobre un dataset de dígitos manuscritos a que provee la librería sklearn
(http://scikit-learn.org)

Descripcioń del dataset
=======================

Cada muestra corresponde a una imagen de 8x8 de un dígito manuscrito

----------------------------------
Clases                          10
Muestras por clase            ~180
Total de muestras             1797
Dimensionalidad                 64
Features              enteros 0-16
----------------------------------

Ejemplo de uso:

  >>> from sklearn.datasets import load_digits
  >>> digits = load_digits()
  >>> print(digits.data.shape)
  (1797, 64)
  >>> print(digits.target.shape)
  (1797,)

Los miembros mas importantes en estructura asociada al dataset son:

  data: Es una matriz de 1797x64 en donde cada fila corresponde a una muestra (imagen vectorizada)

  images: Es un arreglo de 1797x8x8 con las imágenes. Contiene escencialmente la misma informacioń
  que "data"

  target: Es un vector de 1797 elementos con las etiquetas asociadas a las muestras de entrenamiento
  (números del 0 al 9)

Se van a tomar 2/3 de las muestras para el entrenamiento y el 1/3 para test.

===================
PARTE 1: Clustering
===================

Se provee una implementación simple del algoritmo k-means. El parámetro K se fijó inicialmente a
10. Se pide:

1) Analizar y entender cada paso del código que se provee. Los comentarios en las definiciones de
funciones contienen información que le será útil a la hora de resolver el problema. Asegúrse de
entender cada paso.

2) Computar el error cuadrático medio (MSE) en el conjunto de test (ver MSE_TEST)

3) En la implementacioń actual, los centroides se inicializan con valores aleatorios (ver
INICIALIZACIÓN). Cambie esta política de forma tal que los centroides se inicializen a muestras
aleatorias tomadas del conjunto de entrenamiento, es decir, elija K muestras de forma aleatoria y
úselas como valores iniciales del algoritmo.

4) Varíe el parámetro K del algoritmo (1, 10, 100) y compute el MSE sobre el conjunto de test para
cada caso. Que observa?. A que se debe?.

   (su respuesta aquí)

Al finalizar el práctico envíe la solución (programa modificado) a: jsanchez@famaf.unc.edu.ar

'''

import numpy as np

from numpy import random, linalg

import pylab as pl

from sklearn.datasets import load_digits

# Calcula el cuadrado de la distancia Euclidea entre dos vectores
def squared_distance(a, b):
    assert(a.shape == b.shape)
    return linalg.norm(a-b, ord=2)

# Computa la matriz de distancias. Para un problema con K centroides y N muestras (datos), la matriz
# de distancias es una matriz de NxK en donde la entrada [i][j] almacena la distancia de la muestra
# "j" al centroide "i". Por lo tanto, si se considera la fila "i" de la matriz es fácil obtener a
# que centroide corresponde la muestra de índice "i" (ver K-MEANS PASO #1 más abajo)
def distance_matrix(centroids, data):
    N = data.shape[0]
    K = centroids.shape[0]
    dm = np.ndarray((N,K), dtype='float')
    for i, data_i in enumerate(data):
        for j, centroid_j in enumerate(centroids):
            dm[i][j] = squared_distance(data_i, centroid_j)
    return dm

# Calcula el error cuadrático medio (MSE) sobre un conjunto de datos. El MSE es una medida de la
# distorsión que se genera al reemplazar cada muestra del conjunto de datos con el centroide más
# cercano.
def mean_square_error(centroids, data):
    N = data.shape[0]
    dm = distance_matrix(centroids, data)
    mse = 0.0
    for i in range(N):
        mse += dm[i].min()
    return mse / float(N)

# programa principal
def main():

    # carga dataset
    digits = load_digits()

    train_data = digits.data[0:600,:]
    train_labels = digits.target[0:600]

    test_data = digits.data[600:,:]
    test_labels = digits.target[600:]

    # si show_samples es True, muestra las primeras 10 muestras del conjunto de entrenamiento
    show_samples = False
    if show_samples:
        pl.gray()
        for i in range(10):
            pl.subplot(2, 5, i)
            img = train_data[i].reshape(8,8)
            pl.matshow(img, fignum=False)
            pl.xlabel(train_labels[i])
        pl.show()

    # N: número de muestras de entrenamiento
    N = train_data.shape[0]

    # D: dimensionalidad
    D = train_data.shape[1]

    # número de clusters
    K = 100

    # INICIALIZACIÓN: K vectores de D dimensiones con valores aleatorios entre el mínimo y el máximo
    # en el conjunto de entrenamiento
    min_val = train_data.min()
    max_val = train_data.max()
    centroids = random.random_integers(min_val, max_val, (K, D)).astype('float')

    # (TAREA #2)

    # calcula error cuadrático medio inicial
    mse_prev = 0.0
    mse_curr = mean_square_error(centroids, train_data)
    mse_eps = 1.0

    niter = 1
    while niter <= 100 and mse_eps > 1e-5:

        # K-MEANS PASO #1: asignación de muestras a centroides
        distances = distance_matrix(centroids, train_data)

        assignments = [0] * N
        for i in range(N):
            assignments[i] = distances[i].argmin()

        # K-MEANS PASO #2: reestimación de centroides
        new_centroids = np.zeros(centroids.shape, dtype='float')

        counts = [0] * K
        for i in range(N):
            j = assignments[i]
            new_centroids[j] += train_data[i]
            counts[j] += 1

        for j in range(K):
            centroids[j] = new_centroids[j] / float(counts[j] + 1e-6)

        mse_prev = mse_curr
        mse_curr = mean_square_error(centroids, train_data)
        mse_eps = (mse_prev - mse_curr) / mse_curr

        print 'iter = %d, mse = %.3f (eps = %.2e)' % (niter, mse_curr, mse_eps)

        niter += 1

    # (MSE_TEST)

if __name__ == "__main__":
    main()