# ---
# 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
# ---

# +
import pandas as pd
import numpy as np
import math
import matplotlib.pyplot as plt
import time
import pickle

pd.set_option('display.max_columns', None)
pd.set_option('display.max_rows', None)
from yellowbrick.cluster import SilhouetteVisualizer
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans
from sklearn import metrics    
from sklearn.cluster import KMeans
from sklearn.metrics.cluster import adjusted_rand_score
from myFunctions import *
# -

# # Getting your data! & Choosing your features

# ## Echonest

#read echonest.csv
echonest=pd.read_csv("data/echonest.csv")
echonest.shape

echonest.tail()

# ### A little description of some features:
# - track_id: the id of the track
# - audio_features_acousticness: a value in [0-1], when the track is closer to 0 it has mainly an electric sound
# - audio_features_danceability: a value in [0-1], it's expressing a combination of tempo,rythm, beat strengt, so how suitable is for dancing
# - audio_features_energy: a value in [0-1], when the track is closer to 1, the song has fast,loud,noisy sounds
# - audio_features_instrumentalness: a value in [0-1], when the track is closer to 1, the song has the high probability that there is no vocals
# - audio_features_liveness: a value in [0-1], it describes how likely the song was performed live
# - audio_features_speechiness: a value in [0-1], it describes how it is exclusively speech-like the recording
# - audio_features_tempo: it describes the number of beats per minute
# - audio_features_valence:a value in [0-1], it describes how much the song is positive
# - ranks_artist_familiarity_rank: a value in [0-1], it is an indication of how well known the artist is 
# - ranks_artist_hotttnesss_rank: a value in [0-1], it gives an indication how much popular the artist is  now
# - ranks_song_hotttnesss_rank: a value in [0-1], it gives an indication how much popular the song is now

# Now, we are going to drop the columns that are not needed.

# make a copy of our dataframa
e=echonest.copy()
not_dropped_columns=["audio_features_acousticness","audio_features_danceability","audio_features_energy","audio_features_instrumentalness","audio_features_liveness","audio_features_speechiness","audio_features_tempo","audio_features_valence"]
e=e[not_dropped_columns]

#  The next step is counting the NA values after the dropping:

e.isnull().sum()

# We use PCA to make the dataset smaller and to remove the redundancy
#
# this choice was taken for all the datasets.
# the reasons are many:
# - this removes the correlated features
# - it improves the algorithm performance
# - it gives the possibility to choose the quantity of information to keep
#
# However it has some disavantages:
# - indenpendent variables become hard to find
# - there is a limited information loss
#
# First step is Standardization, to make that each variable  contributes equally to the analysis.
# $$
#   z=\frac {value - mean } {standardDeviation}\
# $$

normalised_columns=["audio_features_acousticness","audio_features_danceability","audio_features_energy","audio_features_instrumentalness","audio_features_liveness","audio_features_speechiness","audio_features_tempo","audio_features_valence"]
for col in normalised_columns:
    mean=e[col].mean()
    std=e[col].std()
    e[col]=(e[col]- mean)/std
e.tail()

e.shape

# I compute the PCA on this dataset, creating principals to get a total variance that is more than 70 %.

pca = PCA(n_components = 0.70)
X_new = pca.fit_transform(e)
x=pca.transform(e)
x.shape

# Here we create a new dataframe and insert again the track_id

e=pd.DataFrame(x)
e= e.add_suffix("_e")
e["track_id"]=echonest["track_id"]
e.tail()

# ## Features.csv

#read features.csv
features=pd.read_csv("data/features.csv")
features.shape

features.tail()

# <div>
# <img src = "geatures.png" width = "1000px">
# </div>

# Here we count the NA values

any(features.isnull().sum())

# There aren't NA value in this dataset. So now we do again the standardization.

for col in features:
    if col!="track_id":
        mean=features[col].mean()
        std=features[col].std()
        features[col]=(features[col]- mean)/std
features.tail()

f=features.loc[:, features.columns != 'track_id']

# We proceed with the PCA on this dataset

pca = PCA(n_components = 0.70)
X_new = pca.fit_transform(f)
x=pca.transform(f)
x.shape

f=pd.DataFrame(x)
f= f.add_suffix("_f")
f["track_id"]=features["track_id"]
f.tail()

# ## Tracks.csv

tracks=pd.read_csv("data/tracks.csv")
tracks.shape

tracks.head()

#Let's check which are the top genre
print(*{x for x in set(tracks["track_genre_top"]) if isinstance(x , str) or not math.isnan(x)}, sep = "\n")

# ### Now a list of all the attributes and their NA values:

tracks.isnull().sum()

# I selected the columns according to columns that can be used with pca and columns that have few NA values.

t=tracks.copy()
not_dropped_columns=["track_bit_rate","track_duration" ]
t=t[not_dropped_columns]
t.head()

t.isnull().sum()

# We use PCA to make the dataset smaller and to remove the redundancy
#
# First step is Standardization, to make that each variable  contributes equally to the analysis.
# $$
#   z=\frac {value - mean } {standardDeviation}\
# $$

cols= ["track_bit_rate","track_duration"]
for col in cols:
    mean=t[col].mean()
    std=t[col].std()
    t[col]=(t[col]- mean)/std
t.tail()

pca = PCA(n_components = 0.70)
X_new = pca.fit_transform(t)
x=pca.transform(t)

t=pd.DataFrame(x)
t= t.add_suffix("_t")
t["track_id"]=tracks["track_id"]
t["track_genres"]=tracks["track_genres"]
t["track_genres_all"]=tracks["track_genres_all"]
t["track_genre_top"]=tracks["track_genre_top"]
t.tail()

# ### Here we merge all the datasets using the track_id

new_df = e.merge(f, on = ["track_id"])
new_df=new_df.merge(t,on=["track_id"])
new_df.tail()

#save our result in a csv file
new_df.to_csv("final_dataset.csv",index=False)


# # Clustering

# Now we are going to apply k_means to this dataset
#
# Our k_means choose randomly k points as the first centroids and then associate each point in the cluster to the cluster whose centroid is the nearest.
# After these steps it computes again new centroids using the mean and it redoes all the distances again.
# this is done until for 2 successive steps we have a difference equal to 0 between the new centroids and the old ones.

def K_scratch(dataframe,k,iters):
    #it moves the keys at the end of the rows
    ids=dataframe["track_id"].copy()
    del dataframe["track_id"]
    dataframe["track_id"]=ids
    dataframe["cluster_#"]=0

    differences=[]
    
    #it takes k rows from the dataframe excluding the string columns
    centroids=dataframe.iloc[:,:-5].sample(n=k)
    first=False
    diff=1
    steps=1
    PastDistance=0
    dataframe["cluster_#"]=0
    while diff!=0 and steps<=iters:
        distanceSum=0
        steps+=1
        #it iterates on each row(point)
        for i in range(len(dataframe)):
            #it keeps the distance of each node to the k clusters
            support=np.zeros(k)
            #it iterates on each centroids
            for j in range(len(centroids)):
                #it fills with all the distances
                r=np.asarray(dataframe.iloc[i,:-5])
                c=np.asarray(centroids.iloc[j,:])
                dist = np.linalg.norm(r-c)
                support[j]=dist
            #it associates the  i point to the cluster that minimize the distance that corrisponds to the index of support
            index=np.argmin(support)
            #it measure the total distance from each point to the centroid that is the nearest
            distanceSum+=support[index]
            dataframe.loc[i,"cluster_#"]=index
            #new centroids
        differences.append(distanceSum)   
        new_centroids=dataframe.groupby(["cluster_#"],group_keys=False).mean().iloc[:,:-1]
        if first:
            v=new_centroids-centroids
            diff=np.sum(np.sum(v,axis=1))
            PastDistance=distanceSum
            centroids=new_centroids.copy()
        else:
            diff=1
            PastDistance=distanceSum
            centroids=new_centroids.copy()
            first=True
    vect=np.asarray(differences)
    return vect


def elbow_method(dataframe):
    distanceMeasures=[]
    for k in range(2,20):
        vec=K_scratch(dataframe,k,10)
        distance=vec[-1]
        distanceMeasures.append(distance)
    plt.plot(distanceMeasures,"b-o")
    plt.grid(visible=True)
    plt.xticks(range(0,18), range(2,20))
    plt.show()


dataset=pd.read_csv("final_dataset.csv")

elbow_method(dataset)

# Using the elbow method seems that 11 is the right number to choose 
#
# Now we use the other chosen method: silhoutte method

# +
from sklearn.metrics import silhouette_score

sil = []
kmax = 19

# dissimilarity would not be defined for a single cluster, thus, minimum number of clusters should be 2
for k in range(2, kmax+1):
    K_scratch(dataset,k,10)
    labels = np.asarray(dataset["cluster_#"])
    sil.append(silhouette_score(X=dataset.iloc[:,:-5],labels=labels, metric='euclidean'))
plt.plot(sil,"b-o")
plt.grid(visible=True)
plt.xticks(range(0,18), range(2,20))
plt.show()
# -

# Both using the elbow method and the silhoutte method we learn that the best number of  cluster is 11!

# ## K-means++ implementation

datasetplus=pd.read_csv("final_dataset.csv")

# we use the kmeans function from sklearn and output the labels
# %%time
kmeans = KMeans(n_clusters=11, random_state=0,init="k-means++",max_iter=10).fit(datasetplus.iloc[:,:-5])

kmeans.labels_

# run also our function
# %%time
vect=K_scratch(dataset,11,10)
plt.plot(vect,"b-o")
plt.xticks(range(0,10),range(1,11))
plt.show()

dataset["cluster_#"].to_numpy()

# to check between the two clustering methods we will use the adjusted_rand_score from sklearn
adjusted_rand_score(kmeans.labels_,dataset["cluster_#"].to_numpy())

# Well, turns out that our algorithm takes roughly 400 seconds to run, while the sklearn one runs in just below a quarter of a second. this makes the sklearn one roughly 2k times faster. Also the adjusted_rand_score is very close to 0 so the two predictions are very far from one another.

dataset.to_csv("cluster_dataset.csv",index=False)

# # Analysing my results!
#
# At first we run once more our own k-means algo with 11 clusters to obtain the results, then we create an inverted dictionary to assign to each cluster the corresponding songs that belong to it.
#
# The plot is describing the inertia in comparison at the number of iterations

# +
# below the inverted index dictionary
diz ={}
df=dataset[["track_id","cluster_#"]]
for index in range (13128) :
    row=df.iloc[index,:]
    track_id=row[0]
    cluster_id=row[1]
    if cluster_id not in diz:
        diz[cluster_id]=[track_id]
    else:
        diz[cluster_id].append(track_id)

with open('data/dictionary_myKmean.pkl','wb') as f:
        pickle.dump(diz,f)
# -

# From here on out we begin the analysis. We decided to use 6 variables
# + audio_features_tempo
# + audio_features_danceability
# + audio_features_liveness
# + audio_features_speechiness
# + audio_features_instrumentalness
# + track_duration 

# #### Bitrate (audio_features_tempo)
# The first table shows that the majority of songs across all clusters is between 80 and 130 bitrates, with values ranging from 45% in clusters #2 and #10 to the 60% encountered in cluster #11. Interesting is also the case of cluster #3 where 1 out of five songs has a bitrate below 80.

# +
with open('data/dictionary_myKmean.pkl','rb') as f:
            diz = pickle.load(f)

printPivotTableBit(diz)
# -

# #### Danceability
# The second table here shows how some clusters are clearly more danceable than others. Clusters #1, #9 and #11 have above 65% of the songs with more that 0.50 danceability, where the top position is kept by #11 at over 80% cumulative. On the other hand clusters such as #2, #3, and #7 have very low values of danceability overall, in cluster #3 infact almost 75% of the songs have a danceability score below 0.50.

printPivotTableDanceability(diz)

# #### Liveness
# The variable may express the relative probability of the song to be live, and there doesn't seem to be much information here, all the cluster have very low liveness, in clusters #4 and #5 for example almost 9 out of 10 songs in the cluster has a minuscule liveness score.

printPivotTableLiveness(diz)

# #### Speechiness
# Here again the variable might express the presence of a determined set of frequencies that make up the human voice. As before all values are flattened in the first category, with a speechiness score below 0.25 that in some clusters grows as high as 98%! particular mention to cluster #11 which has overall higher values for this parameter.

printPivotTableSpeechiness(diz)

# #### Instrumentalness
# The instrumentalness score on the other hand revolves around the quantity of instrumentals in the song, we get the highest values in cluster #3 and the lowest in clusters #8 and #9. peculiar is the case of cluster #9 where songs seem to be either incredibly instrumental or not at all.

printPivotTableInstrumentalness(diz)

# #### Duration
# We divided the categories in songs shorter than 2 minutes, between 2 and 8, 8 and 12, and longer than 12. Of course the majority of the songs is in the second category but we can still derive some insight. For example cluster #10 has the highest number of short songs, while cluster #8 contains the longest ones.

printPivotTableDuration(diz)

# #### Genres
# *sidenote: since many songs had NaN values on the genre_top variable the following percentages were calculated according to the songs which had a non-null value as the maximum*
# After all the different variables seen in the last slides we come here to see the truth behind these clusters. Cluster #3 with low bitrate, danceability and speech is mostly made up of historical and classical music, that checks out. Clusters #1 and #9 are very similar, both have high values for rock, hip hop and electronic music. Clusters #4 and #7 are unique as a combination of folk-rock. Clusters #2 and #6 are also similar as both have very high values of rock music. Cluster #11 is filled with electronic and hip-hop, meaning that songs that mix the two genres are grouped here here.

printPivotTableGenres(diz)

# # K-means++ again
#
# Now we performed again the k-means++ analysis without using the variables from the file echonest.csv, so we take in account only datas from features.csv and tracks.csv. Let's create our dataset:

# +
#load file features
features=pd.read_csv("data/features.csv")

#normalize all the values
for col in features:
    if col!="track_id":
        mean=features[col].mean()
        std=features[col].std()
        features[col]=(features[col]- mean)/std

f=features.loc[:, features.columns != 'track_id']

#execute the PCA on features
pca = PCA(n_components = 0.70)
X_new = pca.fit_transform(f)
x=pca.transform(f)

f=pd.DataFrame(x)
f= f.add_suffix("_features")
f["track_id"]=features["track_id"]
f.head()

# +
#load file tracks
tracks=pd.read_csv("data/tracks.csv")
t=tracks.copy()
cols=["track_bit_rate","track_duration","track_listens", 
      "track_interest", "track_favorites", "track_duration",
      "track_bit_rate", "album_tracks",                     
      "album_listens", "artist_favorites"]
t=t[cols]

#normalize all the values
for col in cols:
    mean=t[col].mean()
    std=t[col].std()
    t[col]=(t[col]- mean)/std
t.head()

#execute the PCA on tracks
Y_new = pca.fit_transform(t)
y=pca.transform(t)

t=pd.DataFrame(x)
t= t.add_suffix("_t")
t["track_id"]=tracks["track_id"]
t["track_genres"]=tracks["track_genres"]
t["track_genres_all"]=tracks["track_genres_all"]
t["track_genre_top"]=tracks["track_genre_top"]
# -

#create the final dataset merging the twos just created
final_df = f.merge(t, on = ["track_id"])
final_df.head()

final_df.to_csv("final_dataset_k++.csv",index=False)
#finally we get a datset with ~100k rows
final_df.shape

# Now that we got our dataset, we perform the k-means++ algorithm and to decidece the optimal value of number of clusters we use the elbow method and the silhouette score:

# +
silhouette_scores = {}
distortions = {}

for k in range(2,20):
    model_kmeans_k = KMeans(n_clusters=k, init='k-means++', n_init=10, max_iter=100, random_state=42)
    model_kmeans_k.fit(final_df.iloc[:,:-5])
    labels_k = model_kmeans_k.labels_
    score_k = metrics.silhouette_score(final_df.iloc[:,:-5], labels_k)
    silhouette_scores[k] = score_k
    distortions[k] = model_kmeans_k.inertia_
    
    print("k = %d\tSS: %5.4f\t distortion: %5.4f\t" % (k, score_k, model_kmeans_k.inertia_))
    
print("Done!")
# -

plt.figure(figsize=(8,4))
plt.plot(distortions.values(), 'bo-')
plt.xticks(range(0,18), distortions.keys())
plt.grid(visible=True)
plt.xlabel('k')
plt.ylabel('Distortion')
plt.title('The Elbow Method showing the optimal k')
plt.show()

plt.figure(figsize = (8,5))
plt.plot(silhouette_scores.values(), 'bo-')
plt.xticks(range(0,18), silhouette_scores.keys())
plt.title("Silhouette Metric")
plt.grid(visible=True)
plt.xlabel("k")
plt.ylabel("Silhouette")
plt.show()

# From graph analysis we can say that the optimal value of number of clusters is: 10

# +
kmeans_model = KMeans(n_clusters = 10)
kmeans_model.fit(final_df.iloc[:,:-5])

final_df["cluster_#"] = kmeans_model.labels_

# +
dizz ={}
df=final_df[["track_id","cluster_#"]]
for index in range (len(final_df)) :
    row=df.iloc[index,:]
    track_id=row[0]
    cluster_id=row[1]
    if cluster_id not in dizz:
        dizz[cluster_id]=[track_id]
    else:
        dizz[cluster_id].append(track_id)
        
with open('data/dictionary_k++.pkl','wb') as f:
        pickle.dump(dizz,f)
# -

# Now that we got our dictionary that contains for each number of cluster the correspondent track_id, we can analyze the results obtained in order to see if we can recognize the different genres.
#
# We decided to analyze the next 5 features:
# - Duration: how much time a song lasts;
# - RMSEE mean: it express the mean effective power of the signal sound;
# - RMSEE max: it express the max effective power of the signal sound;
# - Spectral Bandwidth mean: the difference between the upper and lower frequencies;
# - Zero Crossing rate mean: is the rate at which a signal changes from positive to zero to negative or from negative to zero to positive.

# #### Duration

# +
with open('data/dictionary_k++.pkl','rb') as f:
            dizz = pickle.load(f)

printPivotTableDuration(dizz)
# -

# #### RMSEE mean

printPivotTableRMSEEmean(dizz)

# #### RMSEE max

printPivotTableRMSEEmax(dizz)

# #### Spectral Bandwidth

printPivotTableSpectral(dizz)

# #### Zero crossing rate

printPivotTableZero(dizz)

# Analyizing all the pivot table we notice that the distribution of the values are too homogenous and we can't tell which cluster corrispond to which genre.

StarcoderdataJupyterScriptsDedupFiltered