In [55]:

    

import mir_eval
import sys
import os
import glob
import numpy as np
from pprint import pprint
import scipy.stats


    

In [56]:

    

ROOTPATH = '/home/bmcfee/git/olda/data/'


    

In [57]:

    

def load_annotations(path):
    
    files = sorted(glob.glob(path))
    
    data = [mir_eval.io.load_annotation(f) for f in files]
    
    return data


    

In [96]:

    

def evaluate_set(SETNAME, agg=True):
    
    truth = load_annotations('%s/truth/%s/*' % (ROOTPATH, SETNAME))
    
    
    algos = map(os.path.basename, sorted(glob.glob('%s/predictions/%s/*' % (ROOTPATH, SETNAME))))
    
    scores = {}
    for A in algos:
        print 'Scoring %s...' % A
        # Load the corresponding predictions
        predictions = load_annotations('%s/predictions/%s/%s/*' % (ROOTPATH, SETNAME, A))
        
        # Scrub the predictions to valid ranges
        for i in range(len(predictions)):
            predictions[i] = mir_eval.util.adjust_intervals(predictions[i][0], 
                                                            labels=predictions[i][1], 
                                                            t_max=truth[i][0][-1, -1])
            
        # Compute metrics
        my_scores = []
        
        for t, p in zip(truth, predictions):
            S = []
            S.extend(mir_eval.segment.boundary_detection(t[0], p[0], window=0.5))
            S.extend(mir_eval.segment.boundary_detection(t[0], p[0], window=3.0))
            S.extend(mir_eval.segment.boundary_deviation(t[0], p[0]))
            #S.extend(mir_eval.segment.frame_clustering_nce(t[0], t[1], p[0], p[1]))
            #S.extend(mir_eval.segment.frame_clustering_pairwise(t[0], t[1], p[0], p[1]))
            #S.extend(mir_eval.segment.frame_clustering_mi(t[0], t[1], p[0], p[1]))
            #S.append(mir_eval.segment.frame_clustering_ari(t[0], t[1], p[0], p[1]))
            my_scores.append(S)
            
        my_scores = np.array(my_scores)
        if agg:
            scores[A] = np.mean(my_scores, axis=0)
        else:
            scores[A] = my_scores
        
    return scores


    

In [59]:

    

METRICS = ['BD.5 P', 'BD.5 R', 'BD.5 F', 
           'BD3 P', 'BD3 R', 'BD3 F', 
           'BDev T2P', 'BDev P2T']
           #'S_O', 'S_U', 'S_F', 
           #'Pair_P', 'Pair_R', 'Pair_F', 
           #'MI', 'AMI', 'NMI', 'ARI']


    

In [60]:

    

def save_results(outfile, predictions):
    
    with open(outfile, 'w') as f:
        f.write('%s,%s\n' % ('Algorithm', ','.join(METRICS)))
        
        for k in predictions:
            f.write('%s,%s\n' % (k, ','.join(map(lambda x: '%.8f' % x, predictions[k]))))


    

In [61]:

    

def plot_score_histograms(data):
    
    figure(figsize=(16,10))
    for i in range(len(METRICS)):
        subplot(6,3, i+1)
        hist(data[:, i], normed=True)
        xlim([0.0, max(1.0, np.max(data[:, i]))])
        legend([METRICS[i]])


    

In [62]:

    

def plot_boxes(data):
    figure(figsize=(10,8))
    for i in range(len(METRICS)):
        subplot(6, 3, i+1)
        my_data = []
        leg = []
        for k in data:
            leg.append(k)
            my_data.append(data[k][:, i])
        my_data = np.array(my_data).T
        boxplot(my_data)
        xticks(range(1, 1+len(data)), leg)
        ylim([0, max(1.0, my_data.max())])
        tight_layout()
        title(METRICS[i])


    

In [63]:

    

def get_top_sig(SETNAME, perfs, idx, p=0.05):
    
    # Pluck out the relevant algorithm
    data = {}
    mean = {}
    best_mean = -np.inf
    best_alg  = None
    n_algs    = len(perfs)
    flip      = np.ones(len(METRICS))
    
    flip[6]   = -1  #Boundary deviation should be sign-flipped
    flip[7]   = -1
    for k in perfs:
        data[k] = flip[idx] * perfs[k][:, idx]
        mean[k] = np.mean(data[k])
        if mean[k] > best_mean:
            best_mean = mean[k]
            best_alg = k
    
    
    
    # Compute pairwise tests against the best
    sigdiff = {}
    for k in perfs:
        if k == best_alg:
            sigdiff[k] = 1.0
            continue
        # Get the p-value
        _z, _p = scipy.stats.wilcoxon(data[best_alg], data[k])
        sigdiff[k] = _p
        
        
    # Print the results
    ordering = [(v, k) for k, v in mean.iteritems()]
    ordering.sort(reverse=True)
    
    print '%s\t%s' % (METRICS[idx], SETNAME)
    
    for (v, k) in ordering:
        print '%.3f\t%10s\t%.3e\t%r' % (v, k, sigdiff[k], sigdiff[k] * (n_algs -1) < p)


    

In [64]:

    

def get_worst_examples(SETNAME, perfs, algorithm, idx, k=10):
    files = sorted(map(os.path.basename, glob.glob('%s/predictions/%s/%s/*' % (ROOTPATH, SETNAME, algorithm))))
    
    
    indices = np.argsort(perfs[algorithm][:, idx])[:k]
    
    print '%s\t%s\t%s' % (METRICS[idx], SETNAME, algorithm)
    for v in indices:
        print '%.3f\t%s' % (perfs[algorithm][v, idx], files[v])


    

In [143]:

    

ind_perfs_beatles = evaluate_set('BEATLES', agg=False)
perfs_beatles = {}
for alg in ind_perfs_beatles:
    perfs_beatles[alg] = np.mean(ind_perfs_beatles[alg], axis=0)


    

    




Scoring fda...
Scoring olda...
Scoring rfda...
Scoring smga1...
Scoring unsup...

In [144]:

    

pprint(perfs_beatles)


    

    




{'fda': array([ 0.30410974,  0.30993012,  0.29630028,  0.55470147,  0.53647853,
        0.52978163,  3.02629609,  3.24893575]),
 'olda': array([ 0.29742039,  0.32024792,  0.29628414,  0.55337542,  0.55049432,
        0.5348491 ,  2.8535838 ,  3.2883743 ]),
 'rfda': array([ 0.31366817,  0.3350548 ,  0.31081947,  0.56989802,  0.56460779,
        0.5493749 ,  2.68242737,  3.13565642]),
 'smga1': array([ 0.14448849,  0.16943113,  0.15311947,  0.62092775,  0.72784284,
        0.65786475,  1.87511341,  3.19086285]),
 'unsup': array([ 0.27758455,  0.27093257,  0.26626857,  0.56002263,  0.50990414,
        0.52206759,  3.46052514,  3.26752235])}

In [150]:

    

save_results('/home/bmcfee/git/olda/data/final_beatles_scores.csv', perfs_beatles)


    

In [151]:

    

del ind_perfs_beatles['rfda']


    

    




---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-151-ea0bb984a61d> in <module>()
----> 1 del ind_perfs_beatles['rfda']

KeyError: 'rfda'

In [152]:

    

for idx in range(len(METRICS)):
    get_top_sig('BEATLES', ind_perfs_beatles, idx=idx)
    print


    

    




BD.5 P	BEATLES
0.304	       fda	1.000e+00	False
0.297	      olda	4.722e-01	False
0.278	     unsup	2.742e-02	False
0.144	     smga1	4.234e-18	True

BD.5 R	BEATLES
0.320	      olda	1.000e+00	False
0.310	       fda	1.584e-01	False
0.271	     unsup	1.512e-04	True
0.169	     smga1	1.046e-12	True

BD.5 F	BEATLES
0.296	       fda	1.000e+00	False
0.296	      olda	9.778e-01	False
0.266	     unsup	9.898e-03	True
0.153	     smga1	1.464e-15	True

BD3 P	BEATLES
0.621	     smga1	1.000e+00	False
0.560	     unsup	5.518e-04	True
0.555	       fda	1.083e-03	True
0.553	      olda	1.825e-04	True

BD3 R	BEATLES
0.728	     smga1	1.000e+00	False
0.550	      olda	4.589e-13	True
0.536	       fda	5.505e-15	True
0.510	     unsup	1.039e-17	True

BD3 F	BEATLES
0.658	     smga1	1.000e+00	False
0.535	      olda	5.482e-10	True
0.530	       fda	2.041e-10	True
0.522	     unsup	8.033e-12	True

BDev T2P	BEATLES
-1.875	     smga1	1.000e+00	False
-2.854	      olda	8.630e-07	True
-3.026	       fda	3.258e-08	True
-3.461	     unsup	2.797e-12	True

BDev P2T	BEATLES
-3.191	     smga1	1.000e+00	False
-3.249	       fda	5.941e-01	False
-3.268	     unsup	3.352e-01	False
-3.288	      olda	3.994e-01	False

In [153]:

    

plot_boxes(ind_perfs_beatles)


    

In [154]:

    

for alg in sorted(ind_perfs_beatles.keys()):
    get_worst_examples('BEATLES', ind_perfs_beatles, alg, 2, 10)
    print


    

    




BD.5 F	BEATLES	fda
0.000	13_-_She_Came_In_Through_The_Bathroom_Window.lab
0.000	01_-_Come_Together.lab
0.000	CD2_-_06_-_Helter_Skelter.lab
0.000	09_-_You_Never_Give_Me_Your_Money.lab
0.000	CD1_-_02_-_Dear_Prudence.lab
0.000	CD1_-_06_-_The_Continuing_Story_of_Bungalow_Bill.lab
0.000	CD1_-_14_-_Don't_Pass_Me_By.lab
0.000	05_-_Dig_It.lab
0.087	CD2_-_07_-_Long_Long_Long.lab
0.091	01_-_Two_of_Us.lab

BD.5 F	BEATLES	olda
0.000	12_-_I_Want_To_Tell_You.lab
0.000	05_-_Dig_It.lab
0.000	CD1_-_06_-_The_Continuing_Story_of_Bungalow_Bill.lab
0.000	CD1_-_02_-_Dear_Prudence.lab
0.000	CD2_-_06_-_Helter_Skelter.lab
0.083	CD2_-_07_-_Long_Long_Long.lab
0.087	04_-_Love_You_To.lab
0.091	01_-_Two_of_Us.lab
0.091	01_-_Taxman.lab
0.095	12_-_Get_Back.lab

BD.5 F	BEATLES	smga1
0.000	CD1_-_02_-_Dear_Prudence.lab
0.000	11_-_Mean_Mr_Mustard.lab
0.000	11_-_I_Wanna_Be_Your_Man.lab
0.000	11_-_Good_Morning_Good_Morning.lab
0.000	04_-_Chains.lab
0.000	07_-_Maggie_Mae.lab
0.000	10_-_Baby_It's_You.lab
0.000	04_-_I_Me_Mine.lab
0.000	04_-_I_Need_You.lab
0.000	04_-_Love_You_To.lab

BD.5 F	BEATLES	unsup
0.000	09_-_Girl.lab
0.000	16_-_The_End.lab
0.000	CD1_-_06_-_The_Continuing_Story_of_Bungalow_Bill.lab
0.000	05_-_Dig_It.lab
0.000	06_-_The_Word.lab
0.000	02_-_With_A_Little_Help_From_My_Friends.lab
0.000	13_-_She_Came_In_Through_The_Bathroom_Window.lab
0.000	08_-_Eight_Days_a_Week.lab
0.000	CD2_-_05_-_Sexy_Sadie.lab
0.000	12_-_Get_Back.lab

In [161]:

    

ind_perfs_salami = evaluate_set('SALAMI', agg=False)
perfs_salami = {}
for alg in ind_perfs_salami:
    perfs_salami[alg] = np.mean(ind_perfs_salami[alg], axis=0)


    

    




Scoring CNMF...
Scoring SIPLCA...
Scoring fda...
Scoring olda...
Scoring rfda...
Scoring smga1...
Scoring unsup...

In [162]:

    

pprint(perfs_salami)


    

    




{'CNMF': array([ 0.10514414,  0.13322015,  0.10971976,  0.45024537,  0.54333524,
        0.46309017,  4.1670073 ,  5.94623858]),
 'SIPLCA': array([  0.2100008 ,   0.10169394,   0.12774307,   0.45112871,
         0.22782295,   0.28610934,  16.93411006,   5.80685512]),
 'fda': array([ 0.29030651,  0.22163841,  0.24333276,  0.55213126,  0.42938182,
        0.46630035,  7.00804464,  4.4818246 ]),
 'olda': array([ 0.26347095,  0.24773284,  0.24690428,  0.50549125,  0.48305071,
        0.47662926,  4.42681576,  5.28404255]),
 'rfda': array([ 0.25941641,  0.23446233,  0.2386172 ,  0.50677901,  0.46951904,
        0.4706314 ,  4.64858833,  5.15115299]),
 'smga1': array([ 0.11493075,  0.17796945,  0.13437579,  0.43375133,  0.6661867 ,
        0.5076943 ,  2.09821367,  7.15111203]),
 'unsup': array([ 0.23306619,  0.18794069,  0.20194842,  0.49217135,  0.4059798 ,
        0.43095289,  6.39818556,  5.87488361])}

In [157]:

    

save_results('/home/bmcfee/git/olda/data/final_salami_scores.csv', perfs_salami)


    

In [158]:

    

for alg in sorted(ind_perfs_salami.keys()):
    get_worst_examples('SALAMI', ind_perfs_salami, alg, 2, 5)
    print


    

    




BD.5 F	SALAMI	CNMF
0.000	1266.lab
0.000	1074.lab
0.000	1392.lab
0.000	1078.lab
0.000	1390.lab

BD.5 F	SALAMI	SIPLCA
0.000	1352.lab
0.000	1386.lab
0.000	1078.lab
0.000	1470.lab
0.000	1230.lab

BD.5 F	SALAMI	fda
0.000	1000.lab
0.000	976.lab
0.000	970.lab
0.000	1498.lab
0.000	1392.lab

BD.5 F	SALAMI	olda
0.000	1334.lab
0.000	1362.lab
0.000	976.lab
0.000	970.lab
0.000	1496.lab

BD.5 F	SALAMI	rfda
0.000	1226.lab
0.000	1362.lab
0.000	1366.lab
0.000	1082.lab
0.000	1376.lab

BD.5 F	SALAMI	smga1
0.000	1338.lab
0.000	1478.lab
0.000	1050.lab
0.000	1170.lab
0.000	1474.lab

BD.5 F	SALAMI	unsup
0.000	1000.lab
0.000	1142.lab
0.000	1162.lab
0.000	1170.lab
0.000	1190.lab

In [165]:

    

del ind_perfs_salami['rfda']


    

In [166]:

    

for idx in range(len(METRICS)):
    get_top_sig('SALAMI', ind_perfs_salami, idx=idx)
    print


    

    




BD.5 P	SALAMI
0.290	       fda	1.000e+00	False
0.263	      olda	6.203e-05	True
0.233	     unsup	1.344e-12	True
0.210	    SIPLCA	1.296e-08	True
0.115	     smga1	8.342e-36	True
0.105	      CNMF	4.832e-35	True

BD.5 R	SALAMI
0.248	      olda	1.000e+00	False
0.222	       fda	7.092e-05	True
0.188	     unsup	1.754e-17	True
0.178	     smga1	9.086e-09	True
0.133	      CNMF	2.307e-18	True
0.102	    SIPLCA	3.383e-27	True

BD.5 F	SALAMI
0.247	      olda	1.000e+00	False
0.243	       fda	9.357e-01	False
0.202	     unsup	1.686e-12	True
0.134	     smga1	4.908e-24	True
0.128	    SIPLCA	1.595e-21	True
0.110	      CNMF	1.022e-27	True

BD3 P	SALAMI
0.552	       fda	1.000e+00	False
0.505	      olda	7.965e-09	True
0.492	     unsup	4.723e-11	True
0.451	    SIPLCA	1.677e-09	True
0.450	      CNMF	5.195e-12	True
0.434	     smga1	3.439e-19	True

BD3 R	SALAMI
0.666	     smga1	1.000e+00	False
0.543	      CNMF	1.977e-14	True
0.483	      olda	7.783e-25	True
0.429	       fda	2.623e-32	True
0.406	     unsup	4.160e-35	True
0.228	    SIPLCA	6.485e-42	True

BD3 F	SALAMI
0.508	     smga1	1.000e+00	False
0.477	      olda	3.872e-03	True
0.466	       fda	4.921e-04	True
0.463	      CNMF	4.183e-05	True
0.431	     unsup	3.691e-11	True
0.286	    SIPLCA	5.717e-34	True

BDev T2P	SALAMI
-2.098	     smga1	1.000e+00	False
-4.167	      CNMF	1.150e-20	True
-4.427	      olda	1.524e-21	True
-6.398	     unsup	1.375e-35	True
-7.008	       fda	4.007e-34	True
-16.934	    SIPLCA	3.570e-43	True

BDev P2T	SALAMI
-4.482	       fda	1.000e+00	False
-5.284	      olda	1.323e-06	True
-5.807	    SIPLCA	4.955e-09	True
-5.875	     unsup	2.557e-11	True
-5.946	      CNMF	8.853e-13	True
-7.151	     smga1	8.048e-21	True

In [167]:

    

plot_boxes(ind_perfs_salami)


    

In [124]:

    

import librosa
import scipy.signal
import functools


    

In [125]:

    

def get_beat_mfccs(filename):
    y, sr = librosa.load(filename)
    
    S = librosa.feature.melspectrogram(y, sr, n_fft=2048, hop_length=64, n_mels=128, fmax=8000)
    
    tempo, beats = librosa.beat.beat_track(y=y, sr=sr, hop_length=64)
    
    M = librosa.feature.mfcc(librosa.logamplitude(S), n_mfcc=32)
    M = librosa.feature.sync(M, beats)
    return M


    

In [141]:

    

def compress_data(X, k):
    e_vals, e_vecs = scipy.linalg.eig(X.dot(X.T))
    #e_vals, e_vecs = scipy.linalg.eig(np.cov(X))
        
    e_vals = np.maximum(0.0, np.real(e_vals))
    e_vecs = np.real(e_vecs)
        
    idx = np.argsort(e_vals)[::-1]
        
    e_vals = e_vals[idx]
    e_vecs = e_vecs[:, idx]
        
    # Truncate to k dimensions
    if k < len(e_vals):
        e_vals = e_vals[:k]
        e_vecs = e_vecs[:, :k]
        
    # Normalize by the leading singular value of X
    Z = np.sqrt(e_vals.max())
        
    if Z > 0:
        e_vecs = e_vecs / Z
        
    return e_vecs.T.dot(X)


    

In [139]:

    

def make_rep_feature_plot(M):
    
    #R = librosa.segment.recurrence_matrix(librosa.segment.stack_memory(M), metric='seuclidean')
    R = librosa.segment.recurrence_matrix(M, metric='seuclidean')
    
    Rskew = librosa.segment.structure_feature(R)
    Rskew = np.roll(Rskew, M.shape[1], axis=0)
    
    
    Rfilt = scipy.signal.medfilt2d(Rskew.astype(np.float32), kernel_size=(1, 7))
    #Rfilt = Rfilt[Rfilt.sum(axis=1) > 0, :]
    
    Rlatent = compress_data(Rfilt, 8)
    #Rlatent = compress_data(Rfilt, R.shape[0])
    
    figure(figsize=(6,6))
    subplot(221)
    librosa.display.specshow(R), title('Self-similarity')
    xlabel('Beat'), ylabel('Beat')
    xticks(range(0, M.shape[1] + 1, M.shape[1] / 6))
    yticks(range(0, M.shape[1] + 1, M.shape[1] / 6))
    
    subplot(222)
    librosa.display.specshow(Rskew, cmap='gray_r'), title('Skewed self-sim.')
    xlabel('Beat'), ylabel('Lag')
    yticks(range(0, Rskew.shape[0] + 1, Rskew.shape[0] / 6), range(-M.shape[1]+2, 1+M.shape[1], Rskew.shape[0] / 6))
    xticks(range(0, M.shape[1] + 1, M.shape[1] / 6))
    
    subplot(223)
    librosa.display.specshow(Rfilt, cmap='gray_r'), title('Filtered self-sim.')
    xlabel('Beat'), ylabel('Lag')
    yticks(range(0, Rskew.shape[0] + 1, Rskew.shape[0] / 6), range(-M.shape[1]+2, 1+M.shape[1], Rskew.shape[0] / 6))
    xticks(range(0, M.shape[1] + 1, M.shape[1] / 6))
    
    subplot(224)
    librosa.display.specshow(Rlatent, origin='upper'), title('Latent repetition')
    xticks(range(0, M.shape[1] + 1, M.shape[1] / 6))
    xlabel('Beat'), ylabel('Factor'), yticks(range(Rlatent.shape[0]))
    tight_layout()
    
    #savefig('/home/bmcfee/git/olda/paper/figs/rep.pdf', format='pdf', pad_inches=0, transparent=True)
    #savefig('/home/bmcfee/git/olda/paper/figs/rep.svg', format='svg', pad_inches=0, transparent=True, dpi=200)


    

In [128]:

    

M = get_beat_mfccs('/home/bmcfee/data/CAL500/mp3/2pac-trapped.mp3')


    

In [142]:

    

make_rep_feature_plot(M)


    

In [131]:

    

make_rep_feature_plot(M[:,40:137])


    

In [69]:

    

import librosa


    

In [116]:

    

model_olda  = np.load('/home/bmcfee/git/olda/data/model_olda_beatles.npy')
figure(figsize=(8,5))
librosa.display.specshow(model_olda, origin='upper')
colorbar()
yticks([])
xticks([0, 32, 44, 76, 108], ['MFCC', 'Chroma', 'R-MFCC', 'R-Chroma', 'Time'], horizontalalignment='left')

tight_layout()
#savefig('/home/bmcfee/git/olda/model_olda_salami_w.png', format='png', pad_inches=0, transparent=True)


    

In [117]:

    

librosa.display.specshow(model_olda.T.dot(model_olda), origin='upper')


    

    
Out[117]:





<matplotlib.image.AxesImage at 0x5968d10>






    

In [118]:

    

v, e = scipy.linalg.eig(model_olda.dot(model_olda.T))

plot(np.sort(v)[::-1])
axis('tight')


    

    
Out[118]:





(0.0, 111.0, 0.74638458696734089, 12.196980612549629)






    

In [119]:

    

model_fda = np.load('/home/bmcfee/git/olda/data/model_fda_salami.npy')
model_rfda = np.load('/home/bmcfee/git/olda/data/model_rfda_salami.npy')
model_olda  = np.load('/home/bmcfee/git/olda/data/model_olda_salami.npy')
figure(figsize=(8,10))
subplot(311)
librosa.display.specshow(model_fda, origin='upper')

ylabel('SALAMI: FDA')
yticks([])
#ylabel("More important $\\rightarrow$")
#xticks([0, 32, 44, 76, 108], ['MFCC', 'Chroma', 'Rep-M', 'Rep-C', 'Time'], rotation=-30, horizontalalignment='left')
xticks([])
#colorbar()

subplot(312)
librosa.display.specshow(model_rfda, origin='upper')
ylabel('SALAMI: RFDA')
#colorbar(orientation='horizontal')
#ylabel("More important $\\rightarrow$")
yticks([])


subplot(313)
librosa.display.specshow(model_olda, origin='upper')
ylabel('SALAMI: OLDA')
colorbar(orientation='horizontal', use_gridspec=True )
#ylabel("More important $\\rightarrow$")
yticks([])
xticks([0, 32, 44, 76, 108], ['MFCC', '$\uparrow$\nChroma', 'R-MFCC', 'R-Chroma', 'Time'], horizontalalignment='left')

tight_layout()
#savefig('/home/bmcfee/git/olda/paper/figs/w.pdf', format='pdf', pad_inches=0, transparent=True)


    

In [120]:

    

model_fda_beatles = np.load('/home/bmcfee/git/olda/data/model_fda_beatles.npy')
model_fda_salami = np.load('/home/bmcfee/git/olda/data/model_fda_salami.npy')
model_olda_beatles  = np.load('/home/bmcfee/git/olda/data/model_olda_beatles.npy')
model_olda_salami = np.load('/home/bmcfee/git/olda/data/model_olda_salami.npy')

figure(figsize=(14,8))
subplot(221)
imshow(model_fda_beatles, aspect='auto', interpolation='none', cmap='PRGn_r')
ylabel('FDA - Beatles')
yticks([])
xticks([])

subplot(222)
imshow(model_fda_salami, aspect='auto', interpolation='none', cmap='PRGn_r')
ylabel('FDA - Salami')
yticks([])
xticks([])

subplot(223)
imshow(model_olda_beatles, aspect='auto', interpolation='none', cmap='PRGn_r')
ylabel('OLDA - Beatles')
yticks([])
xticks([0, 32, 44, 76, 108], ['MFCC', '$\uparrow$\nChroma', 'R-MFCC', 'R-Chroma', 'Time'], horizontalalignment='left')

subplot(224)
imshow(model_olda_salami, aspect='auto', interpolation='none', cmap='PRGn_r')
ylabel('OLDA - Salami')
yticks([])
xticks([])
xticks([0, 32, 44, 76, 108], ['MFCC', '$\uparrow$\nChroma', 'R-MFCC', 'R-Chroma', 'Time'], horizontalalignment='left')

tight_layout()
#savefig('/home/bmcfee/git/olda/paper/figs/fda-vs-olda.pdf', format='pdf', pad_inches=0, transparent=True)


    

In [132]:

    

def rep_feature_svd(M):
    
    myimshow = functools.partial(imshow, aspect='auto', interpolation='none', origin='lower', cmap='gray_r')
    
    R = librosa.segment.recurrence_matrix(M, k=2*np.sqrt(1./M.shape[1]), sym=False).astype(np.float)
    
    Rskew = librosa.segment.structure_feature(R, pad=True)
    Rskew = np.roll(Rskew, M.shape[1], axis=0)
    
    
    Rfilt = scipy.signal.medfilt2d(Rskew, kernel_size=(1, 7))
    #Rfilt = Rfilt[Rfilt.sum(axis=1) > 0, :]
    
    Rlatent = compress_data(Rfilt, 8)
    #Rlatent = compress_data(Rfilt, R.shape[0])
    
    figure(figsize=(12,4))
    subplot(131)
    myimshow(R), title('Self-similarity')
    xlabel('Beat'), ylabel('Beat')
    xticks(range(0, M.shape[1] + 1, M.shape[1] / 6))
    yticks(range(0, M.shape[1] + 1, M.shape[1] / 6))
    
    subplot(132)
    myimshow(Rskew), title('Skewed self-sim.')
    xlabel('Beat'), ylabel('Lag')
    yticks(range(0, Rskew.shape[0] + 1, Rskew.shape[0] / 6), range(-M.shape[1]+1, M.shape[1], Rskew.shape[0] / 6))
    xticks(range(0, M.shape[1] + 1, M.shape[1] / 6))
    
    subplot(133)
    myimshow(Rfilt), title('Filtered self-sim.')
    xlabel('Beat'), ylabel('Lag')
    yticks(range(0, Rskew.shape[0] + 1, Rskew.shape[0] / 6), range(-M.shape[1]+1, M.shape[1], Rskew.shape[0] / 6))
    xticks(range(0, M.shape[1] + 1, M.shape[1] / 6))
    tight_layout()
    
    # Do the SVD
    U, sigma, Vh = scipy.linalg.svd(Rfilt)
    
    D_little = 32
    
    figure(figsize=(12,4))
    subplot(131)
    myimshow(U[:,:D_little], cmap='PRGn'), title('U')
    xlabel('Factor'), ylabel('Lag')
    #yticks(range(0, Rskew.shape[0] + 1, Rskew.shape[0] / 6), range(-M.shape[1]+1, M.shape[1], Rskew.shape[0] / 6))
    #xticks(range(0, M.shape[1] + 1, M.shape[1] / 6))
    
    subplot(132)
    plot(sigma / sigma[0]), axis('tight'), title('Normalized spectrum $\sigma/\sigma_1$')
    vlines([D_little], 0, 1)
    
    subplot(133)
    myimshow(Vh[:D_little], cmap='PRGn'), title('V\'')
    xlabel('Beat'), ylabel('Factor')
    #yticks(range(0, Rskew.shape[0] + 1, Rskew.shape[0] / 6), range(-M.shape[1]+1, M.shape[1], Rskew.shape[0] / 6))
    #xticks(range(0, M.shape[1] + 1, M.shape[1] / 6))
    tight_layout()
    # Reconstruct
    S_hat = np.zeros(U.shape[1])
    S_hat[:D_little] = sigma[:D_little]
    
    R_reconst = U.dot(np.diag(S_hat))
    R_reconst = R_reconst[:, :Vh.shape[0]].dot(Vh)
    
    figure(figsize=(12,6))
    subplot(121)
    myimshow(Rfilt), title('Original')
    subplot(122)
    myimshow(R_reconst), title('Reconstruction (d=%2d)' % D_little)
    tight_layout()
    
    figure(figsize=(16,32))
    
    for i in range(D_little):
        subplot(1+np.ceil(np.sqrt(D_little)), np.floor(np.sqrt(D_little)), i+1)
        myimshow(np.outer(U[:,i], Vh[i])), title('$U_{%d} V^\mathsf{T}_{%d}$' % (i, i))
    
    tight_layout()


    

In [134]:

    

rep_feature_svd(M[:,40:137])