annotate experiment-reverb/code/training3.py @ 2:c87a9505f294 tip

Added LICENSE for code, removed .wav files

# -*- coding: utf-8 -*-
"""
Created on Mon Mar 21 15:29:06 2016

@author: Emmanouil Theofanis Chourdakis
"""


from glob import glob
from sys import argv
from essentia.standard import YamlInput
#from sklearn.covariance import EllipticEnvelope
from sklearn.naive_bayes import GaussianNB
from sklearn.linear_model import LogisticRegression
from sklearn import svm
from joblib import Parallel, delayed
from collections import Counter
from models import ReverbModel
from sklearn.cluster import KMeans
from time import time
import tempfile
import os
import shutil
from scipy.signal import gaussian


from sklearn.metrics import accuracy_score, f1_score
from sklearn.cross_validation import cross_val_score, cross_val_predict
from sklearn.covariance import EllipticEnvelope

from sklearn.covariance import EmpiricalCovariance, MinCovDet

from sklearn.svm import LinearSVC
from hmmlearn import hmm

from mapping import *
from numpy import *
from pca import *
from numpy.random import randint


from warnings import filterwarnings
filterwarnings("ignore")
from sklearn.base import BaseEstimator
mse = lambda A,B: ((array(A)-array(B)) ** 2).mean()

def compress_parameters_vector(parameters_vector, tol=0.000001):
    new_parameters_vector = zeros(parameters_vector.shape)
    
    for i in range(0, parameters_vector.shape[0]):
    
        p = parameters_vector[i,:]
        p_set = matrix(unique(p)).T
        
        scores = []
        
        score = -1000
        
        for n in range(1, len(p_set)):
            old_score = score    
            
            
            k_means = KMeans(n_clusters=n)
            k_means.fit(p_set)
            score = k_means.score(p_set)
            
            
            scores.append(score)
            
            if abs(score-old_score) < tol:
                n_opt = n
                break
        
        predicted_p = k_means.predict(matrix(p).T)
        new_parameters_vector[i,:] = array(ravel(k_means.cluster_centers_[predicted_p]))
    
    return new_parameters_vector
def smooth_matrix_1D(X,W=100):
    window = gaussian(W+1,4)
    window = window/sum(window)
    intermx = zeros((X.shape[0],X.shape[1]+2*W))
    intermx[:, W:-W] = X

    for m in range(0, X.shape[0]):
       # print intermx.shape
        intermx[m,:] = convolve(ravel(intermx[m,:]), window,'same')

    return intermx[:,W:-W]  
  #  return intermx

#def smooth_matrix_1D(X):
#    window = gaussian(11,4)
#    window = window/sum(window)
#    intermx = zeros((X.shape[0],X.shape[1]+20))
#    intermx[:, 10:-10] = X
#
#    for m in range(0, X.shape[0]):
#       # print intermx.shape
#        intermx[m,:] = convolve(ravel(intermx[m,:]), window,'same')
#
#    return intermx[:,10:-10] 
class HMMsvm:
    def __init__(self, svm_):
        self.svm = svm_

    def fit(self, features_list, states, flr=1e-13):

        # TODO: Concatenate features, train
        #self.svm.fit(X, states, flr)
        #self.prior = self.svm.prior
        #self.transmat = self.svm.transmat


        features = None

        for f in features_list:
            if features is None:
                features = f
            else:
                features = append(features, f, 0)

        fullstates = None
        for s in states:
            if fullstates is None:
                fullstates = s
            else:
                fullstates = append(fullstates, s)






        self.n_classes = self.svm.n_classes
        n_classes =  self.n_classes

        h = histogram(fullstates, n_classes)[0].astype(float)
        self.prior = h/sum(h)

        # Add a very small probability for random jump

        self.prior += flr
        self.prior = self.prior/sum(self.prior)

        transmat = zeros((n_classes, n_classes))
        transitions = zeros((n_classes, ))
        for s in states:
            for i in range(1, len(s)):
                prev = s[i-1]
                cur = s[i]
                transmat[prev,cur] += 1
                transitions[prev] += 1

        transitions[transitions == 0] = 1

        for i in range(0, transmat.shape[0]):
            transmat[i,:] = transmat[i,:]/transitions[i]

        self.transmat = transmat

        # Add a very small probability for random jump to avoid zero values

        self.transmat += flr

        for i in range(0, self.transmat.shape[0]):
            self.transmat[i,:] = self.transmat[i,:]/sum(self.transmat[i,:])


        A = zeros((n_classes, n_classes))

        Aold = self.transmat

        while mse(A,Aold)>10*size(A)*flr:
            Aold = A.copy()
            A = zeros((n_classes, n_classes))
            transitions = zeros((n_classes, ))

            for i in range(0,len(features_list)):
                f = features_list[i]
                s,p = self.logviterbi(f)
                for j in range(1, len(s)):
                    prev = s[j-1]
                    cur = s[j]
                    A[prev,cur] += 1
                    transitions[prev] += 1

            transitions[transitions == 0] = 1


            for i in range(0, A.shape[0]):
                A[i,:] = A[i,:]/transitions[i]

            A += flr



        self.transmat = A

        for i in range(0, A.shape[0]):
            if sum(A[i,:]) > 0:
                A[i,:] = A[i,:]/sum(A[i,:])

    def estimate_emission_probability(self, x, q):
        post = self.svm.estimate_posterior_probability_vector(x)
        prior = self.prior
        p = post/prior
        p = p/sum(p)

        return p[q]




    def predict(self, X):
        return self.logviterbi(X)[0]


    def logviterbi(self, X):
        # Number of states

        N = self.n_classes

        # Number of observations



        T = X.shape[0]



        transmat = self.transmat

        #1. Initalization

        a1 = self.prior

        delta = zeros((T,N))
        psi = zeros((T,N))

        for i in range(0, N):
            delta[0, i] = log(self.transmat[0,i]) + log(self.estimate_emission_probability(X[0,:],i))


        #2. Recursion

        for t in range(1, T):
            x = X[t, :]
            for j in range(0, N):
                delta[t,j] = max(delta[t-1,:] + log(transmat[:,j])) + log(self.estimate_emission_probability(x,j))
                psi[t,j] = argmax(delta[t-1,:] + log(transmat[:,j]))


        # 3. Termination

        p_star = max(delta[T-1,:] + log(transmat[:,N-1]))
        q_star = argmax(delta[T-1,:] + log(transmat[:, N-1]))


        # 4. Backtracking

        q = zeros((T,))
        p = zeros((T,))

        q[-1] = q_star
        p[-1] = p_star
        for t in range(1, T):
            q[-t-1] = psi[-t,   q[-t]]
            p[-t-1] = delta[-t, q[-t]]


        return q,p

    def viterbi(self, X):
        # Number of states

        N = self.n_classes

        # Number of observations

        T = X.shape[0]

        transmat = self.transmat

        #1. Initialization

        a1 = self.prior

        delta = zeros((N, ))
        psi = zeros((N, ))

        for i in range(0, N):
            delta[i] = a1[i]*self.estimate_emission_probability(X[0,:],i)





        #2. Recursion

        for t in range(1, T):
            delta_old = delta.copy()
            x = X[t,:]
            for j in range(0, N):
                delta[j] = max(delta_old*transmat[:,j])*self.estimate_emission_probability(x, j)
                psi[j] = argmax(delta_old*transmat[:,j])

        #3. Termination

        p_star = max(delta*transmat[:,N-1])
        q_star = argmax(delta*transmat[:,N-1])



        #4. Backtracking

        q = zeros((T,))
        q[-1] = q_star

        for t in range(1, T):
            q[-t-1] = psi[q[-t]]

        return q


    def _log_likelihood_matrix(self, X):
        N = self.n_classes
        ll = zeros((X.shape[0], N))

        for i in range(0, X.shape[0]):
            ll[i,:] = self._forward(i, X)

        return ll

    def compute_ksus(self, X):
        N = self.n_classes
        T = X.shape[0]
        print "[II] Computing gammas... "

        gammas = self._forward_backward(X)

        # Initialize

        aij = self.transmat






    def _not_quite_ksu(self, t, i, j, X):
        alpha = exp(self.forward(X[0:t+1,:]))[i]
        beta = exp(self.backward(X[t+1:,:]))[j]

        aij = self.transmat[i,j]
        bj = self.estimate_emission_probability(X[t+1,:], j)

        return alpha*beta*aij*bj

    def _ksu(self, t, i, j, X):
        alphaT = exp(self.forward(X[0:t+1,:]))[0]

        return self._not_quite_ksu(t,i,j,X)


    def _forward_backward(self, X):
        T = X.shape[0]
        N = self.n_classes
        fv = zeros((T, N))
        sv = zeros((T, N))

        b = None
        for t in range(1, T+1):

            fv[t-1,:] = self._forward_message(fv[t-2,:], X[0:t, ])

        for t in range(1, T+1):
            b = self._backward_message(b, X[-t: , :])
            sv[-t,:] = lognorm(fv[t-1,:]*b)

        return sv


    def _backward(self, t, X):
        N = self.n_classes
        a = ones((N,))/N

        beta = zeros((N, ))
        transmat = self.transmat
        T = X.shape[0]
        x = X[t, :]
        if t == T-1:
            beta = log(a)

            return beta
        else:
            tosum = zeros((N, ))
            bt = self._backward(t+1, X)
            btnew = zeros((N, ))
            for j in range(0, N):
                for i in range(0, N):
                    tosum[i] = log(self.estimate_emission_probability(x, j)) + bt[i] + log(transmat[i,j])

                btnew[j] = logsumexp(tosum)
            btnew = lognorm(btnew)
            return btnew


    def score(self, X):
        return self.forward(X)

    def forward(self, X):
        T = X.shape[0]
        f = None
        for t in range(1, T+1):
            f = self._forward_message(f, X[0:t, :])

        return f

    def backward(self, X):
        T = X.shape[0]
        b = None
        for t in range(1,T+1):
            b = self._backward_message(b, X[-t:, :])

        return b

    def _backward_message(self, b, X):
        N = self.n_classes


        beta = zeros((N, ))
        transmat = self.transmat

        x = X[0, :]

        if X.shape[0] == 1:
            beta = log(ones((N,))/N)
            return beta
        else:
            tosum = zeros((N, ))
            btnew = zeros((N, ))
            for j in range(0, N):
                for i in range(0, N):
                    tosum[i] = log(self.estimate_emission_probability(x,j)) + b[i] + log(transmat[i,j])

                btnew[j] = logsumexp(tosum)
            return btnew

    def _forward_message(self, f, X):
        N = self.n_classes
        a = self.prior

        alpha = zeros((N, ))
        transmat = self.transmat

        x = X[-1, :]

        if X.shape[0] == 1:
            for i in range(0, N):
                alpha[i] = log(a[i]) + log(self.estimate_emission_probability(x, i))
            return alpha

        else:
            tosum = zeros((N,))
            ftnew = zeros((N,))
            for j in range(0, N):
                for i in range(0, N):
                    tosum[i] = f[i] + log(transmat[i,j])
                Sum = logsumexp(tosum)
                bj = self.estimate_emission_probability(x, j)
                ftnew[j] = Sum+log(bj)

            return ftnew

    def _forward(self, t, X):
        N = self.n_classes
        a = self.prior
        alpha = zeros((N, ))
        transmat = self.transmat
        x = X[t, :]

        if t == 0:
            for i in range(0, N):
                alpha[i] = log(a[i]) + log(self.estimate_emission_probability(x, i))
            alpha = lognorm(alpha)

            return alpha
        else:
            tosum = zeros((N, ))
            ft = self._forward(t-1, X)
            ftnew = zeros((N, ))
            for j in range(0, N):
                for i in range(0, N):
                    tosum[i] = ft[i] + log(transmat[i,j])

                Sum = logsumexp(tosum)
                bj = self.estimate_emission_probability(x, j)
                ftnew[j] = Sum+log(bj)
            ftnew = lognorm(ftnew)

            return ftnew
            
           
class NBClassifier:
    def __init__(self):
        print "[II] Gaussian Naive Bayes Classifier"
        self.name = "Naive Bayes"
        self.smallname = "gnbc"
        self.gnb = GaussianNB()

    def getName(self):
        return self.name

    def score(self, features, parameters):
        predicted_states = self.predict(features)
        return accuracy_score(parameters, predicted_states)

    def fit(self, X, states):
        self.gnb.fit(X, states)

    def predict(self, X):
        return self.gnb.predict(X)
        
class SVMClassifier(LinearSVC):
    def __init__(self,**arguments):
        LinearSVC.__init__(self,dual=False)
        self.smallname = "svmc"
        

           
class SinkHoleClassifier(BaseEstimator):
    def __init__(self,name="SinkholeClassifier", N=2, n_components=1):
        self.chain_size = N
        self.n_components = n_components
        self.classifierNB = NBClassifier()
       # self.classifierSVM = MyAVAClassifier()
        self.classifierSVM = LinearSVC(dual=False)
        self.classifierHMM = HmmClassifier(N=N, n_components=n_components)
        self.classifierHMMSVM = HMMsvmClassifier(N=N)
        self.name = name
        self.smallname = "sinkholec"

    def score(self, features, parameters):
        predicted_states = self.predict(features)
        return accuracy_score(parameters, predicted_states)

    def getName(self):
        return self.name

    def get_params(self, deep=True):
        return {"N":self.chain_size, "n_components":self.n_components}
        
    def set_params(self, **parameters):
        for parameter, value in parameters.items():
            self.setattr(parameter, value)
        return self
        
    def fit(self, X, y):
        self.n_classes = max(unique(y))+1

        self.classifierNB.fit(X,y)
        self.classifierSVM.fit(X,y)
        self.classifierHMM.fit(X,y)
        self.classifierHMMSVM.fit(X,y)


    def predict(self, X):
        predictedNB = self.classifierNB.predict(X)
        predictedSVM = self.classifierSVM.predict(X)
        predictedHMM = self.classifierHMM.predict(X)
        predictedHMMSVM = self.classifierHMMSVM.predict(X)
        



        predicted = zeros((X.shape[0], ))

        for i in range(0, len(predicted)):
            candidates = [predictedNB[i], predictedSVM[i], predictedHMM[i], predictedHMMSVM[i]]
            
            c = Counter(candidates)

            most_common = c.most_common()

            # If there is an equal voting, select something at random
            if len(unique([k[1] for k in most_common])) == 1:
                predicted[i] = most_common[randint(len(most_common))][0]
            else:
                predicted[i] = most_common[0][0]

        return predicted            
        
class MyAVAClassifier:

    def __init__(self):
        self.classifiers = {}
        self.name  = "Linear SVM Classifier"
        self.smallname = "svc-ava"


    def getName(self):
        return self.name
    def fit(self, X, y, flr = 0, C=0.7):

        n_classes = max(unique(y)) + 1

        if len(unique(y)) == 1:
            self.only_one_class = True
            self.n_classes = 1
            self.one_class_label = y[0]
            return
        elif len(unique(y)) == 2:

            self.n_classes = n_classes
            self.svm = svm.SVC(degree=2,probability = True, kernel='rbf', gamma=2, C = C)
            self.svm.fit(X,y)
            classes_ = unique(y)
            self.classifiers[(classes_[0], classes_[1])] = self.svm
            self.only_two_classes = True
            self.only_one_class = False

            return
        else:
            self.only_two_classes = False
            self.only_one_class = False


        classes = arange(0, n_classes)
        self.n_classes = n_classes

        h = histogram(y, n_classes)[0].astype(float)
        self.prior = h/sum(h)

        transmat = zeros((n_classes, n_classes))

        for i in range(1, len(y)):
            prev = y[i-1]
            cur = y[i]
            transmat[prev,cur] += 1

        transmat = transmat/sum(transmat)

        self.transmat = transmat

        # Add a very small probability for random jump to avoid zero values

        self.transmat += flr
        self.transmat = self.transmat/sum(self.transmat)

        for i in range(0, n_classes):
            for j in range(0, n_classes):
                if i != j and (i,j) not in self.classifiers and (j,i) not in self.classifiers:

                    idx_ = bitwise_or(y == classes[i], y == classes[j])

                    X_ = X[idx_, :]

                    y_ = y[idx_]

                    if len(unique(y_)) > 1:
                        svm_ = svm.SVC(probability = True, kernel='poly', gamma=2, C = C)

                        svm_.fit(X_, y_)
                        self.classifiers[(i,j)] = svm_


    def estimate_pairwise_class_probability(self, i, j, x):


        if (i,j) not in self.classifiers and (j,i) in self.classifiers:
            return  self.classifiers[(j,i)].predict_proba(x)[0,1]
        elif (i,j) not in self.classifiers and (j,i) not in self.classifiers:
            return 0.0
        else:
            return self.classifiers[(i,j)].predict_proba(x)[0,0]

    def estimate_posterior_probability(self, i, x):
        mus = zeros((self.n_classes,))
        for j in range(0, self.n_classes):
            if i != j:
                pcp = self.estimate_pairwise_class_probability(i,j,x)
                pcp += 1e-18
                mus[j] = 1/pcp
        S = sum(mus) - (self.n_classes - 2)
        return 1/S

    def estimate_posterior_probability_vector(self, x):
        posterior = zeros((self.n_classes,))
        for i in range(0, len(posterior)):
            posterior[i] = self.estimate_posterior_probability(i, x)

        return posterior


    def predict(self, X):
        predicted = zeros((X.shape[0],))

        if self.only_one_class == True:
            return ones((X.shape[0],))*self.one_class_label
        elif self.only_two_classes == True:
            return self.svm.predict(X)


        for i in range(0, X.shape[0]):
            x = X[i,:]
            P = zeros((self.n_classes,))


            for c in range(0, len(P)):
                P[c] =  self.estimate_posterior_probability(c, x)

            pred = argmax(P)
            predicted[i] = pred

        return predicted
        
class HMMsvmClassifier(BaseEstimator):
    def __init__(self, N=2):
        self.classifiers = {}
        self.name = "HMM-SVM Classifier"
        self.smallname ="svmhmm"
        self.obs = MyAVAClassifier()
        self.chain_size = N

    def get_params(self, deep=True):
        return {"N":self.chain_size}
        
    def set_params(self, **parameters):
        for parameter, value in parameters.items():
            self.setattr(parameter, value)
        return self
        
    def getName(self):
        return self.name

    def score(self, features, parameters):
        predicted_states = self.predict(features)
        return accuracy_score(parameters, predicted_states)
        
    def fit(self, X, y):
        self.n_classes = max(unique(y))+1

        self.obs.fit(X,y)
        self.hmm = HMMsvm(self.obs)
        self.hmm.fit([X],[y])

    def predict(self, X):
        return self.hmm.predict(X)

    def confidence(self, x, q):
        return self.hmm.estimate_emission_probability(x, q)

        
class HmmClassifier(BaseEstimator):
    def __init__(self, N=2,n_components = 1):
        self.name = "HMM (%d time steps, %d components)" % (N, n_components)
        self.smallname = "ghmm"
        self.n_components = n_components
        self.chain_size = N
   

    def get_params(self, deep=True):
        return {"N":self.chain_size, "n_components":self.n_components}
        
    def set_params(self, **parameters):
        for parameter, value in parameters.items():
            self.setattr(parameter, value)
        return self

    def getName(self):
        return self.name

    def score(self, features, parameters):
        predicted_states = self.predict(features)
        return accuracy_score(parameters, predicted_states)
        
    def fit(self, features, parameters):

        n_classes = max(unique(parameters)) + 1

        if n_classes == 1:
            self.only_one_class = True
            return
        else:
            self.only_one_class = False

        hmms = [None]*n_classes
 
        chain_size = self.chain_size
        obs = [None]*n_classes

        for i in range(chain_size, len(parameters)):
            class_ = parameters[i]
            seq = features[i-chain_size:i,:]


            if obs[class_] is None:
                obs[class_] = [seq]
            else:
                obs[class_].append(seq)



        for i in range(0, len(obs)):
       
            if obs[i] is not None and len(obs[i]) != 0:
                hmm_ = hmm.GaussianHMM(n_components=self.n_components, covariance_type='diag')
                obs_ = concatenate(obs[i])
                hmm_.fit(obs_, [self.chain_size]*len(obs[i]))

                hmms[i] = hmm_

        self.hmms = hmms

        return obs

    def predict(self, features, mfilt=20):

        if self.only_one_class == True:
            return zeros((features.shape[0], ))

        chain_size = self.chain_size
        hmms = self.hmms
        predicted_classes = zeros((features.shape[0],))


        for i in range(chain_size, features.shape[0]):
            scores = zeros((len(hmms),))

            seq = features[i-chain_size:i, :]

            for j in range(0, len(hmms)):
                if hmms[j] is not None:
                    scores[j] = hmms[j].score(seq)
                else:
                    scores[j] = -infty

            predicted_classes[i] = argmax(scores)


        return predicted_classes




def vector_to_states(X):
    """
    Input: a vector MxN with N samples and M variables
    Output: a codeword dictionary `parameters_alphabet',
    state_seq, inverse `parameters_alphabet_inv' """


    parameters_alphabet = {}
    n = 0

    for i in range(0, X.shape[1]):
        vec = tuple(ravel(X[:,i]))
        if vec not in parameters_alphabet:
            parameters_alphabet[vec] = n
            n += 1

    parameters_alphabet_inv = dict([(parameters_alphabet[m],m) for m in parameters_alphabet])

    state_seq = array([parameters_alphabet[tuple(ravel(X[:,m]))] for m in range(0, parameters_vector.shape[1])]    )

    return state_seq, parameters_alphabet, parameters_alphabet_inv


def states_to_vector(predicted, parameters_alphabet_inv):
    estimated = matrix(zeros((len(parameters_alphabet_inv[0]), len(predicted))))
    for i in range(0, len(predicted)):
     #   print "[II] Predicted: ", predicted[i]
        estimated[:, i] = matrix(parameters_alphabet_inv[predicted[i]]).T

    return estimated

if __name__=="__main__":
    if len(argv) != 2:
        print "[EE] Wrong number of arguments"
        print "[II] Correct syntax is:"
        print "[II] \t%s <trainingdir>"
        print "[II] where <trainingdir> is the directory containing the features files in the format *_features.yaml"
        print "[II] and the parameters in the format *_parameters.yaml"
        exit(-1)
        
    totalstart = time()
    
    traindir = argv[1]
    
    songs_in_dir = glob("%s/*_features.yaml" % traindir)

    
    # Shuffle to avoid using the same songs in the same sets in cross-validation
    
    #random.shuffle(songs_in_dir)


    features_vector = None
    parameters_vector = None
    index_vector = None

  
    idx = 0    
    
    for f_ in songs_in_dir:
        infile = f_
        paramfile = "%s_parameters.yaml" % f_.split('_features.yaml')[0]

        print "[II] Using features file: %s" % infile
        print "[II] and parameters file: %s" % paramfile


        features_pool = YamlInput(filename = infile)()
        parameters_pool = YamlInput(filename = paramfile)()

        d1 = parameters_pool['parameters.d1']
        da = parameters_pool['parameters.da']
        g1 = parameters_pool['parameters.g1']
        gc = parameters_pool['parameters.gc']
        G = parameters_pool['parameters.G']    
        
        feature_captions = features_pool.descriptorNames()
        parameter_captions = ['d1','da','g1','gc','G']
        
        
        for c in features_pool.descriptorNames():
            if c.split('.')[0] == 'metadata' or c.split('.')[0] == 'parameter':
                feature_captions.remove(c)
    
    
        print "[II] %d Features Available: " % len(feature_captions)
        print str(feature_captions).replace("', ","\n").replace('[','').replace("'","[II]\t ")[:-7]
        
        nfeatures_in = len(feature_captions)
        nparameters_in = 5
        features_vector_one = zeros((nfeatures_in, len(features_pool[feature_captions[0]])))
        index_vector_one = ones((len(features_pool[feature_captions[0]]),))
        parameters_vector_one = ones((nparameters_in, len(features_pool[feature_captions[0]])))


        for i in range(0, nfeatures_in):
            features_vector_one[i, :] = features_pool[feature_captions[i]].T


        parameters_vector_one[0, :] = d1*parameters_vector_one[0,:]
        parameters_vector_one[1, :] = g1*parameters_vector_one[1,:]
        parameters_vector_one[2, :] = da*parameters_vector_one[2,:]
        parameters_vector_one[3, :] = gc*parameters_vector_one[3,:]
        parameters_vector_one[4, :] = G*parameters_vector_one[4,:]    
        # Stores example index number

        index_vector_one *= idx

        print "[II] %d parameters used:" % len(parameter_captions)
        print str(parameter_captions).replace("', ","\n").replace('[','').replace("'","[II]\t ")[:-7].replace('parameter.','')

        if features_vector is None:
            features_vector = features_vector_one
        else:
            features_vector = append(features_vector, features_vector_one, 1)

        if parameters_vector is None:
            parameters_vector = parameters_vector_one
        else:
            parameters_vector = append(parameters_vector, parameters_vector_one, 1)

        if index_vector is None:
            index_vector = index_vector_one
        else:
            index_vector = append(index_vector, index_vector_one)


        idx += 1        
    
    print "[II] Marking silent parts"
    #features_full_nontransformed_train = features_vector.copy()

    silent_parts = zeros((1, features_vector.shape[1]))

    rms = features_vector[feature_captions.index('rms'), :]

    # Implementing Hysteresis Gate -- High threshold is halfway between
    # the mean and the max and Low is halfway between the mean dn the min

    rms_threshold_mean = mean(rms)

    rms_threshold_max = max(rms)
    rms_threshold_min = min(rms)

    rms_threshold_high = 0.1 * rms_threshold_mean
    rms_threshold_low = 0.01 * rms_threshold_mean

    for n in range(1,  len(rms)):
        prev = rms[n-1]
        curr = rms[n]

        if prev >= rms_threshold_high:
            if curr < rms_threshold_low:
                silent_parts[0,n] = 1
            else:
                silent_parts[0,n] = 0
        elif prev <= rms_threshold_low:
            if curr > rms_threshold_high:
                silent_parts[0,n] = 0
            else:
                silent_parts[0,n] = 1
        else:
            silent_parts[0,n] = silent_parts[0,n-1]


    if silent_parts[0,1] == 1:
        silent_parts[0, 0] = 1


    ## START COMMENTING HERE FOR KEEPING ACTIVE PARTS ONLY

    active_index = invert(silent_parts.flatten().astype(bool))
    # Keep only active parts

    # Uncomment this
    features_vector = features_vector[:, active_index]

    # Keep this for debugging purposes
    parameters_vector = parameters_vector[:, active_index]    
    
    
  
    
    ## UP TO HERE
    parameters_vector[0,:] = (parameters_vector[0,:] - d1_min)/d1_max #d1
    parameters_vector[1,:] = (parameters_vector[1,:] - g1_min)/g1_max #g1
    parameters_vector[2,:] = (parameters_vector[2,:] - da_min)/da_max #g1
    parameters_vector[3,:] = (parameters_vector[3,:] - G_min)/G_max #G
    parameters_vector[4,:] = (parameters_vector[4,:] - gc_min)/gc_max #gc    

    # Store moments    
    moments_vector = zeros((features_vector.shape[0], 2))
    features_vector_original = features_vector.copy()    
    
    print "[II] Storing moments vector"
    for i in range(0, features_vector.shape[0]):
        mean_ = mean(features_vector[i,:])
        std_ = std(features_vector[i,:], ddof=1)
        moments_vector[i,0] = mean_
        moments_vector[i,1] = std_

        features_vector[i,:] = (features_vector[i,:] - mean_)/std_

    features_vector_old_train = features_vector.copy()    
    
    

    

    
    
    print "[II] Extracting PCA configuration "
    kernel, q, featurelist = extract_pca_configuration_from_data(features_vector)   

    q = kernel.shape[0]

    
    print "[II] Optimal number of PCs to keep: %d" % q
    
    
    
    feature_captions_array = array(feature_captions)
    features_to_keep = list(feature_captions_array[featurelist])

    
    print "[II] Decided to keep %d features:" % len(features_to_keep)
    print  str(features_to_keep).replace("', ","\n").replace('[','').replace("'","[II]\t ")[:-7]

    featurelist_bak = featurelist
    features_kept_data = features_vector[featurelist,:]
    features_kept_data_bak = features_kept_data.copy()
    features_vector = (kernel.T*features_kept_data)[0:q,:]
    features_vector_new_train = features_vector.copy()    
    
    
    features_vector_unsmoothed = features_vector
    features_vector = smooth_matrix_1D(features_vector)    
    

    # Construct parameters alphabet, each symbol is going to be a different column vector
    # in parameter code matrix
    
    parameters_state, parameter_state_alphabet, parameter_state_alphabet_inv = vector_to_states(parameters_vector)
    
    
    
  
    CV = 10
    print "[II] %d-fold cross validation" % CV
    
    
    print "[II] GaussianNB..."    
    gnbc = GaussianNB()
    predictsGNBC = cross_val_predict(gnbc,features_vector.T,parameters_state,cv=CV)    
    scoresGNBC = f1_score(parameters_state, predictsGNBC, average='weighted')
    parametersGNBC = states_to_vector(predictsGNBC, parameter_state_alphabet_inv)
    msesGNBC = mse(parametersGNBC, parameters_vector)
    print "[II] f1-weighted: %s" % scoresGNBC    
    print "[II] mse: %s" % msesGNBC
    
    print "[II] o-v-r SVM..."
    ovrc = LinearSVC(dual=False)
    predictsOVRC = cross_val_predict(ovrc,features_vector.T,parameters_state,cv=CV)   
    scoresOVRC = f1_score(parameters_state, predictsOVRC, average='weighted')
    parametersOVRC = states_to_vector(predictsOVRC, parameter_state_alphabet_inv)
    msesOVRC = mse(parametersOVRC, parameters_vector)    
    print "[II] %s" % scoresOVRC
    print "[II] mse: %s" % msesOVRC
    

    print "[II] GaussianHMM..."  
    start = time()
   

    
    
    tempfolder = tempfile.mkdtemp()
 
    _scores_cs = memmap(os.path.join(tempfolder,'scores'), dtype=float,shape=98, mode='w+')  
    _mses_cs = memmap(os.path.join(tempfolder,'mses'), dtype=float,shape=98, mode='w+')    
    
    _chain_sizes = memmap(os.path.join(tempfolder,'cs'), dtype=int,shape=98, mode='w+')    

    
    print "[II] Calibrating GaussianHMM please wait..."
    def iter_(scores, cs, n):
        _n = n
        
        try:
            hmmc = HmmClassifier(N=_n, n_components=1)
            predictsHMMC = cross_val_predict(hmmc,features_vector.T,parameters_state,cv=CV) 
       # scoresHMMC = cross_val_score(hmmc,features_vector.T,parameters_state,cv=CV, scoring='f1_weighted')
            scores[n-2] = f1_score(parameters_state, predictsHMMC, average='weighted')
            parametersHMMC = states_to_vector(predictsHMMC, parameter_state_alphabet_inv) 
            _mses_cs[n-2] = mse(parametersHMMC, parameters_vector)
            
        except: 
            scores[n-2] = 0
        cs[n-2] = _n
        
        


    print "[II] Calibrating GaussianHMM please wait..."

        
    def iter_(scores, nc, n):
        try:
            hmmc = HmmClassifier(N=opt_N, n_components=n)
            predictsHMMC = cross_val_predict(hmmc,features_vector.T,parameters_state,cv=CV) 
        
     #   scoresHMMC = cross_val_score(hmmc,features_vector.T,parameters_state,cv=CV,n_jobs=2, scoring='f1_weighted')
            scores[n-1] = f1_score(parameters_state, predictsHMMC, average='weighted')
            parametersHMMC = states_to_vector(predictsHMMC, parameter_state_alphabet_inv) 
            _mses_nc[n-1] = mse(parametersHMMC, parameters_vector)            
        except:
            scores[n-1] = 0
        nc[n-1] = n

        
    print "[II] Configuring using cross-validation disabled for performance reasons, please ignore results or uncomment the relevant code."


    

    finish = time()

        
    
    outp = ""
    
    outp+= "[II] Results (F1-Scores):\n"

    outp+= "[II] Estimator\tMean\n"
    outp+= "[II] %s\t%.2f\t%.5f\n" % ("GaussianNB", scoresGNBC, msesGNBC)
    outp+= "[II] %s\t%.2f\t%.5f\n" % ("O-V-R SVM", scoresOVRC, msesOVRC)


    
    totalfinish = time()
    totalduration = totalfinish - totalstart
    
    outp+= "[II] Total time: %dm:%ds\n" % (int(totalduration) / 60, int(totalduration) % 60)
    print outp
    shutil.rmtree(tempfolder)

    gnbc = NBClassifier()
    gnbc.fit(features_vector.T, parameters_state)
    
    outp+= "[GNBC] Testing against full data:\n"    
    predictsGNBC = gnbc.predict(features_vector.T)
    scoresGNBC = f1_score(parameters_state, predictsGNBC, average='weighted')
    parametersGNBC = states_to_vector(predictsGNBC, parameter_state_alphabet_inv)
    msesGNBC = mse(parametersGNBC, parameters_vector)
    outp+= "[II] f1-weighted: %s\n" % scoresGNBC    
    outp+= "[II] mse: %s\n" % msesGNBC
    

        
    

    ovrc = SVMClassifier()
    ovrc.fit(features_vector.T, parameters_state)
    
    outp+= "[OVRC ] Testing against full data:\n"    
    predictsOVRC  = ovrc.predict(features_vector.T)
    scoresOVRC  = f1_score(parameters_state, predictsOVRC , average='weighted')
    parametersOVRC  = states_to_vector(predictsOVRC , parameter_state_alphabet_inv)
    msesOVRC  = mse(parametersOVRC , parameters_vector)
    outp+= "[II] f1-weighted: %s\n" % scoresOVRC    
    outp+= "[II] mse: %s\n" % msesOVRC        
    
    

    model_gnb = ReverbModel("Naive Bayes Classifier Model", kernel, q, features_to_keep, gnbc, parameter_dictionary=parameter_state_alphabet_inv, moments_vector=moments_vector  )
    model_svm = ReverbModel("OVR LinearSVM Classifier Model", kernel, q, features_to_keep, ovrc, parameter_dictionary=parameter_state_alphabet_inv, moments_vector=moments_vector  )


    model_gnb.save("model_gnb.dat")
    model_svm.save("model_svm.dat")


    print outp
    f = open('results.txt','w')
    f.write(outp)
    f.close()