# -*- coding: utf-8 -*-
"""HW6_Q4.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1zR1N3xFVakUv7q8it72PAPzxg0xXOd_d
"""

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.decomposition import PCA
import math

"""a. Load the data and construct a table with 9 columns containing the numerical ratings. (Ignore the last 5 columns – they consist auxiliary information such as longitude/latitude, state, etc.)"""

data = pd.read_csv('places.txt', delim_whitespace=True)
# ignore last 5 columns
data_set = pd.DataFrame(data, columns=['Climate', 'HousingCost', 'HlthCare', 'Crime',  'Transp', 'Educ', 'Arts', 'Recreat', 'Econ'])
data_set.head(6)

"""b. Replace each value in the matrix by its base-10 logarithm. (This pre-processing is done for convenience since the numerical range of the ratings is large.) You should now have a data matrix X whose rows are 9-dimensional vectors representing the different cities."""

mat = np.array(data_set).astype('float')
data_z_score = (mat - np.mean(mat, axis=0)) / np.var(mat, axis=0) ** 0.5
data_z_score = data_z_score - np.mean(data_z_score, axis=0)
# do z-score
print(f'data after z-score process is:\n{data_z_score}\n')

for i in range(len(mat)):
  for j in range(len(mat[0])):
    mat[i][j] = math.log10(mat[i][j])
# print(mat)

mean_vec = []
center_mat = mat[:]
for col in range(len(mat[0])):
  sum = 0
  for row in range(len(mat)):
    sum += mat[row][col]
  mean_vec.append(sum / len(mat))

for row in range(len(mat)):
  for col in range(len(mat[0])):
    center_mat[row][col] = center_mat[row][col] - mean_vec[col]

print(f'data after log and center process: \n{center_mat}')

"""c. Perform PCA on the data. Remember to center the data points first by computing the mean data vector μ and subtracting it from every point. With the centered data matrix, do an SVD and compute the principal components.

pca with svd
"""

svd = PCA(n_components=2, svd_solver='auto')
new_data_svd = svd.fit_transform(center_mat)
# print(new_data_svd)
plt.scatter(new_data_svd[:, 0], new_data_svd[:, 1])
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.title('Projection on two principal components')
plt.show()

"""d. Write down the first two principal components v1 and v2. Provide a qualitative interpreta- tion of the components. Which among the nine factors do they appear to correlate the most with?"""

pca = PCA(n_components=2)
new_data = pca.fit_transform(center_mat)
pca_components = pca.components_.T

print(pca.explained_variance_ratio_)
print(pca_components)
print('\nArts most correlated to v1 and Hlthcare is most correlated to v2')

plt.scatter(new_data[:, 0], new_data[:, 1])
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.title('Projection on two principal components')
plt.show()

"""e. Project the data points onto the first two principal components. (That is, compute the highest 2 scores of each of the data points.) Plot the scores as a 2D scatter plot. Which cities correspond to outliers in this scatter plot?"""

# get city names
cities = []
city_name = data_set.index.values
for i in range(len(city_name)):
  cities.append(city_name[i])

pca_res = center_mat.dot(pca_components)
# print(pca_res)
city_res = []
for i, j in enumerate(pca_res):
  r1, r2 = j[0], j[1]
  # according to graph above, we can set up bounds to get strange pointers
  if r1 > 1.5 or r2 < -0.6 or r2 > 0.7:
    city_res.append(cities[i])
print(f'the outlier cities are \n{city_res}')

"""f. Repeat Steps 2-5, but with a slightly different data matrix – instead of computing the
base-10 logarithm, use the z-scores. (The z-score is calculated by computing the mean μ and standard deviation σ for each feature, and normalizing each entry x by x−μ ). How do your answers change?
"""

pca_z_score = PCA(n_components=2, svd_solver='auto')
pca_z_score.fit_transform(data_z_score)
pca_z_score_components = pca_z_score.components_.T
pca_z_score_res = data_z_score.dot(pca_z_score_components)
# print(len(pca_z_score_res))

plt.scatter(pca_z_score_res[:, 0], pca_z_score_res[:, 1])
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.title('Projection on two principal components')
plt.show()

"""inorder to get strange points in the graph, we can set bound to pc1 > 6 or pc2 < -2.5"""

city_res_z_score = []
for i, j in enumerate(pca_z_score_res):
  r1, r2 = j[0], j[1]
  # according to graph above, we can set up bounds to get strange pointers
  if r1 > 6 or r2 < -2.5:
    city_res_z_score.append(cities[i])
print(f'the outlier cities are \n{city_res_z_score}')

"""Analyse two methods"""

print(list(set(city_res).intersection(set(city_res_z_score))))
print(list(set(city_res).difference(set(city_res_z_score))))

"""we can see that the city New-york shows up in both methods"""