annotate dml-cla/python/similarity.py @ 0:718306e29690 tip

commiting public release

# Part of DML (Digital Music Laboratory)
# Copyright 2014-2015 Daniel Wolff, City University
 
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public
# License along with this library; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA

# -*- coding: utf-8 -*-
__author__='wolffd'

# this script derives all pairwise similarity measures for the chroma vectors provided.
# as a first experiment, only the mean chroma vectors per piece are compared
# using euclidean distance


# parameters to be forwarded to API

# similarity type:
#   euclidean, compression
# simtype = 'euclidean'

# parallelisation 
num_cores = 10


#min_clusters = 40## unused
#max_clusters = 256## unused

#set_clusters = 40
#max_clips = 50
encoding = 'binary'
#compressor = 'zxd'
mds_init_tries = 4
mds_max_iter = 100
mfccbins = 12

# resample chroma / timbre values at this fraction for compression distance. 0 to switch off
# we want a vector just every second.
# The standard sample rate and window size are 44100 / 1024 for chroma / timbre
# this is dependent on the sim_downsample parameter
resample_factor = 44100/1024;


from rdflib import RDF, RDFS
from csvutils import *
from aggregate import *
from n3Parser import get_rdf_graph_from_n3 

# numpy, scipy
import numpy as np
from scipy.spatial import distance
from sklearn.metrics.pairwise import pairwise_distances
from scipy.signal import resample

#scikitlearn
from sklearn.datasets import make_blobs
from sklearn.cluster import KMeans
# from sklearn.metrics import silhouette_samples, silhouette_score
from sklearn import manifold


# chord processing
from chord_seq_key_relative import chords_from_csv, keys_from_csv, chord2function, fun2txt, fun2num,most_frequent_key,chord_roots,type_labels

# subprocess, command line and threading
import os, tempfile
import subprocess, threading

# for system / compression calls
import zlib


def chroma_from_csv(filename):
    # we assume CSV: time, chroma 1, ...  chroma 12]
    # return (time, [chroma 1-12])
    return csv_map_rows(filename,13, lambda row:(float(row[0]),np.array(row[1:12],dtype=float)))
    
def mfcc_from_csv(filename):
    # we assume CSV: time, mfcc 1, ...  mfcc 20]
    # return (time, [chroma 1-12])
    return csv_map_rows(filename,21, lambda row:(float(row[0]),np.array(row[1:20],dtype=float)))

    
chroma_parser_table = { 'csv':chroma_from_csv}
mfcc_parser_table = { 'csv':mfcc_from_csv}

# in chord_seq_relative
key_parser_table = { 'csv':keys_from_csv }
chord_parser_table = { 'csv':chords_from_csv }

## generate global dict of chord_keys
chord_keys = []
for chordnum in range(1,12+1):
    for typenum in range(1,11+1):
        chord_keys.append("%02d%02d" % (chordnum,typenum))

def per_file(inputs,opts={}):
    chromas = []
    chromas_idx = []
    mfccs = []
    mfccs_idx = []
    chords = []
    chords_idx = []
    uris = []

    # get options from API
    # print_status(str(opts))
    simtype = opts['sim_type']
    set_clusters = opts['sim_clusters'] # def 40
    downsample = opts['sim_downsample'] # def 1
    limit = opts['sim_reclimit'] # def 50
    compressor = opts['sim_compressor'] # def 'zlib'
    
    # parse feature list
    features = opts['sim_features'].split(',') # features, def: chroma
    use_chromagram = any(ext in 'chromagram' for ext in features)
    use_mfcc = any(ext in 'mfcc' for ext in features)
    use_chords = any(ext in 'chords' for ext in features)
    
    # check number of inputs
    if len(inputs) > limit:
        #return { 'error': ''}
        print_status('Similarity: Too many inputs, truncating collection')
        inputs = inputs[0:limit]
        

    # accumulation for euclidean just gets the mean values over the whole clips
    # todo: add std and other statistics?
    def accum_euclidean(item):
        # accumulate chroma vectors for this piece
        if use_chromagram:
            chroma = [ res[1] for res in decode_tagged(chroma_parser_table,item['chromagram'])]
            # print_status('Chroma Raw Data' + str(chroma))
            # get mean chroma vector
            chroma_mean = np.mean(np.array(chroma), axis = 0)
            #print_status('Chroma Means' + str(chroma_mean))

            # add vector to chromas table
            chromas.append(chroma_mean) 
            
        if use_mfcc:
            mfcc = [ res[1] for res in decode_tagged(mfcc_parser_table,item['mfcc'])]
            mfcc_mean = np.mean(np.array(mfcc), axis = 0)
            mfccs.append(mfcc_mean) 
            
        if use_chords:
            # get duration and normalised frequency for all tuning pitches (A3,A4,A5)
            keys = decode_tagged(key_parser_table,item['keys'])
            # get most frequent key
            key,mode = most_frequent_key(keys)
            relchords = []      
            for (time,chord) in  decode_tagged(chord_parser_table,item['chords']):

                # get chord function
                (root,fun,typ, bfun) = chord2function(chord, key,mode)

                # translate into text
                #txt = fun2txt(fun,typ, bfun, mode)
                #print_status('Chord: ' + chord + ', function: ' + txt)
                num = fun2num(fun,typ, bfun, mode)
                if num > 0:
                    # add to chords of this clip
                    #relchords.append((time,key,mode,fun,typ,bfun))
                    
                    # get the root note of the chord and chord type
                    # ignore mode and base note
                    # format of num [1x mode, 2x function, 2x type, 2x base note]
                    relchords.append(str(num)[1:5])
                    
            # append histogram of all chords for this recording   
            hist = chord_histogram(relchords)
            chords.append(hist)
         
        
        # add uri if everything went well
        uris.append(item['list'])
        
    # accumulation for compression: 
    # save all chroma vectors
    # possibly build a codebook
    # otherwise compare for quantisation
    
    def accum_compression(item):
    
        # get chromas
        if use_chromagram:
            # accumulate chroma vectors for this piece
            chroma = [ res[1] for res in decode_tagged(chroma_parser_table,item['chromagram'])]
            # print_status('Chroma Raw Data' + str(chroma))
            
            # downsample if necessary
            if downsample == 1:
                #chroma = resample(chroma, len(chroma)//resample_factor, axis=0, window=None)
                #chroma = [chroma[i] for i in np.random.randint(0,len(chroma),len(chroma)//resample_factor)]
                chroma = [chroma[i*resample_factor] for i in range(0,len(chroma)//resample_factor)]
            
            chromas.extend(chroma)
            chromas_idx.append(len(chromas))
            
                
        if use_mfcc:
            mfcc = [ res[1] for res in decode_tagged(mfcc_parser_table,item['mfcc'])]
           
            if downsample == 1:
                # mfcc = np.random.randint(0,len(mfcc),len(mfcc)//resample_factor)]
                mfcc = [mfcc[i*resample_factor] for i in range(0,len(mfcc)//resample_factor)]
            mfccs.extend(mfcc)
            mfccs_idx.append(len(mfccs))
            
        if use_chords:
            # get duration and normalised frequency for all tuning pitches (A3,A4,A5)
            keys = decode_tagged(key_parser_table,item['keys'])
            # get most frequent key
            key,mode = most_frequent_key(keys)
            relchords = []      
            for (time,chord) in  decode_tagged(chord_parser_table,item['chords']):

                # get chord function
                (root,fun,typ, bfun) = chord2function(chord, key,mode)

                # translate into text
                #txt = fun2txt(fun,typ, bfun, mode)
                #print_status('Chord: ' + chord + ', function: ' + txt)
                num = fun2num(fun,typ, bfun, mode)
                if num > 0:
                    # add to chords of this clip
                    #relchords.append((time,key,mode,fun,typ,bfun))
                    
                    # get the root note of the chord and chord type
                    # ignore mode and base note
                    # format of num [1x mode, 2x function, 2x type, 2x base note]
                    relchords.append(int(str(num)[1:5]))
                    
            # append histogram of all chords for this recording   
            #hist = chord_histogram(relchords)
            chords.extend(relchords)
            chords_idx.append(len(chords))
            
        # add uri if everything went well
        uris.append(item['list'])
    
            
    
    # ---
    # this is the euclidean distance
    # ---
    if (simtype == 'euclidean'):
        # accumulate over all inputs        
        st=for_each(inputs,accum_euclidean)
    
        # concatenate feature input for all features
        arr = np.empty((len(uris),0), float)
        # concatenate data to nparray for euclidean distance
        if use_chromagram:
            arr = np.append(arr, np.array(chromas), axis=1)

        if use_mfcc:
            arr = np.append(arr, np.array(mfccs), axis=1)
            
        if use_chords:
            # get chord dictionaries
            #print(str(np.array(chords).shape))
            arr = np.append(arr,np.array(chords) , axis=1)

        #dist = distance.pdist(chromas, 'euclidean')
        dist = pairwise_distances(arr, metric = 'euclidean', n_jobs = num_cores)

        # return to non-condensed matrix for simplicity.
        # this can be reversed using the very same function for data
        # efficiency
        #dist = distance.squareform(dist)
        
    # ---
    # this is the normalised compression distance
    # ---
    elif (simtype == 'compression'):
        # accumulate over all inputs    
        print_status('Similarity Module: Accumulating')
        st=for_each(inputs,accum_compression)
        
        dist = np.zeros((len(uris),len(uris)))
        count = 0
        if use_chromagram:
            print_status('Similarity Module: Chroma Quantisation')
            chromas_coded = vector_quantisation(np.array(chromas), set_clusters,num_cores)
            print_status('Similarity Module: Chroma Compression Results')
            dist += similarity_by_mask(chromas_coded,chromas_idx,compressor,encoding)
            count +=1
            
        if use_mfcc:
            print_status('Similarity Module: MFCC Quantisation')
            mfccs_coded = vector_quantisation(np.array(mfccs), set_clusters,num_cores)
            print_status('Similarity Module: MFCC Compression Results')
            dist += similarity_by_mask(mfccs_coded,mfccs_idx,compressor,encoding)
            count +=1
        
        if use_chords:
            print_status('Similarity Module: Chord Compression Results')
            dist += similarity_by_mask(np.array(chords),chords_idx,compressor,encoding)
            count +=1
            
        dist = dist / count

        
        # get rid of zeros in between            
        #for idx1 in range(0,len(chromas_idx)): 
        #    dist[idx1][idx1] = 1
            
        print_status('dist' + str(dist))
        
    # Do MDS scaling with precomputed distance
    mds = manifold.MDS(n_components = 2, max_iter=mds_max_iter, n_init=mds_init_tries, dissimilarity='precomputed')
    
    coordinates = mds.fit_transform(dist)
                    
    return { 'result': { 'list': uris, 'mds': coordinates.tolist()}, 'stats' : st }
    # return { 'result': { 'list': uris, 'distance': dist.tolist(), 'mds': coordinates.tolist()}, 
    #          'stats' : st }
             
             
             
def vector_quantisation(data, set_clusters,num_cores):
    # ---
    # build codebook!    
    # ---
    # --- 1 quantise chroma data
    # --- 1a use scikit-learn k-means
    # http://scikit-learn.org/stable/modules/clustering.html
    
    # quick try
    clusterer = KMeans(n_clusters=set_clusters,n_jobs = num_cores)
    
    # --- 2 get compression distance
    # get all single compressed sizes? 
    data_coded = clusterer.fit_predict(data)
    #print_status('Chromas Coded' + str(chromas_coded))
    # print_status('Coding Histogram' + str(np.histogram(chromas_coded)))
    return data_coded 



def similarity_by_mask(data_coded,data_idx,compressor,encoding):

    # idx is expected to start with the first chroma index of the second piece
    # TODO: check indexing (starts at 0 or 1?)
    lengths = []
    start_idx = [0] + data_idx[:-1]
    dist = np.zeros((len(data_idx),len(data_idx)))
    for idx1 in range(0,len(data_idx)):     
        for idx2 in range(0,len(data_idx)):
            if (idx2 < idx1):
                # select encoded chromas for the clips
                data1_mask = np.zeros(len(data_coded), dtype=bool) 
                data1_mask[start_idx[idx1]:data_idx[idx1]-1] = True  
                
                data2_mask = np.zeros(len(data_coded), dtype=bool) 
                data2_mask[start_idx[idx2]:data_idx[idx2]-1] = True

                a_coded = encode(data_coded[data1_mask],format = encoding)
                b_coded = encode(data_coded[data2_mask],format = encoding)
                # get compression lengths
                if compressor == 'zlib':
                    (a,b,ab) = compressed_length(a_coded,b_coded,compressor)
                    
                else:                 
                    # get complement chroma set
                    ref_mask = ~data1_mask & ~data2_mask
                    ref_coded = encode(data_coded[ref_mask],format = encoding)
                    (a,b,ab) = delta_compressed_length(a_coded,b_coded,ref_coded,compressor)
                
                #NCD(z - min(v, w))/ max(v, w);
                dist[idx1][idx2] = (ab - min(a,b))/float(max(a,b))
                
                # the above normalised compression distance is symmetric
                # this is required by the nds routine below
                dist[idx2][idx1] = dist[idx1][idx2]
                
    return dist
             
def encode(data, format = 'string'):

    # Encoding
    if format == 'binary':
        data_coded = data.tostring()
    elif format == 'string':
        data_coded = str(data)
    return data_coded       
     
def compressed_length(a_coded,b_coded, type = 'zlib'):

    # Compression
    if type == 'zlib':
        # zlib is quite helpful https://docs.python.org/2/library/zlib.html#module-zlib
        a = len(zlib.compress(a_coded, 9))
        b = len(zlib.compress(a_coded, 9))
        ab = len(zlib.compress(a_coded + b_coded, 9))
    
    return (a,b,ab)
    
def delta_compressed_length(a_coded,b_coded,ref_coded, type = 'zxd'):
    # Compression
    # zbs - use bsdiff
    # zxd - uses xdelta3
    # zvcd - uses open-vcdiff
    # zvcz - uses vczip
    # zdiff - converts binary to text and uses diff to produce an ed script

    if type == 'zxd' or type == 'zbs' or type == 'zvcz' or type == 'zdiff' or type == 'zvcd':

        freference = tempfile.NamedTemporaryFile(delete=False)
        freference.write(ref_coded)
        freference.close()
        #print_status('Ref File: ' + freference.name)        
        
        # to be optimised with bufs later
        # get length of a regarding reference
        command = '/home/dml/src/hg/dml-cliopatria/cpack/dml/scripts/compression/%s encode %s | /home/dml/src/hg/dml-cliopatria/cpack/dml/scripts/compression/length' % (type, freference.name)
        # print_status(command) 
        p1 = subprocess.Popen(command, stdin=subprocess.PIPE, stdout=subprocess.PIPE,shell=True)
        output,err = p1.communicate(input=a_coded)
        a = int(output)
        
        # get length of b regarding reference
        p1 = subprocess.Popen(command, stdin=subprocess.PIPE, stdout=subprocess.PIPE,shell=True)
        output,err = p1.communicate(input=b_coded)
        b = int(output)
        
        # get length of a,b regarding reference
        p1 = subprocess.Popen(command, stdin=subprocess.PIPE, stdout=subprocess.PIPE,shell=True)
        output,err = p1.communicate(input=a_coded + b_coded)
        ab = int(output)
   

    #print_status('Compressed Output' + compressed)
    #print_status('Compressed Size' + str(len(compressed)))
    os.remove(freference.name)
    return (a,b,ab)
    
# histogram of the last entry in a list
# returns the most frequently used key
def chord_histogram(chordstr = []):
    global chord_keys
    # build histogram 

    histo = dict.fromkeys(chord_keys,0)
    for chord in chordstr:
        histo[chord] = histo.get(chord,0) + 1 
    #print_status(str(histo.keys()))
    
    counts = np.array(histo.values(),float)
    if max(counts) > 0:
        counts = counts / max(counts)
    return (counts)