<file_sep># Import packages
import tensorflow as tf
import pandas as pd
import numpy as np
import os
os.chdir('C://Users//vnino//Desktop//CE888 - Assignment 2//Python code//mushrooms//data')
import random

# Fetch the dataset
names=['label', 'cap_shape', 'cap_surface', 
       'cap_color', 'bruises', 'odor', 
       'gill_attachment', 'gill_spacing', 'gill_size',
       'gill_color', 'stalk_shape', 'stalk_root', 
       'stalk_surface_above', 'stalk_surface_below', 
       'stalk_color_above', 'stalk_color_below', 'veil_type',
       'veil_color', 'ring_number', 'ring_type',
       'spore_print_color', 'population','habitat']
df = pd.read_csv("agaricus-lepiota.csv", names=names)

# Filter missing values
df=df[df.stalk_root!='?']

df=df.drop(columns=['veil_type', 'veil_color'])

data_df=pd.get_dummies(df, prefix=['label', 'cap_shape', 'cap_surface', 
       'cap_color', 'bruises', 'odor', 
       'gill_attachment', 'gill_spacing', 'gill_size',
       'gill_color', 'stalk_shape', 'stalk_root', 
       'stalk_surface_above', 'stalk_surface_below', 
       'stalk_color_above', 'stalk_color_below',
       'ring_number', 'ring_type',
       'spore_print_color', 'population','habitat'])

data_df=data_df.drop(columns=['label_e'])
data_df=data_df.rename(index=str, columns={'label_p': 'label'})

#Defining labels and features
data_df.reset_index(inplace=True, drop=True)
target=data_df['label']
data=np.asarray(data_df.iloc[:,21:])
n_samples= data.shape[0]
n_clusters=2

# Get the split between training/test set and validation set
random.seed(0)
test_indices = random.sample(range(0, 5644), 1120)
validation_indices = [x for x in range(0, 5644) if x not in test_indices]

# Auto-encoder architecture
input_size = data.shape[1]
hidden_1_size = 500
hidden_2_size = 500
hidden_3_size = 2000
embedding_size = n_clusters
dimensions = [hidden_1_size, hidden_2_size, hidden_3_size, embedding_size, # Encoder layer dimensions
              hidden_3_size, hidden_2_size, hidden_1_size, input_size] # Decoder layer dimensions
activations = [tf.nn.relu, tf.nn.relu, tf.nn.relu, None, # Encoder layer activations # None is linear activation
               tf.nn.relu, tf.nn.relu, tf.nn.relu, None] # Decoder layer activations # None is linear activation
names = ['enc_hidden_1', 'enc_hidden_2', 'enc_hidden_3', 'embedding', # Encoder layer names
         'dec_hidden_1', 'dec_hidden_2', 'dec_hidden_3', 'output'] # Decoder layer names

activations_soft = [tf.nn.relu, tf.nn.relu, tf.nn.relu, tf.nn.softmax, # Encoder layer activations # Uses softmax in last enocder layer
               tf.nn.relu, tf.nn.relu, tf.nn.relu, None]


#data_df['label'].value_counts()<file_sep># Import packages
import tensorflow as tf
import pandas as pd
import numpy as np
import os
os.chdir('C://Users//vnino//Desktop//CE888 - Assignment 2//Python code//HAR//data')
import random

# Fetch the dataset
train_df= pd.read_csv("train.csv")
test_df= pd.read_csv("test.csv")
data_df=pd.concat([train_df,test_df], ignore_index=True)

#Defining labels and features
target=data_df['Activity']
data=np.asarray(data_df.iloc[:,0:560])
n_samples= data.shape[0]
n_clusters=6

# Get the split between training/test set and validation set
random.seed(0)
test_indices = random.sample(range(0, 10299), 2060)
validation_indices = [x for x in range(0, 10299) if x not in test_indices]

# Auto-encoder architecture
input_size = data.shape[1]
hidden_1_size = 500
hidden_2_size = 500
hidden_3_size = 2000
embedding_size = n_clusters
dimensions = [hidden_1_size, hidden_2_size, hidden_3_size, embedding_size, # Encoder layer dimensions
              hidden_3_size, hidden_2_size, hidden_1_size, input_size] # Decoder layer dimensions
activations = [tf.nn.relu, tf.nn.relu, tf.nn.relu, None, # Encoder layer activations # None uses linear act.
               tf.nn.relu, tf.nn.relu, tf.nn.relu, None] # Decoder layer activations # None uses linear act.
names = ['enc_hidden_1', 'enc_hidden_2', 'enc_hidden_3', 'embedding', # Encoder layer names
         'dec_hidden_1', 'dec_hidden_2', 'dec_hidden_3', 'output'] # Decoder layer names

activations_soft = [tf.nn.relu, tf.nn.relu, tf.nn.relu, tf.nn.softmax, # Encoder layer activations # Uses softmax in last enocder layer
               tf.nn.relu, tf.nn.relu, tf.nn.relu, None]

#data_df['ActivityName'].value_counts()
<file_sep>INSTRUCTIONS FOR RUNNING THE CLUSTERING
-----------------------------------------------------------------------

1. SETTING PATH FOR INPUT DATA: 

The input data sets are stored in the input folders: 
- digits
- HAR
- mushrooms

Before running the algorithms, you will need to edit the specs files 
(HAR_specs.py, mushrooms_specs.py) inside these folders to change the
path of the text files. The text files are stored in the same folders 
as the specs files.  

2. RUNNING CLUTERING:

Each clustering technique is programmed in a separate .py file: 
- kmeans.py
- autoencoder.py
- autoencodersoftmax.py
- deepkmeans.py

You can run each of these programs on your command windows or using 
a python IDE. For each run the program will ask to insert the path 
where the specs files are stored and to specify the dataset you want 
to process.

3. OUTPUT: 

The run has 2 files as output written in the folder "results" that can 
be found inside of each the input folders. The files are: a .txt file 
containing the results (performance scores) of the run and a .png file 
containing the confusion matrix for the clusters obtained.

The folders already contain the results of running the codes. Re-running
it, will over write the files.

ADDITIONAL FILES
-----------------------------------------------------------------------

Inside the input folders, you can also found the code used in assignment 1
to make the previous EDA of the three data sets.
<file_sep>import pandas as pd
import numpy as np

def power(sample1, sample2, reps, perm_iterations, alpha):
    tobs=np.mean(sample1)-np.mean(sample2)
    m=0
    for r in range(reps):
        n=0
        for i in range(perm_iterations):
            join_sample=np.append(sample1, sample2)
            resample1=np.random.choice(join_sample, size=sample1.shape[0], replace=True)
            resample2=np.random.choice(join_sample, size=sample2.shape[0], replace=True)
            if (np.mean(resample1)-np.mean(resample2)) > tobs:
                n+=1
        p_value=n/perm_iterations
        print(p_value)
        if p_value<alpha:
            m+=1
    return(m/reps)

sample2=np.arange(0,1000,1)
sample1=np.arange(0,1000,1)

power(sample1,sample2,100,100,0.05)         
    
np.mean(sample1)-np.mean(sample2)  

join_sample=np.append(sample1, sample2)
resample1=np.random.choice(join_sample, size=sample1.shape[0], replace=True)
resample2=np.random.choice(join_sample, size=sample2.shape[0], replace=True)

np.mean(resample1)-np.mean(resample2)

    
<file_sep># LAB 2

## **Exercise 1 - Plotting**

The code for this exercise can be found in the folder Exercise 1 - Plotting, in the file vehicles.py.

First, we create an histogram on the Current Fleet and compute the main descriptive statistics:

![logo](./HistogramCurrentFleet.png?raw=true)

Mean: 20.144578
Median: 19.000000
Var: 40.983113
std: 6.401805
MAD: 4.000000

Second, we create an histogram on the New Fleet and compute the main descriptive statistics:

![logo](./HistogramNewFleet.png?raw=true)

Mean: 30.481013
Median: 32.000000
Var: 36.831918
std: 6.068931
MAD: 4.000000

A scatterplot is also computed:

![logo](./scaterplot.png?raw=true)

## **Exercise 2 - Bootstrap**

The code for this exercise can be found in the folder Exercise 2 - Bootstrap, in the file bootstrap.py.

First, the function boostrap() is created. 
Then this is used to compare the mean values of the MPG for each Fleet at a 95% CI.

The results of this analysis are: 

Current Fleet Mean: 20.161739
Current Fleet Lower: 19.321185
Current Fleet Upper: 20.960040

New Fleet Mean: 30.481823
New Fleet Lower: 29.176899
New Fleet Upper: 31.887025
<file_sep># For creating this code the following sources were used: 
# 1. Building autoencoders in Keras. Source: https://blog.keras.io/building-autoencoders-in-keras.html
# 2. Implementation of "Deep k-Means: Jointly Clustering with k-Means and Learning Representations" by <NAME>, <NAME>, and <NAME>. (https://github.com/MaziarMF/deep-k-means)
# 3. In Depth: k-Means Clustering (https://jakevdp.github.io/PythonDataScienceHandbook/05.11-k-means.html)

# IMPORTS PACKAGES ############################################################################################################################
###############################################################################################################################################

import os
import math
import numpy as np
import tensorflow as tf
import sklearn
from sklearn.metrics import  silhouette_score, accuracy_score, completeness_score, homogeneity_score
from sklearn.cluster import KMeans
from scipy.stats import mode
import seaborn as sns
sns.set(font_scale=3)
import matplotlib.pyplot as plt
TF_FLOAT_TYPE = tf.float32

from keras.optimizers import SGD, Adam
from keras.layers import Dense, Input
from keras.models import Model

from keras.backend import clear_session


# FUNCTIONS DEFINITION ########################################################################################################################
###############################################################################################################################################

def cluster_acc(y_true, y_pred, n_clusters):
    labels = np.zeros_like(y_pred)
    for i in range(n_clusters):
        mask = (y_pred == i)
        labels[mask] = mode(y_true[mask])[0]
    if y_true is not None:
        acc = np.round(accuracy_score(y_true, labels), 5)
    return acc

def save_confusion_matrix(y_true, y_pred, n_clusters, name):
    labels = np.zeros_like(y_pred)
    for i in range(n_clusters):
        mask = (y_pred == i)
        labels[mask] = mode(y_true[mask])[0]
    confusion_matrix = sklearn.metrics.confusion_matrix(y_true, labels)
    plt.figure(figsize=(16, 14))
    sns.heatmap(confusion_matrix, annot=True, fmt="d", annot_kws={"size": 20});
    plt.title("Confusion matrix", fontsize=30)
    plt.ylabel('True label', fontsize=25)
    plt.xlabel('Clustering label', fontsize=25)
    plt.savefig('./results/'+ name + '.png')

def next_batch(num, data):
    indices = np.arange(0, data.shape[0])
    np.random.shuffle(indices)
    indices = indices[:num]
    batch_data = np.asarray([data[i, :] for i in indices])
    return indices, batch_data

def create_autoencoder(dims):
    input_layer = Input(shape=(dims[0],), name='input')
    encoded = Dense(dims[1], activation='relu', name='enc_0')(input_layer)
    encoded = Dense(dims[2], activation='relu', name='enc_1')(encoded)
    encoded = Dense(dims[3], activation='relu', name='enc_2')(encoded)
    encoded = Dense(dims[4], activation='softmax', name='enc_3')(encoded) ##, activation='softmax'
    
    decoded = Dense(dims[3], activation='relu', name='dec_3')(encoded)
    decoded = Dense(dims[2], activation='relu', name='dec_2')(decoded)
    decoded = Dense(dims[1], activation='relu', name='dec_1')(decoded)
    decoded = Dense(dims[0], activation='sigmoid', name='dec_0')(decoded) # activation='sigmoid'
    
    return Model(input_layer, outputs=decoded, name='AE'), Model(input_layer, outputs=encoded, name='encoder')

# MAIN PROGRAM ################################################################################################################################
###############################################################################################################################################

# This program executes DEEP K-means for a chosen dataset. 
# It runs the algorithm 10 times, each with different seed.
# It computes the final perfomance metrics selecting the iteration with the 
# maximum Accuracy. The performance metrics are printed in the screen as 
# they are computed. Also, a txt file is generated with these results
# and 2 png file with the confusion matrices are created.

if __name__=='__main__':
    
# STEP 0: Setting the directory -------------------------------------------------------------------------------------------------------------
    directory = input('Please enter the directory where you have the datasets: ')
    os.chdir(directory)
    
# STEP 1 Importing dataset ------------------------------------------------------------------------------------------------------------------    
    dataset = input('Select the dataset for which you want to run DEEP K-Means, digits/HAR/mushrooms [d/h/m]: ')
    if dataset=='d':
        import digits_specs as specs 
    elif dataset=='h':
        import HAR_specs as specs 
    elif dataset=='m':
        import mushrooms_specs as specs 
    else: 
        print('Enter a valid dataset [d/h/m]: ')

# STEP 2: Running DEEP K-Means ------------------------------------------------------------------------------------------------------------------
    
    # 2.1. Setting initial parameters------------------------------------------------------------------------------------------------------------
    
  
    n_train_epochs = 50
    batch_size = 256 # Size of the mini-batches used in the stochastic optimizer
    seeded = True # Specify if runs are seeded
    seeds = [8905, 9129, 291, 4012, 1256, 6819, 4678, 6971, 1362, 575]  # Use a fixed seed for this run, as defined in the seed list 
                                                                        # There is randomness in centroids initialization and mini-batch selection
    n_runs = 10
    
    dims = [specs.input_size, 500, 500, 2000, specs.embedding_size]

    target = specs.target
    data = specs.data

    dict_acc = dict()
    dict_sc = dict()
    dict_cs = dict()
    dict_hs = dict()
    dict_clusterassign= dict()

    gpu_options = tf.GPUOptions(per_process_gpu_memory_fraction=0.05)
    config = tf.ConfigProto(gpu_options=gpu_options)

    # 2.2. Running k-means and creating the output files-----------------------------------------------------------------------------------------
    
    if os.path.exists('./results/AEsoft_results.txt'):
        os.remove('./results/AEsoft_results.txt') # Removes the text file if exists
    with open('./results/AEsoft_results.txt', 'x') as f: 
        
        # 2.2.1. Running k-means n_run times-------------------------------------------------------------------------------------------------------
       
        for run in range(n_runs):
            clear_session()
            if seeded:
                tf.reset_default_graph()
                tf.set_random_seed(seeds[run])
                np.random.seed(seeds[run])
            print("Run", run)
            f.write("RUN "+ str(run) +'\n'+'\n') 
            autoencoder, encoder = create_autoencoder(dims)
            autoencoder.compile(optimizer='adadelta', loss='binary_crossentropy')
            #autoencoder.compile(optimizer='adam', loss='kld')
            autoencoder.fit(data, data, batch_size=batch_size, epochs=n_train_epochs)
            # Second, run k-means++ on the trained embeddings
            print("Running k-means on the learned embeddings...")
            f.write("1. K-MEANS ON THE EMBEDDED FEATURES"+'\n'+'\n')
            #kmeans_model = KMeans(n_clusters=specs.n_clusters, init="k-means++").fit(embeddings) # n_init=10 initial centroids
            kmeans_model = KMeans(n_clusters=specs.n_clusters, init="k-means++").fit(encoder.predict(data))
            # Evaluate the clustering performance using the ground-truth labels
            sc = silhouette_score(data, kmeans_model.labels_)
            print("SC", sc)
            f.write("SC: "+ str(sc) +'\n')
        
            acc = cluster_acc(target, kmeans_model.labels_, specs.n_clusters)
            print("ACC", acc)
            f.write("ACC: "+ str(acc) +'\n')
            
            cs = completeness_score(target, kmeans_model.labels_)
            print("CS", cs)
            f.write("CS: "+ str(cs) +'\n')
            
            hs = homogeneity_score(target, kmeans_model.labels_)
            print("HS", hs)
            f.write("HS: "+ str(hs) +'\n')
                     
        # Record the clustering performance for the run
            dict_acc.update({run : acc})
            dict_sc.update({run : sc})
            dict_cs.update({run : cs})
            dict_hs.update({run : hs})
            dict_clusterassign.update({run : kmeans_model.labels_ })
            f.write('\n'+'-------------------------------------------------------------'+'\n') 
        
        max_acc = max(dict_acc, key=dict_acc.get)
        print('Iteration with maximum ACC: '+ str(max_acc))
        print("ACC: {:.3f}".format(dict_acc[max_acc]))
        print("SC: {:.3f}".format(dict_sc[max_acc]))
        print("CS: {:.3f}".format(dict_cs[max_acc]))
        print("HS: {:.3f}".format(dict_hs[max_acc]))
        f.write('Run with maximum ACC: '+ str(max_acc) +'\n')
        f.write("ACC: {:.3f}".format(dict_acc[max_acc]) +'\n')
        f.write("SC: {:.3f}".format(dict_sc[max_acc]) +'\n')
        f.write("CS: {:.3f}".format(dict_cs[max_acc]) +'\n')
        f.write("HS: {:.3f}".format(dict_hs[max_acc]) +'\n')
        save_confusion_matrix(target, dict_clusterassign[max_acc], specs.n_clusters, 'AEsoft_CM') 
   





<file_sep># For creating this code the following sources were used: 
# 1. Implementation of "Deep k-Means: Jointly Clustering with k-Means and Learning Representations" by <NAME>, <NAME>, and <NAME>. (https://github.com/MaziarMF/deep-k-means)
# 2. In Depth: k-Means Clustering (https://jakevdp.github.io/PythonDataScienceHandbook/05.11-k-means.html)

# IMPORTS PACKAGES ############################################################################################################################
###############################################################################################################################################

import os
import math
import numpy as np
import tensorflow as tf
import sklearn
from sklearn.metrics import  silhouette_score, accuracy_score, completeness_score, homogeneity_score
from sklearn.cluster import KMeans
from scipy.stats import mode
import seaborn as sns
sns.set(font_scale=3)
import matplotlib.pyplot as plt
TF_FLOAT_TYPE = tf.float32

# FUNCTIONS DEFINITION ########################################################################################################################
###############################################################################################################################################

def cluster_acc(y_true, y_pred, n_clusters):
    labels = np.zeros_like(y_pred)
    for i in range(n_clusters):
        mask = (y_pred == i)
        labels[mask] = mode(y_true[mask])[0]
    if y_true is not None:
        acc = np.round(accuracy_score(y_true, labels), 5)
    return acc

def save_confusion_matrix(y_true, y_pred, n_clusters, name):
    labels = np.zeros_like(y_pred)
    for i in range(n_clusters):
        mask = (y_pred == i)
        labels[mask] = mode(y_true[mask])[0]
    confusion_matrix = sklearn.metrics.confusion_matrix(y_true, labels)
    plt.figure(figsize=(16, 14))
    sns.heatmap(confusion_matrix, annot=True, fmt="d", annot_kws={"size": 20});
    plt.title("Confusion matrix", fontsize=30)
    plt.ylabel('True label', fontsize=25)
    plt.xlabel('Clustering label', fontsize=25)
    plt.savefig('./results/'+ name + '.png')
    
def next_batch(num, data):
    indices = np.arange(0, data.shape[0])
    np.random.shuffle(indices)
    indices = indices[:num]
    batch_data = np.asarray([data[i, :] for i in indices])
    return indices, batch_data

def fc_layers(input, specs):
    [dimensions, activations, names] = specs
    for dimension, activation, name in zip(dimensions, activations, names):
        input = tf.layers.dense(inputs=input, units=dimension, activation=activation, name=name, reuse=False)
    return input

def autoencoder(input, specs):
    [dimensions, activations, names] = specs
    mid_ind = int(len(dimensions)/2)
    # Encoder
    embedding = fc_layers(input, [dimensions[:mid_ind], activations[:mid_ind], names[:mid_ind]])
    # Decoder
    output = fc_layers(embedding, [dimensions[mid_ind:], activations[mid_ind:], names[mid_ind:]])
    return embedding, output

def f_func(x, y):
    return tf.reduce_sum(tf.square(x - y), axis=1)

def g_func(x, y):
    return tf.reduce_sum(tf.square(x - y), axis=1)

# COMPUTATION GRAPH FOR DEEP K-MEANS ##########################################################################################################
###############################################################################################################################################

class DkmCompGraph(object):
    def __init__(self, ae_specs, n_clusters, val_lambda):
        input_size = ae_specs[0][-1]
        embedding_size = ae_specs[0][int((len(ae_specs[0])-1)/2)]

        # Placeholder tensor for input data
        self.input = tf.placeholder(dtype=TF_FLOAT_TYPE, shape=(None, input_size))

        # Auto-encoder loss computations
        self.embedding, self.output = autoencoder(self.input, ae_specs)  # Get the auto-encoder's embedding and output
        rec_error = g_func(self.input, self.output)  # Reconstruction error based on distance g

        # k-Means loss computations
        ## Tensor for cluster representatives
        minval_rep, maxval_rep = -1, 1 # Clusters rep (centroids) are initialized randomly using U(-1,1)
        self.cluster_rep = tf.Variable(tf.random_uniform([n_clusters, embedding_size],
                                                    minval=minval_rep, maxval=maxval_rep,
                                                    dtype=TF_FLOAT_TYPE), name='cluster_rep', dtype=TF_FLOAT_TYPE)

        ## First, compute the distance f between the embedding and each cluster representative
        list_dist = []
        for i in range(0, n_clusters):
            dist = f_func(self.embedding, tf.reshape(self.cluster_rep[i, :], (1, embedding_size)))
            list_dist.append(dist)
        self.stack_dist = tf.stack(list_dist)

        ## Second, find the minimum squared distance for softmax normalization
        min_dist = tf.reduce_min(list_dist, axis=0)

        ## Third, compute exponentials shifted with min_dist to avoid underflow (0/0) issues in softmaxes
        self.alpha = tf.placeholder(dtype=TF_FLOAT_TYPE, shape=())  # Placeholder tensor for alpha
        list_exp = []
        for i in range(n_clusters):
            exp = tf.exp(-self.alpha * (self.stack_dist[i] - min_dist))
            list_exp.append(exp)
        stack_exp = tf.stack(list_exp)
        sum_exponentials = tf.reduce_sum(stack_exp, axis=0)

        ## Fourth, compute softmaxes and the embedding/representative distances weighted by softmax
        list_softmax = []
        list_weighted_dist = []
        for j in range(n_clusters):
            softmax = stack_exp[j] / sum_exponentials
            weighted_dist = self.stack_dist[j] * softmax
            list_softmax.append(softmax)
            list_weighted_dist.append(weighted_dist)
        stack_weighted_dist = tf.stack(list_weighted_dist)

        # Compute the full loss combining the reconstruction error and k-means term
        self.ae_loss = tf.reduce_mean(rec_error)
        self.kmeans_loss = tf.reduce_mean(tf.reduce_sum(stack_weighted_dist, axis=0))
        self.loss = self.ae_loss + val_lambda * self.kmeans_loss

        # The optimizer is defined to minimize this loss
        optimizer = tf.train.AdamOptimizer() 
        #optimizer = tf.train.AdadeltaOptimizer()
        self.pretrain_op = optimizer.minimize(self.ae_loss) # Pretrain the autoencoder before starting DKM
        self.train_op = optimizer.minimize(self.loss) # Train the whole DKM model

# MAIN PROGRAM ################################################################################################################################
###############################################################################################################################################

# This program executes DEEP K-means for a chosen dataset. 
# It runs the algorithm 10 times, each with different seed.
# It computes the final perfomance metrics selecting the iteration with the 
# maximum Accuracy. The performance metrics are printed in the screen as 
# they are computed. Also, a txt file is generated with these results
# and 2 png file with the confusion matrices are created.

if __name__=='__main__':
    
# STEP 0: Setting the directory -------------------------------------------------------------------------------------------------------------
    directory = input('Please enter the directory where you have the datasets: ')
    os.chdir(directory)
    
# STEP 1 Importing dataset ------------------------------------------------------------------------------------------------------------------    
    dataset = input('Select the dataset for which you want to run DEEP K-Means, digits/HAR/mushrooms [d/h/m]: ')
    if dataset=='d':
        import digits_specs as specs 
    elif dataset=='h':
        import HAR_specs as specs 
    elif dataset=='m':
        import mushrooms_specs as specs 
    else: 
        print('Enter a valid dataset [d/h/m]: ')

# STEP 2: Running DEEP K-Means ------------------------------------------------------------------------------------------------------------------
    
    # 2.1. Setting initial parameters------------------------------------------------------------------------------------------------------------
    
    n_pretrain_epochs = 50
    n_finetuning_epochs = 5
    lambda_ = 0.1 # Optimal lambda tested on USPS (similar database)
    batch_size = 256 # Size of the mini-batches used in the stochastic optimizer
    n_batches = int(math.ceil(specs.n_samples / batch_size)) # Number of mini-batches
    seeded = True # Specify if runs are seeded
    max_n = 20  # Number of alpha values to consider (constant values are used here), 100 epochs in total 20*5
    alphas = 1000*np.ones(max_n, dtype=float) # alpha is constant, all epochs take constant alpha of 1000
    seeds = [8905, 9129, 291, 4012, 1256, 6819, 4678, 6971, 1362, 575]  # Use a fixed seed for this run, as defined in the seed list 
                                                                        # There is randomness in centroids initialization and mini-batch selection
    n_runs = 10

    target = specs.target
    data = specs.data
    
    dict_acc = dict()
    dict_sc = dict()
    dict_cs = dict()
    dict_hs = dict()
    dict_clusterassign= dict()

    gpu_options = tf.GPUOptions(per_process_gpu_memory_fraction=0.05)
    config = tf.ConfigProto(gpu_options=gpu_options)

    # 2.2. Running k-means and creating the output files-----------------------------------------------------------------------------------------
    
    if os.path.exists('./results/DKM_results.txt'):
        os.remove('./results/DKM_results.txt') # Removes the text file if exists
    with open('./results/DKM_results.txt', 'x') as f: 
        
        # 2.2.1. Running k-means n_run times-------------------------------------------------------------------------------------------------------
       
        for run in range(n_runs):
            if seeded:
                tf.reset_default_graph()
                tf.set_random_seed(seeds[run])
                np.random.seed(seeds[run])
            print("Run", run)
            f.write("RUN "+ str(run) +'\n'+'\n') 
        # Define the computation graph for DKM
        # We need a computation graph that performs the simultaneous minimization of clusters and low-dim representation, so we have to program a specific low-level graph 
            cg = DkmCompGraph([specs.dimensions, specs.activations, specs.names], specs.n_clusters, lambda_)
        # DkmCompGraph([[500, 500, 2000, 10,
        #                2000, 500, 500, 64], 
        #               [tf.nn.relu, tf.nn.relu, tf.nn.relu, None, 
        #                tf.nn.relu, tf.nn.relu, tf.nn.relu, None], 
        #               ['enc_hidden_1', 'enc_hidden_2', 'enc_hidden_3', 'embedding', 
        #                'dec_hidden_1', 'dec_hidden_2', 'dec_hidden_3', 'output'] ], 
        #               10, 1.0)
            
        
        # Run the computation graph
            with tf.Session(config=config) as sess:
            # Initialization
                init = tf.global_variables_initializer()
                sess.run(init)
            
            # Variables to save tensor content
                distances = np.zeros((specs.n_clusters, specs.n_samples))
                print("Starting autoencoder pretraining...")
                
                # Variables to save pretraining tensor content
                embeddings = np.zeros((specs.n_samples, specs.embedding_size), dtype=float)
                
                # First, pretrain the autoencoder
                ## Loop over epochs
                for epoch in range(n_pretrain_epochs):
                    print("Pretraining step: epoch {}".format(epoch))
                    
                    # Loop over the samples
                    for _ in range(n_batches):
                        # Fetch a random data batch of the specified size
                        indices, data_batch = next_batch(batch_size, data)

                        # Run the computation graph until pretrain_op (only on autoencoder) on the data batch
                        _, embedding_, ae_loss_ = sess.run((cg.pretrain_op, cg.embedding, cg.ae_loss),
                                                       feed_dict={cg.input: data_batch})

                        # Save the embeddings for batch samples
                        for j in range(len(indices)):
                            embeddings[indices[j], :] = embedding_[j, :]
                            
                            #print("ae_loss_:", float(ae_loss_))

                # Second, run k-means++ on the pretrained embeddings
                print("Running k-means on the learned embeddings...")
                f.write("1. K-MEANS ON THE PRETRAINED FEATURES"+'\n'+'\n')
                kmeans_model = KMeans(n_clusters=specs.n_clusters, init="k-means++").fit(embeddings) # n_init=10 initial centroids
                
                    # Evaluate the clustering performance using the ground-truth labels
                sc = silhouette_score(data, kmeans_model.labels_)
                print("SC", sc)
                f.write("SC: "+ str(sc) +'\n')
                
                acc = cluster_acc(target, kmeans_model.labels_, specs.n_clusters)
                print("ACC", acc)
                f.write("ACC: "+ str(acc) +'\n')
                
                cs = completeness_score(target, kmeans_model.labels_)
                print("CS", cs)
                f.write("CS: "+ str(cs) +'\n')
                
                hs = homogeneity_score(target, kmeans_model.labels_)
                print("HS", hs)
                f.write("HS: "+ str(hs) +'\n')
                    
                # The cluster centers are used to initialize the cluster representatives in DKM
                sess.run(tf.assign(cg.cluster_rep, kmeans_model.cluster_centers_))
                
            # Train the full DKM model
                if (len(alphas) > 0):
                    print("Starting DKM training...")
                    f.write('\n') 
                    f.write("2. DEEP K-MEANS TRAINING"+'\n'+'\n') 
                ## Loop over alpha (inverse temperature), from small to large values
                for k in range(len(alphas)):
                    print("Training step: alpha[{}]: {}".format(k, alphas[k]))
                    f.write("TRAINING STEP: {}".format(k)+'\n'+'\n')
                    # Loop over epochs per alpha
                    for _ in range(n_finetuning_epochs):
                        # Loop over the samples
                        for _ in range(n_batches):
                            #print("Training step: alpha[{}], epoch {}".format(k, i))
                            
                            # Fetch a random data batch of the specified size
                            indices, data_batch = next_batch(batch_size, data)
                            
                            #print(tf.trainable_variables())
                            #current_batch_size = np.shape(data_batch)[0] # Can be different from batch_size for unequal splits
                            
                            # Run the computation graph on the data batch
                            _, loss_, stack_dist_, cluster_rep_, ae_loss_, kmeans_loss_ =\
                            sess.run((cg.train_op, cg.loss, cg.stack_dist, cg.cluster_rep, cg.ae_loss, cg.kmeans_loss),
                                     feed_dict={cg.input: data_batch, cg.alpha: alphas[k]})
                            
                            # Save the distances for batch samples
                            for j in range(len(indices)):
                                distances[:, indices[j]] = stack_dist_[:, j]
                                
                    # Evaluate the clustering performance every print_val alpha and for last alpha
                    print_val = 1
                    if k % print_val == 0 or k == max_n - 1:
                        print("loss:", loss_)
                        f.write("Total loss: "+ str(loss_) +'\n') 
                        print("ae loss:", ae_loss_)
                        f.write("AE loss: "+ str(ae_loss_) +'\n') 
                        print("kmeans loss:", kmeans_loss_)
                        f.write("KMeans loss: "+ str(kmeans_loss_) +'\n') 
                                    
                        # Infer cluster assignments for all samples
                        cluster_assign = np.zeros((specs.n_samples), dtype=float)
                        for i in range(specs.n_samples):
                            index_closest_cluster = np.argmin(distances[:, i])
                            cluster_assign[i] = index_closest_cluster
                        cluster_assign = cluster_assign.astype(np.int64)
                        
                        # Evaluate the clustering performance using the ground-truth labels
                        sc = silhouette_score(data, kmeans_model.labels_)
                        print("SC", sc)
                        f.write("SC: "+ str(sc) +'\n')
                            
                        acc = cluster_acc(target, kmeans_model.labels_, specs.n_clusters)
                        print("ACC", acc)
                        f.write("ACC: "+ str(acc) +'\n')
                        
                        cs = completeness_score(target, kmeans_model.labels_)
                        print("CS", cs)
                        f.write("CS: "+ str(cs) +'\n')
                
                        hs = homogeneity_score(target, kmeans_model.labels_)
                        print("HS", hs)
                        f.write("HS: "+ str(hs) +'\n')
                        f.write('\n') 
                            
        # Record the clustering performance for the run
                dict_acc.update({run : acc})
                dict_sc.update({run : sc})
                dict_cs.update({run : cs})
                dict_hs.update({run : hs})
                dict_clusterassign.update({run : kmeans_model.labels_ })
            f.write('\n'+'-------------------------------------------------------------'+'\n') 
        
            
        max_acc = max(dict_acc, key=dict_acc.get)
        print('Iteration with maximum ACC: '+ str(max_acc))
        print("ACC: {:.3f}".format(dict_acc[max_acc]))
        print("SC: {:.3f}".format(dict_sc[max_acc]))
        print("CS: {:.3f}".format(dict_cs[max_acc]))
        print("HS: {:.3f}".format(dict_hs[max_acc]))
        f.write('Run with maximum ACC: '+ str(max_acc) +'\n')
        f.write("ACC: {:.3f}".format(dict_acc[max_acc]) +'\n')
        f.write("SC: {:.3f}".format(dict_sc[max_acc]) +'\n')
        f.write("CS: {:.3f}".format(dict_cs[max_acc]) +'\n')
        f.write("HS: {:.3f}".format(dict_hs[max_acc]) +'\n')
        save_confusion_matrix(target, dict_clusterassign[max_acc], specs.n_clusters, 'DKM_CM') 
   



<file_sep>import pandas as pd
import matplotlib.pyplot as plt
import os
os.chdir('C://Users//vnino//Desktop//CE888 - Assignment 2//Python code//mushrooms')

names=['label', 'cap_shape', 'cap_surface', 
       'cap_color', 'bruises', 'odor', 
       'gill_attachment', 'gill_spacing', 'gill_size',
       'gill_color', 'stalk_shape', 'stalk_root', 
       'stalk_surface_above', 'stalk_surface_below', 
       'stalk_color_above', 'stalk_color_below', 'veil_type',
       'veil_color', 'ring_number', 'ring_type',
       'spore_print_color', 'population','habitat']
df = pd.read_csv("agaricus-lepiota.csv", names=names)

#Checking shape and head of the data
df.shape
df.head()

# Descriptive stats
df.describe()

#Distribution of categories of the features
for var in df.columns:
    print(df[var].value_counts()/8124)
    
# Plot of classes distribution
plot1=df['label'].value_counts().plot(kind='barh', color='0.75', title='Number of instances by edibility label', figsize=(7,3));

for i, v in enumerate(df['label'].value_counts()):
   plot1.text(v + 2, i-0.15 , str(v), fontsize='large')
plot1<file_sep># Import packages
import os
import tensorflow as tf
from sklearn import datasets
os.chdir('C://Users//vnino//Desktop//CE888 - Assignment 2//Python code//digits')
import random

# Fetch the dataset
dataset = datasets.load_digits()
data = dataset.data
target = dataset.target
n_samples = data.shape[0] # Number of obs in the dataset
n_clusters = 10 # Number of clusters to obtain

# Pre-process the dataset
data = data / 16.0 # Normalize the levels of grey between 0 and 1

# Get the split between training/test set and validation set
random.seed(0)
test_indices = random.sample(range(0, 1797), 360)
validation_indices = [x for x in range(0, 1797) if x not in test_indices]

# Auto-encoder architecture
input_size = data.shape[1]
hidden_1_size = 500
hidden_2_size = 500
hidden_3_size = 2000
embedding_size = n_clusters
dimensions = [hidden_1_size, hidden_2_size, hidden_3_size, embedding_size, # Encoder layer dimensions
              hidden_3_size, hidden_2_size, hidden_1_size, input_size] # Decoder layer dimensions
activations = [tf.nn.relu, tf.nn.relu, tf.nn.relu, None, # Encoder layer activations # None is linear activation
               tf.nn.relu, tf.nn.relu, tf.nn.relu, None] # Decoder layer activations # None is linear activation
names = ['enc_hidden_1', 'enc_hidden_2', 'enc_hidden_3', 'embedding', # Encoder layer names
         'dec_hidden_1', 'dec_hidden_2', 'dec_hidden_3', 'output'] # Decoder layer names

activations_soft = [tf.nn.relu, tf.nn.relu, tf.nn.relu, tf.nn.softmax, # Encoder layer activations # Uses softmax in last enocder layer
               tf.nn.relu, tf.nn.relu, tf.nn.relu, None]<file_sep>import matplotlib
matplotlib.use('Agg')
import pandas as pd
import seaborn as sns
import numpy as np

## PART 1: Creating the boostrap function

def boostrap(sample, sample_size, iterations,ciparam):
    sample_matrix=np.zeros((iterations,sample_size))
    data_mean_mx=np.zeros(iterations)
    for i in range(iterations):
        sample_matrix[i, :]=np.random.choice(sample, size=sample_size, replace=True)
        data_mean_mx[i]=np.mean(sample_matrix[i, :])
    lower=np.percentile(data_mean_mx,ciparam/2)
    upper=np.percentile(data_mean_mx,100-ciparam/2)
    data_mean=np.mean(data_mean_mx)
    return data_mean, lower, upper

## PART 2 Applying the boostrap function on the vehicles data

vehicles_df = pd.read_csv('./vehicles.csv')
vehicles_df.columns=['Currentfleet','Newfleet']
vehicles_df2=vehicles_df[vehicles_df.Newfleet.notnull()]
current_fleet_data = vehicles_df.values.T[0]
new_fleet_data = vehicles_df2.values.T[1]

cf=boostrap(current_fleet_data, current_fleet_data.shape[0],1000,5)

print(("Current Fleet Mean: %f")%(cf[0]))
print(("Current Fleet Lower: %f")%(cf[1]))
print(("Current Fleet Upper: %f")%(cf[2]))

nf=boostrap(new_fleet_data, new_fleet_data.shape[0],1000,5)

print(("New Fleet Mean: %f")%(nf[0]))
print(("New Fleet Lower: %f")%(nf[1]))
print(("New Fleet Upper: %f")%(nf[2]))

## The New fleet has a higher MPG mean, the 95% CI for the mean of the new fleet 
# is between 29.1 MPG and 31.7 MPG, while for the Current Fleet it is between 19.4 and 20.9.

## Plots on salaries using boostrap iterations

if __name__ == "__main__":
	df = pd.read_csv('./salaries.csv')

	data = df.values.T[1]
	boots = []
	for i in range(100, 100000, 1000):
		boot = boostrap(data, data.shape[0], i, 5)
		boots.append([i, boot[0], "mean"])
		boots.append([i, boot[1], "lower"])
		boots.append([i, boot[2], "upper"])

	df_boot = pd.DataFrame(boots, columns=['Boostrap Iterations', 'Mean', "Value"])
	sns_plot = sns.lmplot(df_boot.columns[0], df_boot.columns[1], data=df_boot, fit_reg=False, hue="Value")

	sns_plot.axes[0, 0].set_ylim(0,)
	sns_plot.axes[0, 0].set_xlim(0, 100000)

	sns_plot.savefig("bootstrap_confidence.png", bbox_inches='tight')
	sns_plot.savefig("bootstrap_confidence.pdf", bbox_inches='tight')


	#print ("Mean: %f")%(np.mean(data))
	#print ("Var: %f")%(np.var(data))
	



	<file_sep>#Source: https://github.com/srvds/Human-Activity-Recognition/blob/master/HAR_EDA.ipynb

# Import packages
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import os
os.chdir('C://Users//vnino//Desktop//CE888 - Assignment 2//Python code//HAR')

# Fetch the dataset
train_df= pd.read_csv("train.csv")
test_df= pd.read_csv("test.csv")
data_df=pd.concat([train_df,test_df], ignore_index=True)

#Defining labels and features
target=data_df['Activity']
data=np.asarray(data_df.iloc[:,0:560])
n_samples= data.shape[0]


#Distribution of the mean of Linear Body Acceleration Euclidean Norm (tBodyAccMagmean)
sns.set_palette("Set1", desat=0.80)
facetgrid = sns.FacetGrid(data_df, hue='ActivityName', size=6,aspect=2)
facetgrid.map(sns.distplot,'tBodyAccMagmean', hist=False)\
    .add_legend()
plt.annotate("Stationary Activities", xy=(-0.956,17), xytext=(-0.9, 23), size=20,\
            va='center', ha='left',\
            arrowprops=dict(arrowstyle="simple",connectionstyle="arc3,rad=0.1"))

plt.annotate("Moving Activities", xy=(0,3), xytext=(0.2, 9), size=20,\
            va='center', ha='left',\
            arrowprops=dict(arrowstyle="simple",connectionstyle="arc3,rad=0.1"))
plt.show()

#Boxplot of the mean of Linear Body Acceleration Euclidean Norm (tBodyAccMagmean)
sns.boxplot(x='ActivityName', y='tBodyAccMagmean', data = data_df, showfliers=False)
plt.xticks(rotation = 40)
plt.axhline(y=-0.22, xmin=0.1, xmax=0.8, dashes=(5,3), c='m')
plt.xticks(rotation=90)
plt.show()

#Boxplot of the mean of Angle Y-Gravity (tBodyAccMagmean)
sns.boxplot(x='ActivityName', y='angleYgravityMean', data = data_df, showfliers=False)
plt.xticks(rotation = 40)
plt.axhline(y=-0.22, xmin=0.1, xmax=0.8, dashes=(5,3), c='m')
plt.xticks(rotation=90)
plt.show()<file_sep>from sklearn import datasets
import pandas as pd
import matplotlib.pyplot as plt

digits = datasets.load_digits()

print(digits.data.shape)
print(digits.data)
print(digits.target)

digits_df=pd.DataFrame(digits.data)
digits_df.describe()

diglab_df=pd.DataFrame(digits.target)
diglab_df.describe()

# Count of the classes
pd.DataFrame(target)[0].value_counts()

# Plot of classes distribution
plot1=diglab_df[0].value_counts().plot(kind='barh', color='0.75', title='Number of instances by digit', figsize=(7,3));
values=diglab_df[0].value_counts()/1797
values_df=values.to_frame()
values_df['Count']= pd.Series(["{0:.2f}%".format(val * 100) for val in values_df[0]], index = values_df.index)
for i, v in enumerate(diglab_df[0].value_counts()):
   plot1.text(v + 1, i-0.25 , str(v), fontsize='x-small')
plot1

# Plot of the images of the digits
plt.matshow(digits.images[9], cmap='binary')

<file_sep># For creating this code the following sources were used: 
# 1. Implementation of "Deep k-Means: Jointly Clustering with k-Means and Learning Representations" by <NAME>, <NAME>, and <NAME>. (https://github.com/MaziarMF/deep-k-means)
# 2. In Depth: k-Means Clustering (https://jakevdp.github.io/PythonDataScienceHandbook/05.11-k-means.html)

# IMPORTS PACKAGES ############################################################################################################################
###############################################################################################################################################

import os
import numpy as np
import sklearn
from sklearn.metrics import silhouette_score, accuracy_score, completeness_score, homogeneity_score
from sklearn.cluster import KMeans
from scipy.stats import mode
import seaborn as sns
sns.set(font_scale=3)
import matplotlib.pyplot as plt

# FUNCTIONS DEFINITION ########################################################################################################################
###############################################################################################################################################

def cluster_acc(y_true, y_pred, n_clusters):
    labels = np.zeros_like(y_pred)
    for i in range(n_clusters):
        mask = (y_pred == i)
        labels[mask] = mode(y_true[mask])[0]
    if y_true is not None:
        acc = np.round(accuracy_score(y_true, labels), 5)
    return acc


def save_confusion_matrix(y_true, y_pred, n_clusters, name):
    labels = np.zeros_like(y_pred)
    for i in range(n_clusters):
        mask = (y_pred == i)
        labels[mask] = mode(y_true[mask])[0]
    confusion_matrix = sklearn.metrics.confusion_matrix(y_true, labels)
    plt.figure(figsize=(16, 14))
    sns.heatmap(confusion_matrix, annot=True, fmt="d", annot_kws={"size": 20});
    plt.title("Confusion matrix", fontsize=30)
    plt.ylabel('True label', fontsize=25)
    plt.xlabel('Clustering label', fontsize=25)
    plt.savefig('./results/'+ name + '.png')

# MAIN PROGRAM ################################################################################################################################
###############################################################################################################################################

# This program executes the traditional K-means for a chosen dataset. 
# It runs the algorithm 10 times, each with different seed.
# It computes the final perfomance metrics selecting the iteration with the 
# maximum Accuracy. The performance metrics are printed in the screen as 
# they are computed. Also, a txt file is generated with these results
# and 2 png file with the confusion matrices are created.

if __name__=='__main__':
    
# STEP 0: Setting the directory -------------------------------------------------------------------------------------------------------------
    directory = input('Please enter the directory where you have the datasets: ')
    os.chdir(directory)

# STEP 1 Importing dataset ------------------------------------------------------------------------------------------------------------------    
    dataset = input('Select the dataset for which you want to run K-Means, digits/HAR/mushrooms [d/h/m]: ')
    if dataset=='d':
        import digits_specs as specs 
    elif dataset=='h':
        import HAR_specs as specs 
    elif dataset=='m':
        import mushrooms_specs as specs 
    else: 
        print('Enter a valid dataset [d/h/m]: ')

# STEP 2: Running K-Means ------------------------------------------------------------------------------------------------------------------
    
    # 2.1. Setting initial parameters------------------------------------------------------------------------------------------------------------

    seeded = True # Specify if runs are seeded
    seeds = [8905, 9129, 291, 4012, 1256, 6819, 4678, 6971, 1362, 575] #Specific seeds used
    n_runs = 10 # How many times should we run the algorithm
    
    dict_acc = dict()
    dict_sc = dict()
    dict_cs = dict()
    dict_hs = dict()
    dict_clusterassign= dict()
    
    target = specs.target
    data = specs.data

    # 2.2. Running k-means and creating the output files-----------------------------------------------------------------------------------------
    
    if os.path.exists('./results/KM_results.txt'):
        os.remove('./results/KM_results.txt') # Removes the text file if exists
    with open('./results/KM_results.txt', 'x') as f: 
        
        # 2.2.1. Running k-means n_run times-------------------------------------------------------------------------------------------------------
        
        for run in range(n_runs):
            np.random.seed(seeds[run])
            # Shuffle the dataset
            print("Run", run)
            f.write("RUN "+ str(run) +'\n'+'\n') 
            
            # Run k-means(++) on the original data

            kmeans_model = KMeans(n_clusters=specs.n_clusters, init="k-means++").fit(data)
                
            # Evaluate the clustering
            
            cluster_assign = np.asarray(kmeans_model.labels_)
                
            # Evaluate the clustering
            
            sc = silhouette_score(data, cluster_assign)
            print("SC", sc)
            f.write("SC: "+ str(sc) +'\n')
                
            acc = cluster_acc(target, cluster_assign, specs.n_clusters)
            print("ACC", acc)
            f.write("ACC: "+ str(acc) +'\n')
                
            cs = completeness_score(target, cluster_assign)
            print("CS", cs)
            f.write("CS: "+ str(cs) +'\n')
                
            hs = homogeneity_score(target, cluster_assign)
            print("HS", hs)
            f.write("HS: "+ str(hs) +'\n')
                
            # Save performance scores and predictions
            
            dict_acc.update({run : acc})
            dict_sc.update({run : sc})
            dict_cs.update({run : cs})
            dict_hs.update({run : hs})
            dict_clusterassign.update({run : cluster_assign})
            
            f.write('\n'+'-------------------------------------------------------------'+'\n') 
                
        # 2.2.2. Getting the iteration with the best Accuracy Score---------------------------------------------------------------------------------
        
        max_acc = max(dict_acc, key=dict_acc.get)
        print('Iteration with maximum ACC: '+ str(max_acc))
        print("ACC: {:.3f}".format(dict_acc[max_acc]))
        print("SC: {:.3f}".format(dict_sc[max_acc]))
        print("CS: {:.3f}".format(dict_cs[max_acc]))
        print("HS: {:.3f}".format(dict_hs[max_acc]))
        f.write('Run with maximum ACC: '+ str(max_acc) +'\n')
        f.write("ACC: {:.3f}".format(dict_acc[max_acc]) +'\n')
        f.write("SC: {:.3f}".format(dict_sc[max_acc]) +'\n')
        f.write("CS: {:.3f}".format(dict_cs[max_acc]) +'\n')
        f.write("HS: {:.3f}".format(dict_hs[max_acc]) +'\n')
        save_confusion_matrix(target, dict_clusterassign[max_acc], specs.n_clusters, 'KM_CM') 
    
 