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

# ## <center>NYU CSCI-UA 9473 Introduction to ML </center>
#    ### <center> Unsupervised Learning (Part II)</center>
# 
# <center>This week we will cover the EM algorithm, spectral clustering as well as some of the linear latent variable models </center>

# #### Exercise 1. White matter, Gray matter and cerebrospinal fluid
# 
# For this exercise, you will need to install nibabel. You can do this by typing the line 'pip install nibabel' in your terminal. 
# 
# We consider the MR image shown below. We want to segment this image into the white matter, gray matter and cerebrospinal fluid. 
# 
#  - Start by plotting the historgram of all pixel values (you can do this with 'numpy.histogram'). How does the distribution look like ? 
# 
# - Use your histogram to initialize the E-M algorithm (means and variance). The plot the resulting distributions on top ofthe histogram. Does your algorithm correctly capture the three clusters?

# In[2]:


import nibabel as nib
import numpy as np
import matplotlib.pyplot as plt

img = nib.load('T1.nii')
plt.imshow(img.dataobj[:,80,:], cmap = 'gray')

np_array = np.asarray(img.dataobj[:,80,:])

plt.show()


# In[ ]:


#### Exercise 2. Spectral clustering





# ### Exercise 2.  Image compression 
# 
# 
# #### 2.1. plotting the eigenvalues
# 
# Consider the simple image given below. We would like to compress this image by means of principal component factorization. For this purpose, we will see the columns as our feature vectors. Build the covariance matrix associated to the columns of the image (don't forget to center them first) and plot the eigenvalues associated to this covariance matrix. What do you notice?   

# In[1]:


import numpy as np
from PIL import Image
from matplotlib import pyplot as plt

img = Image.open('shapes.png')
data1 = np.array(img)

plt.imshow(data1)
plt.show()


# #### 2.2. projection onto the principal vectors
# 
# Given your result from the previous question, compress the image by computing its principal component factorization for increasing numbers of components. Compare the reconstructions.

# In[ ]:


# put your code here



# #### Exercise 3: Face recognition
# 
# In this exercise, we consider the 'Labeled Faces in the Wild' dataset from UM Amherst. We would like to train a Maximum margin classifier on this dataset (scikit learn 'svm') with a Gaussian kernel. In its original form, the dimensionality of the data is too high (if we used the raw images as our feature vectors, each vector would encode 75$\times$56 features). Although it can sometimes be a blessing, high dimensionality can also be a curse. Working in a high dimensional space means that the feature vectors will be very distant from each other which could result in a meaningless classifier. To avoid this, we will start by computing a principal component decomposition of our set of images. The SVM classifier will then be applied on the PCA representation of the training images. 

# In[6]:


import numpy as np

from sklearn.datasets import fetch_lfw_people
from sklearn.model_selection import train_test_split

lfw_people = fetch_lfw_people(min_faces_per_person=70, resize=.6)

n_samples, h, w = lfw_people.images.shape

images = lfw_people.data
targets = lfw_people.target

training_images, test_images, train_targets, test_targets = \
train_test_split(images,targets, test_size = .20, random_state=10)

plt.imshow(training_images[0,:].reshape(h,w), cmap='gray')
plt.axis('off')
plt.show()


# #### Question 3.1. PCA decomposition
# 
# - Extract the first 50, 100, 200 and 300 principal components from the image dataset and display the corresponding reconstruction for one of the faces. 
# 
# 
# - Display the plot of the reconstruction error (deviation between the true image and its PCA reconstruction for K=1 to h*w)
# 
# 
# - Finally, plot each of the principal images for K=1,...10

# In[ ]:


import numpy as np
from sklearn.decomposition import PCA

# put your code here 


# #### Question 3.2. 
# 
# Now that we have a PCA decomposition, we can learn our Max Margin classifier. We will do this with the svc model of scikit learn. finding the optimal parameters is a little tricky because of the dimension. One approach is to carry out cross validation over a parameter grid. To spare the search for the parameter values, this step is implemented below. 
# 
# - Complete the code below with the training of the Max Margin classifier (remember to split your dataset between a training and validation part).
# 
# 
# - Once you have trained the SVC, apply it on the test set and display the predictions vs actual targets. 

# In[ ]:


from sklearn.svm import SVC
from sklearn.model_selection import GridSearchCV

param_grid = {'C': [1e3, 5e3, 1e4, 5e4, 1e5],
              'gamma': [0.0001, 0.0005, 0.001, 0.005, 0.01, 0.1], }

clf = GridSearchCV(SVC(kernel='rbf', class_weight='balanced'), param_grid)


# ### II. Independent component analysis: The coktail party 

# In this exercise, we will use another approach to dimensionality reduction, known as independent component analysis (ICA). ICA is particularly useful in speech or more generally, source separation. In the classical version of this problem, known as "The coktail party problem", one is interested in recovering two distinct signals from their mixing. 

# <img src="ConversationMovie.jpg" alt="Drawing" style="width: 400px;"/>
# 
# image credit:  [The conversation](https://en.wikipedia.org/wiki/The_Conversation) 

# __Exercise II.4.a.__ Using the FastICA transform from scikit-learn, recover the two speeches from the mixed1.wav and mixed1.wav files which are given on github. 
# 
# (Hint: start by storing the two signals 'samples1' and 'samples1' into a single matrix, then pass this matrix as an input to the FastICA method of Scikit-learn) 

# In[ ]:


print(__doc__)

import os
import wave
import pylab
import matplotlib

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.io import wavfile

from sklearn.decomposition import FastICA, PCA

###############################################################################


# read data from wav files
sample_rate1, samples1 = wavfile.read('mixed1.wav')
sample_rate2, samples2 = wavfile.read('mixed2.wav')

fig, axs = plt.subplots(2)
axs[0].plot(samples1)
axs[1].plot(samples2)
plt.show()

# add your code here



# __Exercise II.4.b.__ Once you have recovered the original signals, plot them as time series. Then use the lines below to store them in new .wav file. Then play them and compare them with their mixings.

# In[ ]:


# write data to wav files

scaled1 = np.int16(S[:,0]/np.max(np.abs(S[:,0])) * 32767)
wavfile.write('recovered-1DemoJune28.wav', sample_rate1, scaled1)

scaled2 = np.int16(S[:,1]/np.max(np.abs(S[:,1])) * 32767)
wavfile.write('recovered-2DemoJune28.wav', sample_rate2, scaled2)