# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %matplotlib inline
import re
import json
import mca
import pandas as pd
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
from ast import literal_eval
from scipy.spatial.distance import cosine

# http://www.gastonsanchez.com/visually-enforced/how-to/2012/10/13/MCA-in-R/  
# https://pypi.org/project/mca/

datadir = '/Users/klarnemann/Documents/Insight/Insight_project/data'
figdir = '/Users/klarnemann/Documents/Insight/Insight_project/docs/figures'

poi_category_df_f = '%s/poi_museums_and_parks.xlsx' % (datadir)
poi_category_df = pd.read_excel(poi_category_df_f, index_col=0)

category_df_f = '%s/poi_by_city_museums_and_parks.xlsx' % (datadir)
category_df = pd.read_excel(category_df_f, index_col=0)

poi_category_df.head()

category_df.head()


# ## Functions

def compute_pairwise_cosine_distance(data):
    n_cities, _ = data.shape
    out_mat = np.zeros(shape=(n_cities, n_cities))
    for i in np.arange(n_cities):
        if i % 200 == 0:
            print(i)
        for j in np.arange(n_cities):
            if i < j:
                out_mat[i,j] = cosine(data[i,:], data[j,:])
    return out_mat


# # K-Means Clustering
#
# https://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html#sklearn.cluster.KMeans
#
# https://scikit-learn.org/stable/modules/generated/sklearn.metrics.silhouette_score.html#sklearn.metrics.silhouette_score
#
# https://www.datanovia.com/en/lessons/determining-the-optimal-number-of-clusters-3-must-know-methods/#elbow-method

from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_samples, silhouette_score

data = np.array(category_df[category_df.columns[3:]])

kmeans_stats = []

for k in np.linspace(5,1000,200, dtype=int):
    kmeans = KMeans(n_clusters=k, random_state=0).fit(data)
    score = silhouette_avg = silhouette_score(data, kmeans.labels_)
    inertia = kmeans.inertia_
    kmeans_stats.append([k, score, inertia])

#kmeans_silhouette = np.array(kmeans_silhouette)
#kmeans_silhouette.shape
kmeans_stats = np.array(kmeans_stats)
kmeans_stats.shape

plt.figure(figsize=(6,4))
plt.plot(kmeans_stats[:,0], kmeans_stats[:,2])
plt.plot([100,100], [0,25000], 'r--')
plt.xlabel('# Clusters', fontsize=14)
plt.ylabel('Within-Cluster\nSum of Squares', fontsize=14)
plt.tight_layout()
plt.savefig('%s/plot_elbow_kmeans_categorical_by_city.png' % (figdir))

kmeans_100 = KMeans(n_clusters=100, random_state=0).fit(data)
kmeans_labels = kmeans_100.labels_

# # Multiple Correspondence Analysis

print(mca.MCA.__doc__)

category_df.columns[3:]

mca_ben = mca.MCA(category_df[category_df.columns[3:]])

fs_cols = mca_ben.fs_c()
print(fs_cols.shape)

fs_cols

fs_rows = mca_ben.fs_r()
fs_rows.shape

fs_cosine_mat = compute_pairwise_cosine_distance(np.array(fs_rows))

fs_cosine_mat[np.isnan(fs_cosine_mat)] = 0.

expl_var_by_row = mca_ben.expl_var()

n_factors= len(expl_var_by_row)

plt.figure(figsize=(8,4))
plt.bar(np.arange(n_factors)[:30], expl_var_by_row[:30]);
plt.xlim(-0.75, 29.75)
plt.xticks(np.arange(n_factors)[:30], np.arange(n_factors)[1:31], rotation=90);
plt.xlabel('Component', fontsize=14)
plt.ylabel('% Variance Explained', fontsize=14);
plt.tight_layout()
plt.savefig('%s/bar_mca_explained_variance_factors.png' % (figdir), dpi=150, transparent=True)

# ~~~
# for i, col in enumerate(poi_by_city.columns[2:]):
#     tmp_data = list(fs_cols[i,:])
#     max_val = np.max(tmp_data)
#     max_col_1 = tmp_data.index(max_val)
#     tmp_data[max_col_1] = 0
#     max_val = np.max(tmp_data)
#     max_col_2 = tmp_data.index(max_val)
#     print(col, max_col_1, max_col_2)
# ~~~

museum_cats = category_df.columns[2:27]
park_cats = category_df.columns[27:]
museums = category_df[museum_cats].sum(axis=1)
parks = category_df[park_cats].sum(axis=1)

n_cities = float(len(category_df.index))

# +
vals = category_df[museum_cats].sum(axis=0).sort_values(ascending=True)/n_cities*100
xlabels = list(category_df[museum_cats].sum(axis=0).sort_values(ascending=True).index)

plt.figure(figsize=(6,7))
plt.barh(np.arange(len(xlabels)), vals);
plt.xlim(0,50);
plt.ylim(-0.75,len(xlabels)-0.5)
plt.yticks(np.arange(len(xlabels)), xlabels, fontsize=14);
plt.xlabel('% of Cities', fontsize=14)
plt.title('Museum Features', fontsize=14)
plt.tight_layout()
plt.savefig('%s/bar_pct_museum_features_by_city.png' % (figdir), dpi=150, transparent=True)

# +
vals = category_df[park_cats].sum(axis=0).sort_values(ascending=True)/n_cities*100
xlabels = list(category_df[park_cats].sum(axis=0).sort_values(ascending=True).index)

plt.figure(figsize=(6,7))
plt.barh(np.arange(len(xlabels)), vals);
plt.xlim(0,30)
plt.ylim(-0.75,len(xlabels)-0.5)
plt.yticks(np.arange(len(xlabels)), xlabels, fontsize=14);
plt.xlabel('% of Cities', fontsize=14)
plt.title('Park Features', fontsize=14)
plt.tight_layout()
plt.savefig('%s/bar_pct_park_features_by_city.png' % (figdir), dpi=150, transparent=True)
# -

category_df.columns

plt.figure(figsize=(5,5))
plt.scatter(fs_data[:,0], fs_data[:,1], c=np.array(parks.values>0), 
            cmap='seismic', alpha=0.75)
plt.title('Presence vs. Absence\nof Park Features', fontsize=14)
plt.xlabel('Component 1\n(from MCA)', fontsize=14)
plt.ylabel('Component 2\n(from MCA)', fontsize=14)
plt.tight_layout()
plt.savefig('%s/scatter_mca_componenet_1_vs_2_parks.png' % (figdir))

plt.figure(figsize=(5,5))
plt.scatter(fs_data[:,0], fs_data[:,1], c=np.array(category_df['Local history']>0), 
            cmap='seismic', alpha=0.75)
plt.title('Presence vs. Absence\nof Local History Museum Feature', fontsize=14)
plt.xlabel('Component 1\n(from MCA)', fontsize=14)
plt.ylabel('Component 2\n(from MCA)', fontsize=14)
plt.tight_layout()
plt.savefig('%s/scatter_mca_componenet_1_vs_2_local_history.png' % (figdir))

# ## K-Means

fs_data = fs_rows

fs_kmeans_stats = []

for k in np.linspace(5,1000,200, dtype=int):
    kmeans = KMeans(n_clusters=k, random_state=0).fit(fs_data)
    score = silhouette_avg = silhouette_score(fs_data, kmeans.labels_)
    inertia = kmeans.inertia_
    fs_kmeans_stats.append([k, score, inertia])

plt.figure(figsize=(6,4))
plt.plot(fs_kmeans_stats[:,0], fs_kmeans_stats[:,2])
plt.plot([90,90], [0,25000], 'r--')
plt.xlabel('# Clusters', fontsize=14)
plt.ylabel('Within-Cluster\nSum of Squares', fontsize=14)
plt.tight_layout()
plt.savefig('%s/plot_elbow_kmeans_categorical_by_city_fs.png' % (figdir))

score_results = []

for k in np.linspace(5,1000,200, dtype=int):
    print(k)
    kmeans = KMeans(n_clusters=k, random_state=0).fit(fs_data)
    fs_score = silhouette_avg = silhouette_score(fs_data, kmeans.labels_)
    score = silhouette_avg = silhouette_score(data, kmeans.labels_)
    score_results.append([fs_score, score])

score_results = np.array(score_results)

plt.figure(figsize=(6,4))
k_values = np.linspace(5,200,40, dtype=int)
plt.plot(k_values, score_results[:40,0], 'k-', label='Post-MCA')
plt.plot(k_values, score_results[:40,1], 'r-', label='Pre-MCA')
plt.plot([25,25],[0,0.8],'k--')
plt.xlabel('# Clusters', fontsize=14)
plt.ylabel('Silhouette Score', fontsize=14)
plt.legend()
plt.tight_layout()
plt.savefig('%s/plot_silhouette_mca_vs_original_by_city_fs.png' % (figdir))

# choose 25 cluster solution using the elbow convention
fs_kmeans = KMeans(n_clusters=25, random_state=0).fit(fs_data)
fs_kmeans_labels = fs_kmeans.labels_

cluster_dict = dict(zip(list(category_df.index), fs_kmeans_labels.astype(str)))
cluster_dict_f = '%s/results/cluster_dict_kmeans_25_mca_21_by_city.json' % (datadir)
with open(cluster_dict_f, 'w') as fp:
    json.dump(cluster_dict, fp, sort_keys=True, indent=4)
#with open(cluster_dict_f, 'r') as fp:
#    cluster_dict = json.load(fp)
#cluster_dict

plt.figure(figsize=(7,5))
fs_corr = pd.DataFrame(fs_data[reorder]).T.corr()
plt.imshow(fs_corr)
plt.clim(-1,1)
cbar = plt.colorbar()
cbar.set_label('Correlation of Features (MCA space)')
plt.xlabel('Cities', fontsize=14)
plt.ylabel('Cities', fontsize=14)
plt.tight_layout()
plt.savefig('%s/corr_mat_fs-mca_fs.png' % (figdir))

plt.figure(figsize=(7,5))
fs_corr = pd.DataFrame(data[reorder]).T.corr()
plt.imshow(fs_corr)
plt.clim(-1,1)
cbar = plt.colorbar()
cbar.set_label('Correlation of Features (MCA space)')
plt.xlabel('Cities', fontsize=14)
plt.ylabel('Cities', fontsize=14)
plt.tight_layout()
plt.savefig('%s/corr_mat_fs-mca.png' % (figdir))

silhouette_score(fs_data, fs_kmeans_labels)

silhouette_score(data, fs_kmeans_labels)

n,_ = fs_cosine_mat.shape
reorder = [x for _,x in sorted(zip(fs_kmeans_labels, np.arange(n)))]
fs_cosine_mat_clusters = fs_cosine_mat.copy()
fs_cosine_mat_clusters = fs_cosine_mat_clusters[reorder,:]
fs_cosine_mat_clusters = fs_cosine_mat_clusters[:,reorder]
plt.imshow(fs_cosine_mat_clusters, interpolation='nearest')
cbar = plt.colorbar()
cbar.set_label('Cosine Distance/n(in MCA space)')

# ### Ward Clustering
#
# https://scikit-learn.org/stable/modules/generated/sklearn.cluster.AgglomerativeClustering.html#sklearn.cluster.AgglomerativeClustering
#
# https://stackabuse.com/hierarchical-clustering-with-python-and-scikit-learn/



# ## Latent Corresponce Analysis

StarcoderdataJupyterScriptsDedupFiltered