%pylab inline

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
import glob

pylab.rcParams['figure.figsize'] = (16, 12)

dataset = pd.read_csv('data/beatles/derived/all_chords_with_keys.tsv', sep='\t')
dataset = dataset.dropna()
dataset

"""
Explanation: Classification of key based on chords
The task is to classify each chord segment with a key.
Inputs  & outputs
Each data point represents a chord segment.
Inputs are chord segments with duration and a set of binary pitch classes contained in the chord. Thus the input dimension is 1+12 = 13.
Output of each data point is the key within the segment represented by its canonical diatonic key root pitch class. So there are 12 possible class values. With one-hot encoding the output dimension is 12.
Assumptions

assumptions to simplify the problem:
we assume that each segments has just a single chord and key
we ignore chord segments with no chord (N) or no key (silence)
we represent the key with its canonical diatonic key root pitch class
ie. no distinction of modes (major, minor, dorian, etc.) or accidentals (sharp, flat)


the inputs should be selected such that each input spans a single key, not multiple ones

The strategy of tasks - from simpler to harder
First we can try to classify each chord segment independently to establish a baseline performance. In practice this is not too realistic, since the key may be distinguisehd by a set of multiple chords.
Then we can try to use as inputs short sequences of adjacent chord segments of fixed size. This might be more realistic and be usable for non-recurrent ML models.
Ideally we should handle the fact that the sequence of chords is of a variable length and use recurrent models and not be constrained by a fixed window.
ML algorithms
As for the ML models we should start with a plain multi-class logistic regression.
Then we could try deep neural networks:
- architecture:
  - input layer: 13
  - hidden layer(s): ???
  - output layer: 12
    - soft-max classifer
  - regularization - dropout, ...
- pre-training with an autoencoder, then fine tuning with labels
- experiment with:
  - convolution layers (http://theanets.readthedocs.org/en/stable/reference.html#module-theanets.layers.convolution)
    - eg. we would like to capture several distance metrics on the pitch class space:
      - skip by 1 (plain pitch class neighborhood), skip by 7 (circle-of-fifths similarity)
    - for pitch spaces we could also use skip by 12 (octave similarity)
  - recurrent layers (http://theanets.readthedocs.org/en/stable/reference.html#module-theanets.layers.recurrent)
    - we could use the Predict output layer
Evaluation
TODO
End of explanation
"""


pcs_columns = ['C','Db','D','Eb','E','F','Gb','G','Ab','A','Bb','B']
# note: we could also try use features like 'root' or 'bass'
input_cols = pcs_columns
output_col = 'key_diatonic_root'

def select_X_y(df):
    X = df[input_cols].astype(np.int32)
    y = df[output_col].astype(np.int32)
    return X, y

X, y = select_X_y(dataset)

len(dataset)

X.head()

y.head()

from sklearn.cross_validation import train_test_split

def split_dataset(X, y, test_percentage=0.2, valid_percentage=0.2, random_state=42):
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=test_percentage, random_state=random_state)
    X_train, X_valid, y_train, y_valid = train_test_split(X_train, y_train, test_size=valid_percentage/(1-test_percentage), random_state=random_state)
    return (X_train, y_train), (X_valid, y_valid), (X_test, y_test)

((X_train, y_train), (X_valid, y_valid), (X_test, y_test)) = split_dataset(X, y)

[d.shape for d in (X_train, y_train, X_valid, y_valid, X_test, y_test)]

from sklearn.linear_model import LogisticRegression

lr_model = LogisticRegression()
lr_model.fit(X_train, y_train.ravel())

from sklearn.metrics import classification_report, confusion_matrix

def accuracy_report(model, X_train, y_train, X_valid, y_valid):
    # mean accuracy on the given test data and labels
    accuracy_train, accuracy_valid = model.score(X_train, y_train), model.score(X_valid, y_valid)
    print('training set accuracy (%):', accuracy_train)
    print('validation set accuracy (%):', accuracy_valid)
    print('difference (% points):', accuracy_train - accuracy_valid)

accuracy_report(lr_model, X_train, y_train, X_valid, y_valid)

# TODO: we could better use http://scikit-learn.org/stable/modules/generated/sklearn.dummy.DummyClassifier.html

random_batches = 100
y_random_uniform = np.random.random_integers(low=0, high=11, size=(len(X_valid), random_batches))
y_pred_lr_valid = lr_model.predict(X_valid)

def prior_class_probabilities(y):
    priors = pd.Series(np.zeros(12))
    hist = (pd.Series(y.ravel()).value_counts() / len(y))
    priors[hist.index] = hist
    return priors

priors_valid = prior_class_probabilities(y_valid)
priors_valid

# random samples of classes based on prior probabilities of each class
y_random_prior = numpy.random.choice(np.arange(12), size=(len(y_valid), random_batches), p=priors_valid)
pd.Series(y_random_prior[0]).value_counts() / len(y_random_prior)

mean_random_uniform_accuracy = ((y_random_uniform == np.outer(np.ones(y_random_uniform.shape[1]), y_valid).T).sum(axis=0) / len(X_valid)).mean() * 100
print('validation set, mean accuracy of uniform random guess (%):', mean_random_uniform_accuracy)
mean_random_prior_accuracy = ((y_random_prior == np.outer(np.ones(y_random_prior.shape[1]), y_valid).T).sum(axis=0) / len(X_valid)).mean() * 100
print('validation set, mean accuracy of random guess with prior class probabilities (%):', mean_random_prior_accuracy)

lr_accuracy = sum(y_pred_lr_valid == y_valid.ravel()) / len(X_valid) * 100
print('validation set, accuracy of LR model (%):', lr_accuracy)

print('the model is better than uniform random guess', lr_accuracy / mean_random_uniform_accuracy, 'times')
print('the model is better than prior random guess', lr_accuracy / mean_random_prior_accuracy, 'times')

"""
Explanation: Simple silly classification of each segment alone
End of explanation
"""


print(classification_report(y_pred_lr_valid, y_valid))

"""
Explanation: We can see that the seemingly good performance of the simple model is caused by the model trained on the skewed class prior probabilities. When we compare the model to random guessing based on the priors, the model is just 3x better than this kind of random guessing.
End of explanation
"""


def classification_results(X, y_true, y_pred):
    return pd.DataFrame(
        np.hstack([X, y_true.reshape(-1, 1), y_pred.reshape(-1, 1)]),
        columns=input_cols+['label_true']+['label_pred'])

classification_valid = classification_results(X_valid, y_valid, y_pred_lr_valid)

classification_valid['is_correct'] = classification_valid['label_true'] == classification_valid['label_pred']

print('accuracy on validation set (%)', sum(classification_valid['is_correct']) / len(classification_valid) * 100)

"""
Explanation: Notes

an accuracy of pure guess should be theoretically 1/12 ~= 8.333 %, a few rounds of random samples give a very close number
even with this simple ML method and context of a single chord the accuracy both both sets is not to bad (~ 5x better than random guess)
the training accuracy is just slightly better than the validation one
it means the model is very biased and is not overfitted, more data would likely not help
the classes are very skewed, so we could generate synthetic data points in order to make classes uniform

Let's look at the actual results of classification on the validation set.
End of explanation
"""


def key_classification_error_report(y_true, y_pred):
    errors = pd.Series(((y_true - y_pred) + 12) % 12)
    errors_hist = errors.value_counts()
    stem(errors_hist / len(errors) * 100)
    gca().set_xticks(np.arange(len(errors_hist)))
    gca().set_xticklabels(errors_hist.index)
    xlabel('error (true - predicted key pitch class)')
    ylabel('%')
    return errors, errors_hist

key_classification_error_report(classification['label_true'], classification['label_pred']);

"""
Explanation: The score() from the LogisticRegression model seems to really corresponds to manually computed accuracy.
End of explanation
"""


dataset_synth = pd.read_csv('data/beatles/derived/all_chords_with_keys_synth.tsv', sep='\t')

dataset_synth = dataset_synth.dropna()
X, y = select_X_y(dataset_synth)

len(dataset_synth), len(X), len(y)

((X_train, y_train), (X_valid, y_valid), (X_test, y_test)) = split_dataset(X, y)

lr_model_synth = LogisticRegression()
lr_model_synth.fit(X_train, y_train)

accuracy_report(lr_model_synth, X_train, y_train, X_valid, y_valid)
lr_synth_accuracy = 100 * lr_model_synth.score(X_valid, y_valid)
print('accuracy of model without synth data', lr_accuracy)
print('the model is better than uniform random guess', lr_synth_accuracy / mean_random_uniform_accuracy, 'times')

"""
Explanation: Synthetic data - all transpositions of the original dataset
This dataset show have all classes with uniform prior probabilities.
It is expected that the model would perform worse since it cannot overfit on the priors.
End of explanation
"""


y_pred_lr_synth_valid = lr_model_synth.predict(X_valid)

key_classification_error_report(classification_results(X_valid, y_valid, y_pred_lr_synth_valid));

"""
Explanation: We can se the model is worse than model trained on data with skewed classes but not dramatically.
End of explanation
"""


from sklearn.grid_search import GridSearchCV

def tune_lr_model_on_grid():
    model = LogisticRegression()
    grid = GridSearchCV(estimator=model, param_grid=dict(C=[0.1, 0.5, 1], penalty=['l2', 'l1']))
    # let the CV do multiple train-valid splits
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
    grid.fit(X_train, y_train)
    return grid

lr_model_synth_grid = tune_lr_model_on_grid()
print(lr_model_synth_grid)
print('best mean accuracy:', lr_model_synth_grid.best_score_)
print(lr_model_synth_grid.best_estimator_.C, lr_model_synth_grid.best_estimator_.penalty)
print(lr_model_synth_grid.grid_scores_)

"""
Explanation: Since there's almost no gap between training and validation accuracy the model seems to be highly biased, just too simple. We can try a few combinations of regularization parameters, but we expect regularization do not help here.
As for the kind of errors, we can see the correct key is often quite close to the true key on the circle of fifths. Moreover, the distribution of errors nicely follows the two-way distance on the circle of fifths. This is an important thing and should be taken into account when choosing/designing the cost function or validation metric.
End of explanation
"""


from sklearn.linear_model import SGDClassifier
lsgd_model_synth = SGDClassifier(loss='log')
lsgd_scores = {'train':[], 'valid':[]}
for i in range(15):
    lsgd_model_synth.partial_fit(X_train, y_train, classes=np.arange(12))
    lsgd_scores['train'].append(lsgd_model_synth.score(X_train, y_train))
    lsgd_scores['valid'].append(lsgd_model_synth.score(X_valid, y_valid))
    print('iter:', i, 'validation score:', lsgd_scores['valid'][-1])
plot(lsgd_scores['train'], label='train')
plot(lsgd_scores['valid'], label='validation')
legend();

"""
Explanation: As we can see regularization didn't help.
Linear model trained with stochastic gradient descent
End of explanation
"""


from sklearn.tree import DecisionTreeClassifier
tree_model = DecisionTreeClassifier()
tree_model = tree_model.fit(X_train, y_train)
print('training score', tree_model.score(X_train, y_train))
print('validation score', tree_model.score(X_valid, y_valid))

"""
Explanation: It seems that the logistic regression trained with SGD (full-batch) obtains a very similar result.
Decision tree
End of explanation
"""


dataset = pd.read_csv('data/beatles/derived/all_chords_with_keys_synth_rolling_16.tsv', sep='\t')
X, y = dataset.iloc[:,:-1].astype(int32), dataset.iloc[:,-1].astype(int32)
((X_train, y_train), (X_valid, y_valid), (X_test, y_test)) = split_dataset(X, y)

[d.shape for d in (X_train, y_train, X_valid, y_valid, X_test, y_test)]

lr_model = LogisticRegression()
lr_model.fit(X_train, y_train)

lr_model.score(X_train, y_train), lr_model.score(X_valid, y_valid)

key_classification_error_report(y_valid, lr_model.predict(X_valid));

"""
Explanation: Also a simple decision tree gives a similar result.
Sequence of multiple chords as an input feature
As we could see trying to classify the key based on a single chord gives not completely bad results, it is very limited. Naturally even humans perceive the key in a context of multiple chords. So we should now move to providing input features that consist of a sequence of multiple chords. We expect much better results.
In order to keep things simple we restrict our problem such that each sequence of adjacent chords lies within a single key, ie. no sequence spans multiple keys. We then classify the key for for the whole sequence.
A question is how the context size affects the accuracy. What is a good tradeoff between accuracy and redundancy?
Next steps

prepare chord blocks
put several adjacent chords within a single key into each data point
chord block are expected to carry enough information to enable more reasonable key classification
use the synthetic unskewed data
try the same model training on the chord block data
LogisticRegression [with cross-validated grid search]
try some tree-based model - eg. RandomForestClassifier
try deep neural networks
explore the effect of context size
try RNNs without the need on the full sequences (no need for explicit windowing)
End of explanation
"""


import theanets
from sklearn.metrics import accuracy_score

def theano_accuracy_report(model):
    y_train_pred = model.predict(X_train.astype(float32))
    y_valid_pred = model.predict(X_valid.astype(float32))
    print('training score:', accuracy_score(y_train, y_train_pred))
    print('validation score:', accuracy_score(y_valid, y_valid_pred))
    key_classification_error_report(y_valid, y_valid_pred)

# no hidden layers, no pretraining or regularization
exp = theanets.Experiment(theanets.Classifier, layers=[X.shape[1], 12])
exp.train(
    (X_train.astype(float32),y_train),
    (X_valid.astype(float32),y_valid),
    algo='sgd', learning_rate=1e-4, momentum=0.9)

theano_accuracy_report(exp.network)

exp.save('data/beatles/models/chords_keys_synth_theano_01.model')

# two hidden layers of size 100, autoencoer pretraining, no regularization
exp02 = theanets.Experiment(theanets.Classifier, layers=[X.shape[1], 100, 12])
exp02.train(
    (X_train.astype(float32),y_train),
    (X_valid.astype(float32),y_valid),
    algo='pretraining', learning_rate=1e-4, momentum=0.9)

theano_accuracy_report(exp02.network)

exp02.save('data/beatles/models/chords_keys_synth_theano_02.model')

# on hidden layer of size 100, dropout regularization
exp03 = theanets.Experiment(theanets.Classifier, layers=[X.shape[1], 100, 12])
exp03.train(
    (X_train.astype(float32),y_train),
    (X_valid.astype(float32),y_valid),
    algo='pretraining', learning_rate=1e-4, momentum=0.9, input_dropout=0.5)

theano_accuracy_report(exp03.network)

exp03.save('data/beatles/models/chords_keys_synth_theano_03.model')

# two hidden layers of size 100, dropout regularization
exp04 = theanets.Experiment(theanets.Classifier, layers=[X.shape[1], 100, 100, 12])
exp04.train(
    (X_train.astype(float32),y_train),
    (X_valid.astype(float32),y_valid),
    algo='pretraining', learning_rate=1e-4, momentum=0.9, input_dropout=0.5)

theano_accuracy_report(exp04.network)

exp04.save('data/beatles/models/chords_keys_synth_theano_04.model')

"""
Explanation: By adding context we can get significant improvemnt in classification performance, eg. from 38% to 71% (!) even with a simple logistic regression model.
Neural networks
End of explanation
"""


