<jupyter_start><jupyter_text>3b_FactorialSearch notebook, Vehciels clustering analysisProject: clustering-analysis-domain-agnostic-features-2018Authors: <NAME>, <NAME>License: BSD 3-ClauseCopyright (c) 2021 Alliance for Sustainable Energy LLC DescriptionThis notebook uses Spark to parallelize a hyperparameter search over the clustering analysis. It produces the violin plots in Figure 7 and the data in Table 4.This notbook must be run with a spark context. The easiet way to do this is to boot up Spark using start_spark_jupyter_notebook.sh<jupyter_code>%%time
%reset -f

%load_ext autoreload
%autoreload 2
%matplotlib inline

import itertools
import pyspark
from pyspark.sql import SparkSession, Row
import json
import pandas as pd
import numpy as np
from sklearn.cluster import KMeans, DBSCAN
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE
import math, random
import pickle
from sklearn.metrics import silhouette_score
from matplotlib import pyplot as plt


sc = SparkSession.builder.getOrCreate().sparkContext

# Load feature dataframes from CSV
features_traditional = pd.read_csv("./data/FleetDNAETL_CoDA_epaprime_traditional_nolimit.csv").set_index("vdir").sort_index()
features_agnostic = pd.read_csv("./data/FleetDNAETL_CoDA_epaprime_agnostic_50klimit.csv").set_index("vdir").sort_index()
cols = [x for x in features_agnostic.columns if x[0]!='g']
features_agnostic = features_agnostic[cols]

# Load metadata from CSV
meta = pd.DataFrame.from_csv("./data/vehicle_specifications_fdna.csv")
m = meta.reset_index()
m["vdir"] = "v_" + m["id"].astype(str)
m = m.set_index("vdir").sort_index()


X_agnostic = features_agnostic
X_traditional = features_traditional

# Preprocessing

X_agnostic = X_agnostic.dropna(axis=1)
X_traditional = X_traditional.dropna(axis=1)

def normalize(df):
    return (df-df.min())/(df.max()-df.min())

X_agnostic = normalize(features_agnostic).dropna(axis=1)
X_traditional = normalize(features_traditional).dropna(axis=1)

X_agnostic = X_agnostic.reindex_axis(sorted(X_agnostic.columns), axis=1)

corr_matrix = X_agnostic.corr().abs()
upper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(np.bool))
to_drop = [column for column in upper.columns if any(upper[column] > 0.90)]
X_agnostic = X_agnostic.drop(to_drop, axis=1)

corr_matrix = X_traditional.corr().abs()
upper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(np.bool))
to_drop = [column for column in upper.columns if any(upper[column] > 0.90)]
X_traditional = X_traditional.drop(to_drop, axis=1)

X_traditional = normalize(X_traditional)
X_agnostic = normalize(X_agnostic)


def clusters(alg, data, index):
    C = alg.fit_predict(data)
    return pd.DataFrame(C).set_index(index)

def do_projection(p,f):
    proj = p.fit_transform(f)
    proj = pd.DataFrame(proj).set_index(f.index)
    proj = normalize(proj) # This helps guard against bias due to scale
    return proj


### Define VI

def variation_of_information_direct(X, Y):
    """
    Compute the variation of information between two cluster labelings
    Should be bounded by 0 <= vi <= ln(n)
    """
    dict_X = {}
    for x in set(X):
        dict_X[x] = []
    for i in range(len(X)):
        dict_X[X[i]].append(i)
    dict_Y = {}
    for y in set(Y):
        dict_Y[y] = []
    for i in range(len(Y)):
        dict_Y[Y[i]].append(i)

    n = float(len(X))

    VI = 0
    for i in dict_X.keys():
        for j in dict_Y.keys():
            Xi = set(dict_X[i])
            Yj = set(dict_Y[j])
            pi = len(Xi) / n
            qj = len(Yj) / n
            rij = len(Xi.intersection(Yj)) / n
            if rij != 0:
                VI = VI + rij * (math.log(rij / pi, 2) + math.log(rij / qj, 2))

    return -VI

def vi(left, right):
    comp = left.merge(right, how='inner', right_index=True, left_index=True)
    return variation_of_information_direct(np.array(comp["0_x"]), np.array(comp["0_y"]))


# Set up the parallel computation

#_bcast_ground_truth = sc.broadcast(ground_truth)
_bcast_features_agnostic = sc.broadcast(X_agnostic)
_bcast_features_traditional = sc.broadcast(X_traditional)
_bcast_metadata = sc.broadcast(m)


def unpack_projection(dict_projection):
    if dict_projection[u"type"] == u"PCA":
        return PCA(n_components=dict_projection[u"n_components"])
    elif dict_projection[u"type"] == u"TSNE":
        return TSNE(perplexity=dict_projection[u"perplexity"])
    else:
        return dict_projection

def unpack_clusteting(dict_clustering):
    if dict_clustering[u"type"] == u"KMeans":
        return KMeans(n_clusters=dict_clustering[u"n"])
    elif dict_clustering[u"type"] == u"DBSCAN":
        return DBSCAN(eps=dict_clustering[u"eps"])
    else:
        return dict_clustering

def unpack_features(dict_features):
    if dict_features == u"traditional":
        return _bcast_features_traditional.value
    elif dict_features == u"agnostic":
        return _bcast_features_agnostic.value
    else:
        raise Exception("Features value was {}".format(dict_features))
    
def unpack_labels(cd):
    return cd["features"], cd["projection"][u"type"], cd["clustering"][u"type"]

def unpack(cd):
    features = unpack_features(cd["features"])
    projection = unpack_projection(cd["projection"])
    clustering = unpack_clusteting(cd["clustering"])
    return features, projection, clustering

def constraints(c):
    if (c[0].value_counts() < 10).any():
        return "some cluster has less than 10"
    if len(set(c[0].to_dict().values())) < 3:
        return "not at least three clusters"
    else:
        return False

def do_computation(computation_descriptor):
    # Unpack argument for this computation from Pyspark cell
    cd = json.loads(computation_descriptor)
    features, projection, clustering = unpack(cd)
    features_label, projection_label, clustering_label = unpack_labels(cd)
    
    # Execute the computation
    p = do_projection(projection, features)
    c = clusters(clustering, p, features.index)
    
    # Enforce constraints
    if constraints(c):
        return None
    
    # Compute metrics
    _silhouette_score = silhouette_score(p, c)
    
    _vi = {}
    m = _bcast_metadata.value
    ground_truth_cols = ['class', 'vocation', 'vehicletype', 'providers']
    c.columns = map(str, c.columns)
    merged = c.merge(m[ground_truth_cols], how='inner', right_index=True, left_index=True)
    for col_name in ground_truth_cols:
        ourMerge = merged[["0", col_name]].dropna()
        vvi = variation_of_information_direct(np.array(ourMerge["0"]), np.array(ourMerge[col_name]))
        _vi["vi_{}".format(col_name)] = vvi
    
    
    ans = {"features":features_label,
            "projection":projection_label,
            "clustering":clustering_label,
            "silhouette_score":_silhouette_score}
    
    ans.update(_vi)
    
    return ans

# Create a collection of experiment configuration parameters
random.seed(2018)
feature_opts = ["traditional", "agnostic"]

clustering_opts_dbscan = [{"type":"DBSCAN", "eps":n} for n in np.linspace(0.01,.1,10)]
clustering_opts_kmeans = [{"type":"KMeans", "n":n} for n in range(2,5)]
clusetring_opts = itertools.chain(clustering_opts_dbscan, clustering_opts_kmeans)

projection_opts_PCA = [{"type":"PCA", "n_components":2}]
projection_opts_TSNE = [{"type":"TSNE", "perplexity":p} for p in np.linspace(5,50,10)]
projection_opts = itertools.chain(projection_opts_PCA, projection_opts_TSNE)

product_tuple = itertools.product(clusetring_opts, projection_opts, feature_opts)
configuration_dict = [{"clustering":_c0, "projection":_c1, "features":_c2} for _c0, _c1, _c2 in product_tuple]
configuration_json = map(json.dumps, configuration_dict)

# Run map over parallelized rdd
configuration_rdd = sc.parallelize(configuration_json, len(configuration_json))
spark_result = configuration_rdd.map(do_computation)<jupyter_output>CPU times: user 2.3 s, sys: 2.04 s, total: 4.34 s
Wall time: 2.76 s<jupyter_text>The following cell will execute a potentially time consuming spark job<jupyter_code>len(configuration_json)
answer = spark_result.collect()
df = pd.DataFrame([a for a in answer if a])<jupyter_output><empty_output><jupyter_text>Factor Analysis<jupyter_code># Table of factor sensitivities
metrics = ["silhouette_score", "vi_class", "vi_vehicletype", "vi_vocation"]

# Clustering
c = df.groupby("clustering").agg(['mean', 'std', 'min'])
f = df.groupby("features").agg(['mean', 'std', 'min'])
p = df.groupby("projection").agg(['mean', 'std', 'min'])


p["type"] = "projection"
c["type"] = "clustering"
f["type"] = "features"

final_df = p.append(c).append(f)
final_df.index.name='subtype'
final_df.set_index('type', append=True, inplace=True)
final_df = final_df.reorder_levels(['type', 'subtype'])
final_df.round(2)
f = open("paper_figs/latex_table.txt", "w")
f.write(final_df.to_latex())
f.close()<jupyter_output><empty_output><jupyter_text>With Silhouette Score<jupyter_code>import seaborn as sns
a = df[["clustering", "silhouette_score"]].rename(columns={"clustering":"bin"})
b = df[["features", "silhouette_score"]].rename(columns={"features":"bin"})
c = df[["projection", "silhouette_score"]].rename(columns={"projection":"bin"})

df_long = a.append(b).append(c)

sns.boxplot(x="bin", y="silhouette_score", data=df_long[df_long["silhouette_score"]>0.1]);
a = df[["clustering", "silhouette_score", "features"]].rename(columns={"clustering":"bin"})
b = df[["projection", "silhouette_score", "features"]].rename(columns={"projection":"bin"})

df_long_split = a.append(b)

fig = plt.figure(figsize=(5,4))

h = sns.violinplot(x="bin", y="silhouette_score", hue="features", data=df_long_split,
               split=True, inner="stick", palette="Set3");
h.set(xlabel='Algorithm', ylabel='Silhouette Score')

df_long_split.groupby(["bin","features"]).agg(['mean', 'std'])

fig.savefig("./paper_figs/violin_plot_silhouette_score.eps")
a = df[["clustering", "silhouette_score", "features"]].rename(columns={"clustering":"bin"})
b = df[["projection", "silhouette_score", "features"]].rename(columns={"projection":"bin"})

df_long_split_ss = a.append(b)


sns.violinplot(x="features", y="silhouette_score", data=df_long_split,
               split=True, inner="stick", palette="Set3");

df_long_split_ss.groupby(["features","bin"]).agg(['mean', 'std'])<jupyter_output><empty_output><jupyter_text>With VI<jupyter_code>a = df[["clustering", "vi_class", "vi_vocation", "vi_vehicletype", "features"]].rename(columns={"clustering":"bin"})
b = df[["projection", "vi_class", "vi_vocation", "vi_vehicletype", "features"]].rename(columns={"projection":"bin"})

df_long_split = a.append(b)


fig = plt.figure(figsize=(5,4))
h = sns.violinplot(x="bin", y="vi_vehicletype", hue="features", data=df_long_split,
               split=True, inner="stick", palette="Set3")
h.set(xlabel='Algorithm', ylabel='Variation of Information to Vehicle Type')
fig.savefig("./paper_figs/violin_plot_vi_type.eps")

fig = plt.figure(figsize=(5,4))
h = sns.violinplot(x="bin", y="vi_vocation", hue="features", data=df_long_split,
               split=True, inner="stick", palette="Set3")
h.set(xlabel='Algorithm', ylabel='Variation of Information to Vehicle Vocation')
fig.savefig("./paper_figs/violin_plot_vi_vocation.eps")

fig = plt.figure(figsize=(5,4))
h = sns.violinplot(x="bin", y="vi_class", hue="features", data=df_long_split,
               split=True, inner="stick", palette="Set3");
h.set(xlabel='Algorithm', ylabel='Variation of Information to Vehicle Class')
fig.savefig("./paper_figs/violin_plot_vi_class.eps")

df_long_split.groupby(["bin","features"]).agg(['mean', 'std']).round(2)<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup