view experiment-reverb/code/#training.py# @ 2:c87a9505f294 tip

Added LICENSE for code, removed .wav files

#!/usr/bin/python2
# -*- coding: utf-8 -*-
"""
Created on Thu Apr 23 11:53:17 2015

@author: mmxgn
"""

# This file does the cluster estimation and the removal of outliers

from warnings import filterwarnings

from sys import argv, exit
from essentia.standard import YamlInput, YamlOutput
from essentia import Pool
from pca import *

from numpy import *
from sklearn import cluster
from sklearn.metrics import pairwise_distances
from sklearn.cluster import KMeans, MiniBatchKMeans
from matplotlib.pyplot  import *
#from sklearn.mixture import GMM
from sklearn.naive_bayes import GaussianNB, MultinomialNB
from scipy.signal import decimate
from sklearn import cross_validation
from sklearn import multiclass
from sklearn.covariance import EllipticEnvelope

#from hmmlearn import hmm
from hmmlearn.hmm import GMM
from hmmlearn import hmm

from sklearn import svm
import copy as cp

from scipy.misc import logsumexp
import pickle

from collections import Counter
#from adpcm import adm, adm_reconstruct


from reshape_observations import reshape_observations

mse = lambda A,B: ((array(A)-array(B)) ** 2).mean()

lognorm = lambda A: A - logsumexp(A)


rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
## for Palatino and other serif fonts use:
#rc('font',**{'family':'serif','serif':['Palatino']})
rc('text', usetex=True)
#logsum = lambda X: logsumexp(log(X))


filterwarnings("ignore")
C = 0.7
NCOMPONENTS=2
CHAINSIZE = 5

close('all')

from numpy.random import randint

class HMMsvm:
    def __init__(self, svm_):
        self.svm = svm_


      #  self.svm = MyAVAClassifier()


    def fit(self, features_list, states, flr=1e-13):
 #       def fit(self, features_list, 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
#            T = features.shape[0]
        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

  #      print fullstates, shape(fullstates)

        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]
  #              print "FEATURES:", f
             #   print f
                s,p = self.logviterbi(f)
            #    print s
                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,:])



                 #   print observations


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

        return p[q]

 #   def estimate_emission_probability(self, x, q):
 #       return self.svm.estimate_emission_probability(q, x)


    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):
        #    delta_old = delta.copy()
            x = X[t, :]
            for j in range(0, N):
              #  print "j: %d, S" % j, delta_old+log(transmat[:,j])
                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]))

    #        print "t: %d, delta: " % t, delta
     #       print "t: %d, psi:" % t, psi


        # 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, ))
         #   print bt
            for j in range(0, N):
                for i in range(0, N):
                  #  print log(self.estimate_emission_probability(x, j)),  bt[i], log(transmat[i,j])
                    tosum[i] = log(self.estimate_emission_probability(x, j)) + bt[i] + log(transmat[i,j])
              #  print tosum

                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):
         #   print t
    #        print t,b
            idx = arange(-t,T)
     #       print idx
            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)
        #    btnew = lognorm(btnew)
            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))
         #   alpha = lognorm(alpha)
            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)

        #    ftnew = lognorm(ftnew)
            return ftnew

    def _forward(self, t, X):
        N = self.n_classes
        a = self.prior
    #    print a
        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))
           # print "--"
          #  print alpha
            alpha = lognorm(alpha)
          #  print alpha
          #  print "xx"
            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):
               #     print ft
                    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 HMMsvm2:
    def __init__(self, svm_, n_components=4):
        self.svm = svm_
        self.n_components = n_components


      #  self.svm = MyAVAClassifier()


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

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


        n_classes = self.n_components
        self.n_classes = n_classes

        A = zeros((n_classes, n_classes))

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

        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]
  #              print "FEATURES:", f
             #   print f
                s,p = self.logviterbi(f)
            #    print s

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


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


                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,:])



                 #   print observations


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

        return p[q]

 #   def estimate_emission_probability(self, x, q):
 #       return self.svm.estimate_emission_probability(q, x)


    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):
        #    delta_old = delta.copy()
            x = X[t, :]
            for j in range(0, N):
              #  print "j: %d, S" % j, delta_old+log(transmat[:,j])
                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]))

    #        print "t: %d, delta: " % t, delta
     #       print "t: %d, psi:" % t, psi


        # 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, ))
         #   print bt
            for j in range(0, N):
                for i in range(0, N):
                  #  print log(self.estimate_emission_probability(x, j)),  bt[i], log(transmat[i,j])
                    tosum[i] = log(self.estimate_emission_probability(x, j)) + bt[i] + log(transmat[i,j])
              #  print tosum

                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):
         #   print t
    #        print t,b
            idx = arange(-t,T)
     #       print idx
            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)
        #    btnew = lognorm(btnew)
            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))
         #   alpha = lognorm(alpha)
            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)

        #    ftnew = lognorm(ftnew)
            return ftnew

    def _forward(self, t, X):
        N = self.n_classes
        a = self.prior
    #    print a
        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))
           # print "--"
          #  print alpha
            alpha = lognorm(alpha)
          #  print alpha
          #  print "xx"
            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):
               #     print ft
                    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 fit(self, X, states):
        self.gnb.fit(X, states)

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

    def cross_validate(self, data, classes):
        print "[II] SVN Classifier Crossvalidating... "
        from copy import deepcopy
        estimator = deepcopy(self)
        estimator_fulldata = deepcopy(self)
        estimator_fulldata.fit(data, classes)

        percents = arange(0.1,0.9,0.1)
        self.percents = percents * 100

        RATIOS = []
        labels = estimator.predict(data)


        print "[II] for full training-testing: %.2f" %  (sum(array(classes==labels).astype(float))/len(labels))

        for p in percents:
            traindata,testdata,trainlabels,testlabels = cross_validation.train_test_split(data,classes,test_size=p,random_state=1)

   #             print "[II] [%.2f-%.2f] Examples used for training: %s" % (p, (1-p), str(traindata))
   #             print "[II] [%.2f-%.2f] Examples used for testing: %s" % (p, (1-p), str(testdata))

            # Training phase



            estimator_ = deepcopy(estimator)
            estimator_.fit(traindata, trainlabels)

            predicted_labels = estimator_.predict(testdata)

            m = sum(array(predicted_labels==testlabels).astype(float))/len(testlabels)

            print "[II] [%.2f-%.2f] 5-fold cross-validation for NB Classifier: %.1f" % (p,  (1-p), m)

            RATIOS.append(m)

        return RATIOS,percents




class HmmClassifier3:
    def __init__(self, N=2,n_components = 1):
        print "[II] HMM Classifier (%d states, %d components)" % (N, n_components)
        self.name = "HMM (%d time steps, %d components)" % (N, n_components)
        self.n_components = n_components

        self.n_components = n_components
        self.smallname = "hmmc"
        self.chain_size = N

    def getName(self):
        return self.name

    def fit(self, features, parameters):


#        print "Parameters:"
#        print parameters
        n_classes = max(unique(parameters)) + 1

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

        print "[II] Number of classes: ", n_classes
        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')
                # print(obs[i])
                #obs_ = reshape_observations(obs[i])
                obs_ = concatenate(obs[i])
                #lengths_ = [o.shape[0] for o in obs[i]]
                #print (lengths_)
            #    print len(obs[i])
            #    print obs[i][0][0].shape
            #    print obs[i]
      #          for s in obs[i]:
      #              if s.ndim > 2:
      #                  print s
                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)

        # Do a median filter at the end

  #      predicted_classes = median_filter(predicted_classes,mfilt)

        return predicted_classes



    def k_fold_cross_validate(self, data, classes, K=5):
        print "[II] HMM Classifier Crossvalidating... "
        print "[II] Validating data: ", data
        print "[II] With classes: ", classes
        from copy import deepcopy
        estimator = deepcopy(self)
        estimator_fulldata = deepcopy(self)
        estimator_fulldata.fit(data, classes)


        labels = estimator_fulldata.predict(data)
        print "[II] for full training-testing: %.2f" %  (sum(array(classes==labels).astype(float))/len(classes))

        # Crossvalidate this.

        example_sequences = []
        example_labels = []

        chain_size = self.chain_size

        percents = arange(0.1,0.9,0.1)
        self.percents = percents * 100

        RATIOS = []



        # Constructing examples and labels


        L = data.shape[0]/K
        M = K
 #       print "[II] Example size:", L

        n = 1
        while L*(n-1) < M*L:
            if L*n > shape(data)[0]:
                example = data[L*(n-1):,:]
                classes_ = classes[L*(n-1):]
            else:
                example = data[L*(n-1):L*n, :]
                classes_ = classes[L*(n-1):L*n]

         #   print example
           # print classes_

            example_sequences.append(example)
            example_labels.append(classes_)
            n+=1

 #       print len(example_sequences)
#        print len(example_labels)
        from sklearn.cross_validation import KFold
        kf = KFold(K, n_folds=K)

        ratio = 0

        for train, test in kf:
            print "[II] Validating examples %s against %s." % (train, test)

            topredict_data = example_sequences[test[0]]
            topredict_labels = example_labels[test[0]]

            tofit_data = None
            tofit_labels = None

            for t in train:
      #          print t
                if tofit_data is None:
                    tofit_data = example_sequences[t]
                    tofit_labels = example_labels[t]
                else:
                    tofit_data = append(tofit_data, example_sequences[t], 0)
                    tofit_labels = append(tofit_labels, example_labels[t], 0)

            estimator_ = deepcopy(estimator)
            estimator_.fit(tofit_data, tofit_labels)

            labels = estimator_.predict(topredict_data)

            m = sum(array(topredict_labels==labels).astype(float))/len(labels)


            print "[II] %d-fold cross-validation for HMM Classifier: %.1f" % (K, m)

            ratio += m/float(len(kf))

        return ratio






        # Splitting to train/test


    def cross_validate(self, data, classes):
        print "[II] HMM Classifier Crossvalidating... "
        from copy import deepcopy
        estimator = deepcopy(self)
        estimator_fulldata = deepcopy(self)
        estimator_fulldata.fit(data, classes)


        labels = estimator_fulldata.predict(data)
        print "[II] for full training-testing: %.2f" %  (sum(array(classes==labels).astype(float))/len(classes))

        # Crossvalidate this.

        example_sequences = []
        example_labels = []

        chain_size = self.chain_size

        percents = arange(0.5,0.9,0.1)
        self.percents = percents * 100

        RATIOS = []



        # Constructing examples and labels

        M = self.chain_size
        L = data.shape[0]/20

        print "[II] Example size:", L

        n = 1
        while L*n < M*L-L:
            example = data[L*(n-1):L*n, :]
         #   print example
            classes_ = classes[L*(n-1):L*n]
           # print classes_

            example_sequences.append(example)
            example_labels.append(classes_)
            n+=1



        # Splitting to train/test

        for p in percents:
            traindata,testdata,trainlabels,testlabels = cross_validation.train_test_split(example_sequences, example_labels, test_size=p, random_state=1)

            print traindata
            print testdata
            # Concatenate traindata

            tofit_data = None
            tofit_labels = None

            topredict_data = None
            topredict_labels = None


            for t in traindata:
                if tofit_data is None:
                    tofit_data = t
                else:
                    tofit_data = append(tofit_data, t, 0)

            for l in trainlabels:
                if tofit_labels is None:
                    tofit_labels = l
                else:
                    tofit_labels = append(tofit_labels, l)

            for t in testdata:
                if topredict_data is None:
                    topredict_data = t
                else:
                    topredict_data = append(topredict_data, t, 0)

            for l in testlabels:
                if topredict_labels is None:
                    topredict_labels = l
                else:
                    topredict_labels = append(topredict_labels, l)


       #     print "[II] Fit data: ", tofit_data
       ##     print "[II] Fit labels: ", tofit_labels
       #     print "[II] Predict data: ", topredict_data
       #     print "[II] Predict labels: ", topredict_labels
        #
            estimator_ = deepcopy(estimator)

            print tofit_labels
            estimator_.fit(tofit_data, tofit_labels)

            labels = estimator_.predict(topredict_data)

            m = sum(array(topredict_labels==labels).astype(float))/len(labels)


            print "[II] [%.2f-%.2f] cross-validation for HMM Classifier: %.1f" % (p,  (1-p), m)

            RATIOS.append(m)

        return RATIOS,percents


class HMMsvmClassifier2:
    def __init__(self, N=2):
        self.classifiers = {}
        self.name = "HMM-SVM Classifier 2"
        self.smallname = "hmmsvm2"
        self.obs = MyAVAClassifier()
        self.chain_size = N

    def getName(self):
        return self.name

    def fit(self, X, y):
        self.n_classes = max(unique(y))+1
        n_classes = self.n_classes
        chain_size = self.chain_size

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

            print "[II] Number of classes: ", n_classes

            hmms = [None]*n_classes
            obs = [None]*n_classes


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



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

    def getName(self):
        return self.name

    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):
   #     chain_size = self.chain_size
     #   predicted = zeros((X.shape[0]),)
    #    for i in range(chain_size, len(predicted)):
    #        predicted[i] = self.hmm.predict(X[i-chain_size:i,:])[-1]



      #  return predicted

        return self.hmm.predict(X)

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



class SinkHoleClassifier:
    def __init__(self):
        self.classifierNB = NBClassifier()
        self.classifierSVM = MyAVAClassifier()
        self.classifierHMM = HmmClassifier3()
        self.classifierHMMSVM = HMMsvmClassifier()



    def getName(self):
        return self.name


    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"
     #   self.probabilities = {}


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

        n_classes = max(unique(y)) + 1

        # Transform to a binary if there are only two classes
        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(probability = True, kernel='poly', 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
#        M = n_classes*(n_classes-1) # Number of classifiers

        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:
                   # print (i,j)
                   # print classes
                    idx_ = bitwise_or(y == classes[i], y == classes[j])
                 #   idx_ = bitwise_or(y == i, y == j)

                    X_ = X[idx_, :]

                    y_ = y[idx_]
                    #print y_

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

                  #      print y_
                        svm_.fit(X_, y_)
                        self.classifiers[(i,j)] = svm_
                #    self.probabilities[(i,j)] = svm_.predict_proba(X)


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

       # print self.classifiers.keys()

        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
    #    print mus
        S = sum(mus) - (self.n_classes - 2)
    #    print S
        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 estimate_emission_probability(self, i, x):
#            post = self.estimate_posterior_probability_vector(x)
#         #   print "post: ", post
#            prior = self.prior
#         #   print "prior: ", prior
#            p = post/prior
#            p = p/sum(p)
#
#            return p[i]

#        def estimate_emissmat(self, X):
#            emissmat = zeros((X.shape[0], self.n_classes))
#            for i in range(0, X.shape[0]):
#                for j in range(0, self.n_classes):
#                    emissmat[i,j] = self.estimate_emission_probability(j, X[i,:])
#
#            return emissmat

    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

def logsum(X):
    if len(X) == 1:
        return log(X[0])
    else:
        return logaddexp(log(X[0]),logsum(X[1:]))




def smooth_matrix_1D(X):
    window = scipy.signal.gaussian(51,4)
    window = window/sum(window)
    intermx = zeros((X.shape[0],X.shape[1]+100))
    intermx[:, 50:-50] = X

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

    return intermx[:,50:-50]

def adm_reconstruct(codeword, h, dmin=.01, dmax=.28):
    x = zeros((1, codeword.shape[1]))

    delta1 = dmin
    delta2 = dmin
    Sum = h

    x[0] = h
    for i in range(0, codeword.shape[1]):
        if codeword[0,i] == 0:
            delta1 = dmin
            delta2 = dmin

        elif codeword[0,i] == 1:
            delta2 = dmin
            Sum += delta1
            delta1 *= 2
            if delta1 > dmax:
                delta1 = dmax

        elif codeword[0,i] == -1:
            delta1 = dmin
            Sum -= delta2
            delta2 *= 2
            if delta2 > dmax:
                delta2 = dmax
        x[0,i] = Sum
    return x

def adm(x, dmin=.01, dmax=.28, tol=0.0001):

    # Adaptive delta modulation adapted by code:
    # (adeltamod.m)
    #
    #   Adaptive Delta Modulator
    #   by Gandhar Desai (gdesai)
    #   BITS Pilani Goa Campus
    #   Date: 28 Sept, 2013

    xsig = x

    Lx = len(x)

    ADMout = zeros((1, Lx))
    codevec = zeros((1, Lx))


    Sum = x[0]
    delta1 = dmin
    delta2 = dmin
    mult1 = 2
    mult2 = 2
    for i in range(0, Lx):
        #print abs(xsig[i] - Sum)
        if (abs(xsig[i] - Sum) < tol):
            bit = 0
            delta2 = dmin
            delta1 = dmin


        elif (xsig[i] >= Sum):
            bit = 1
            delta2 = dmin
            Sum += delta1
            delta1 *=  mult1
            if delta1 > dmax:
                delta1 = dmax


        else:
            bit = -1
            delta1 = dmin
            Sum -= delta2
            delta2 *= mult2
            if delta2 > dmax:
                delta2 = dmax



        ADMout[0, i] = Sum
        codevec[0, i]= bit

    return ADMout,codevec, x[0]

def median_filter(v, K):
    v2 = zeros((len(v),))
    for i in range(K, len(v)):
        seq = v[i-K:i]
        m = median(seq)
        v2[i-K:i] = m

    return v2




from models import ReverbModel



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)


    n_clusters = 3
    UpsamplingFactor = 1
    dmin = 0.001
    dmax = 0.28
    tol = 0.001

    d1min = 0.01
    d1max = 0.1

    g1min = 0.01
    g1max = 0.99

    damin = 0.006
    damax = 0.012

    Gmin = 0.01
    Gmax = 0.99

    gcmin = 0.01
    gcmax = 0.99

    # For searching the directory
    from glob import glob

    traindir = argv[1]

    songs_in_dir = glob("%s/*_features.yaml" % traindir)

    print "[II] Using files: %s" % songs_in_dir


    chain_length=1


   # asdsd

   # fullfeatures_pool = Pool()

    features_vector = None
    parameters_vector = None
    index_vector = None

    idx = 0
    # Shuffle segments
    random.shuffle(songs_in_dir)
    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


      #  paramfile = infile.split(')

        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']


        svmhmmstr = ""


        for c in features_pool.descriptorNames():
            if c.split('.')[0] == 'metadata' or c.split('.')[0] == 'parameter':
                feature_captions.remove(c)


    #    close('all')

    #    print "[II] Loaded training data from %s (%s) " % (infile, features_pool['metadata.date'][0])
        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]])))
       # paramet


        for i in range(0, nfeatures_in):
            features_vector_one[i, :] = features_pool[feature_captions[i]].T
       # for i in range(0, nparameters_in):
       #     parameters_vector[i, :] = features_pool[parameter_captions[0]].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, len(features_pool[feature_captions[i]].T)))

    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



    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]
  #  index_vector = index_vector[active_index]

    # Label examples for a chain of chain_length

#    idx = 0
#    for i in range(0, len(index_vector)):
#        index_vector[i] = idx
#        if i % chain_length == 0:
#            idx += 1
#    number_of_examples = idx



    # Scale parameters to [0,1]


    parameters_vector[0,:] = (parameters_vector[0,:] - d1min)/d1max #d1
    parameters_vector[1,:] = (parameters_vector[1,:] - g1min)/g1max #g1
    parameters_vector[2,:] = (parameters_vector[2,:] - damin)/damax #g1
    parameters_vector[3,:] = (parameters_vector[3,:] - Gmin)/Gmax #G
    parameters_vector[4,:] = (parameters_vector[4,:] - gcmin)/gcmax #gc







    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()
  #  moments_vector_train = moments_vector


    print "[II] Extracting PCA configuration "

    kernel, q, featurelist = extract_pca_configuration_from_data(features_vector)


    #q = 1


   # q = q + 1

    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_old_train = features_kept_data
    features_vector = (kernel.T*features_kept_data)[0:q,:]
    features_vector_new_train = features_vector.copy()


#    example_features = None
#    example_parameters = None
#    example_idx = None
#
#    for idx in range(0, shape(parameters_vector)[1]-chain_length):
#        example_features_ = features_vector[:, idx:idx+chain_length]
#      #  print example_features_
#        example_parameters_ = parameters_vector[:, idx:idx+chain_length]
#        example_idx_ = ones((shape(example_parameters_)[1],))
#
#        if example_features is None:
#            example_features = example_features_
#        else:
#            #print "appending"
#            example_features = append(example_features, example_features_, 1)
#
#        if example_parameters is None:
#            example_parameters = example_parameters_
#        else:
#            example_parameters = append(example_parameters, example_parameters_, 1)
#
#        if example_idx is None:
#            example_idx = example_idx_*idx
#        else:
#            example_idx = append(example_idx, example_idx_*idx, 1)
#
#
#
#
#    features_vector = example_features
#    parameters_vector = example_parameters
#    index_vector = example_idx

    print "[II] Trying ADM-coded parameters"
    print "[II] Upsampling features and parameters by a factor of %d" % UpsamplingFactor


    # Upsampled features and parameters
    features_vector_upsampled = smooth_matrix_1D(repeat(features_vector,UpsamplingFactor, axis=1))
 #   features_vector_upsampled = repeat(features_vector, UpsamplingFactor, axis=1)


  #  features_vector_upsampled = smooth_matrix_1D(repeat(features_vector_old_train,UpsamplingFactor, axis=1))

   # feature_labels_upsampled = repeat(features_clustering_labels,UpsamplingFactor, axis=0)
    parameters_vector_upsampled = repeat(parameters_vector,UpsamplingFactor, axis=1)


 #   km = KMeans(init='k-means++', n_init=10, max_iter=300, tol=0.0000001, verbose = 0, n_clusters=12).fit(parameters_vector_upsampled.T)
  #  centers_ = km.cluster_centers_
 #   labels__ = km.labels_
    # Remove the following two lines in order to restore non-kmeans version
 #   parameters_vector_kmeans_upsampled = array(centers_[labels__])
#
  #  parameters_vector_upsampled = parameters_vector_kmeans_upsampled.T
#
#
  #  parameters_vector_upsampled = smooth_matrix_1D(parameters_vector_upsampled)

    parameters_vector_upsampled_adm = matrix(zeros(shape(parameters_vector_upsampled)))
    parameters_vector_upsampled_code = matrix(zeros(shape(parameters_vector_upsampled)))
    parameters_vector_upsampled_firstval = matrix(zeros((parameters_vector_upsampled.shape[0],1)))

    # Reconstructed parameters



    parameters_vector_upsampled_reconstructed = matrix(zeros(shape(parameters_vector_upsampled)))



    # ADM HERE
    #

    def adm_matrix(X, dmin=0.001,dmax=0.28,tol=0.001):

        out = matrix(zeros(shape(X)))
        code = matrix(zeros(shape(X)))
        firstval = matrix(zeros((X.shape[0], 1)))

        for i in range(0, X.shape[0]):
            out[i,:], code[i,:], firstval[i,0] = adm(X[i,:],dmin=dmin,dmax=dmax,tol=tol)

        return out,code,firstval

#        parameters_vector_upsampled_reconstructed[i,:] = adm_reconstruct(parameters_vector_upsampled_code[i,:],parameters_vector_upsampled_firstval[i,0], dmin=dmin,dmax=dmax)

    def adm_matrix_reconstruct(code, firstval, dmin=0.001, dmax=0.28):
        X = matrix(zeros(shape(code)))
        for i in range(0, code.shape[0]):
            X[i,:] = adm_reconstruct(code[i,:], firstval[i,0], dmin=dmin, dmax=dmax)

        return X


    parameters_vector_upsampled_adm, parameters_vector_upsampled_code, parameters_vector_upsampled_firstval = adm_matrix(parameters_vector_upsampled, dmin, dmax, tol)
  #  parameters_vector_upsampled = parameters_vector_upsampled_code


    def diff_and_pad(X):
        return concatenate((
                zeros((
                    shape(X)[0],
                    1
                )),
                diff(X, axis=1)),
                axis=1)


    # Change for adm,

 #   parameters_vector_upsampled = parameters_vector_upsampled_code
    print "[II] Clustering features."
#
    features_clustering = GMM(n_components = n_clusters, covariance_type='diag')
#
    features_clustering.fit( features_vector_upsampled.T, y=parameters_vector_upsampled_code)
#
    features_clustering_means = features_clustering.means_
    features_clustering_labels = features_clustering.predict(features_vector_upsampled.T)
    features_clustering_sigmas = features_clustering.covars_
#
    features_vector_upsampled_estimated = zeros(shape(features_vector_upsampled))
#
#
    for n in range(0, len(features_vector_upsampled_estimated[0])):
        features_vector_upsampled_estimated[:,n] = features_clustering_means[features_clustering_labels[n]]
#
#
    print "[II] Features MSE for %d clusters: %.3f" % (n_clusters, mse(features_vector_upsampled, features_vector_upsampled_estimated))



    def happy_curve_classification(data, classes, estimator, Nd=1):
        print "[II] Generating Happy Curve "
        from copy import deepcopy
        estimator_fulldata = deepcopy(estimator)
        estimator_fulldata.fit(data, classes)
        labels = estimator.predict(data)

        # 1. Split data in two, training and testing

        Ntr =  int(round(data.shape[0]/2)) # Training dataset size
        Nts = data.shape[0] - Ntr     # Testing dataset size

        ratios = []

        for n in range(Nd, Ntr):
            train, test, trainl, testl = cross_validation.train_test_split(data, classes, test_size = 0.5, random_state = 0)
            train = train[0:n,:]
            trainl = trainl[0:n]
         #   print "trainl"
         #   print trainl
            estimator_ = deepcopy(estimator)
            estimator_.fit(train,trainl)
            labels = estimator_.predict(test)

            ratio = sum(array(testl==labels).astype(float))/len(labels)

            ratios.append(ratio)


        return ratios


    def cross_validate_classification(data, classes, estimator):
        print "[II] Crossvalidating... "
        from copy import deepcopy
        estimator_fulldata = deepcopy(estimator)
        estimator_fulldata.fit(data, classes)

        percents = arange(0.1,0.9,0.1)
        MSEs = []
        labels = estimator.predict(data)

        print "[II] for full training-testing: %.2f" %  (sum(array(classes==labels).astype(float))/len(labels))

        for p in percents:
            train,test,trainlabels,testlabels = cross_validation.train_test_split(data,classes,test_size=p,random_state=0)
            estimator_ = deepcopy(estimator)
            estimator_.fit(train, trainlabels)
            labels = estimator_.predict(test)
            print "[II] for training(%.2f)-testing(%.2f): %.2f" %  ((1-p),p,sum(array(testlabels==labels).astype(float))/len(labels))

        return MSEs

    def cross_validate_clustering(data, estimator):
        print "[II] Crossvalidating... "
        estimator_fulldata = estimator
        estimator_fulldata.fit(data)

     #   labels = estimator_fulldata.predict(data)
        means = estimator_fulldata.means_
    #    print means

        percents = arange(0.1,0.6,0.1)
        MSEs = []
        reconstructed = zeros(shape(data))
        labels = estimator.predict(data)
        for n in range(0, len(reconstructed)):
            reconstructed[n,:] = means[labels[n]]

        MSEs.append(mse(data,reconstructed))
        for p in percents:
            train,test = cross_validation.train_test_split(data,test_size=p,random_state=0)
            train = matrix(train)
            test = matrix(test)

            estimator.fit(train)
            means = estimator.means_
            labels = estimator.predict(test)
            reconstructed = zeros(shape(test))
            for n in range(0, len(reconstructed)):
                reconstructed[n,:] = means[labels[n]]

            m = mse(test,reconstructed)

            print "[II] MSE for clustering crossvalidated data %.2f-%.2f: %.5f" % ((1-p), p, m)
            MSEs.append(m)

            print "[II] Crossvalidation complete"

        return MSEs




    # Construct parameters alphabet, each symbol is going to be a different column vector
    # in parameter code matrix


    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_upsampled_code.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

   # state_seq, parameters_alphabet, parameters_alphabet_inv = vector_to_states(parameters_vector_upsampled_code)


  #  parameters_change_variable = matrix(diff_and_pad(parameters_vector_upsampled)!=0).astype(int)

   # changes_state_seq, changes_parameters_alphabet, changes_parameters_alphabet_inv = vector_to_states(parameters_change_variable)


    # This is an hmm that just codes the changes"
    # We have only two states, change and stay the same.

    # Uncomment that here

 #   parameters_vector_upsampled = parameters_vector_upsampled_code
  #  parameters_state_orig, parameter_state_alphabet_orig, parameter_state_alphabet_inv_orig = vector_to_states(parameters_vector)
    # TODO: HA
  #  parmaeters_vector_upsampled = parameters_vector_upsampled_code
    parameters_state, parameter_state_alphabet, parameter_state_alphabet_inv = vector_to_states(parameters_vector_upsampled)


    # Remove outliers

    ee = EllipticEnvelope()

 #   inliers = ee.fit(features_vector_upsampled.T).predict(features_vector_upsampled.T) == 1

    inoutliers = zeros((len(parameters_state),))
#
    for p in unique(parameters_state):
        ee = EllipticEnvelope(contamination=0.1)
        inoutliers[parameters_state == p] = ee.fit(features_vector_upsampled[:, parameters_state == p].T).predict(features_vector_upsampled[:, parameters_state == p].T)
#
#
    # Label inliers and outliers

    inliers = inoutliers > 0
    outliers = inoutliers < 0
#
#
#

    # Do pre-processing and re-label the outliers

    #completeavac = SinkHoleClassifier()
    #completeavac.fit(features_vector_upsampled[:, inliers].T, parameters_state[inliers])

    #parameters_state_original = parameters_state.copy()
    #parameters_state[outliers] = completeavac.predict(features_vector_upsampled[:,outliers].T)





   # features_vector_upsampled = features_vector_upsampled[:, inliers]
  #  parameters_state = parameters_state[inliers]

   # print "UNIQUE PARAMETERS STATE"
  # print unique(parameters_state)
#    asdasdasd
    print "[II] Testing Gaussian Naive Bayes Classifier"
    gnb = GaussianNB()
    gnb.fit(features_vector_upsampled.T, parameters_state)



    parameters_state_estimated = gnb.predict(features_vector_upsampled.T)

    output = states_to_vector(parameters_state_estimated, parameter_state_alphabet_inv)

    figure()
    subplot(211)
    plot(parameters_vector_upsampled.T)
    title('Parameter value')# % UpsamplingFactor)
    ylabel('value')
    xlabel('frame \#')
    subplot(212)
    #plot(smooth_matrix_1D(output).T)
    plot(output.T)
    ylabel('value')
    xlabel('frame \#')

    errorvector_nbc = [mse(output[:,i], parameters_vector_upsampled[:,i]) for i in range(0, len(parameters_state_estimated))]




    #cross_validate_classification(features_vector_upsampled.T, parameters_state, gnb)
   # hc = happy_curve_classification(features_vector_upsampled.T, parameters_state, gnb)

#
#    figure()
#    plot(hc)
#    figure()

    print "[II] Trying Multinomial HMM"

    # In order to do classification with HMMs, we need to:
    # 1. Split the parameters into classes
    # 2. Train one model per class
    # 3. Feed our data to all the models
    # 4. Check which has a better score and assig,n to M
    class SVMClassifier:
        def __init__(self,**kwargs):
            print "[II] SVM Classifier "
       #     self.clf = svm.SVC(**kwargs)
            self.name = "SVM"
            self.smallname = "svm"
            self.clf = multiclass.OneVsOneClassifier(svm.SVC(**kwargs) )

        def fit(self, X, classes):
            n_classes = max(unique(classes))+1

            self.clf.fit(X,classes)

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

        def getName(self):
            return self.name

        def cross_validate(self, data, classes):
            print "[II] SVN Classifier Crossvalidating... "
            from copy import deepcopy
            estimator = deepcopy(self)
            estimator_fulldata = deepcopy(self)
            estimator_fulldata.fit(data, classes)

            percents = arange(0.1,0.9,0.1)
            self.percents = percents * 100

            RATIOS = []
            labels = estimator.predict(data)


            print "[II] for full training-testing: %.2f" %  (sum(array(classes==labels).astype(float))/len(labels))

            for p in percents:
                traindata,testdata,trainlabels,testlabels = cross_validation.train_test_split(data,classes,test_size=p,random_state=1)


         #       print "[II] [%.2f-%.2f] Examples used for training: %s" % (p, (1-p), str(traindata))
          #      print "[II] [%.2f-%.2f] Examples used for testing: %s" % (p, (1-p), str(testdata))
                # Training phase



                estimator_ = deepcopy(estimator)
                estimator_.fit(traindata, trainlabels)

                predicted_labels = estimator_.predict(testdata)

                m = sum(array(predicted_labels==testlabels).astype(float))/len(testlabels)

                print "[II] [%.2f-%.2f] 5-fold cross-validation for SVN Classifier: %.1f" % (p,  (1-p), m)

                RATIOS.append(m)

            return RATIOS,percents







    class HmmClassifier2:
        def __init__(self, N=2, chain_size=4, n_components = 1):
            print "[II] HMM Classifier (%d states, %d components)" % (N, n_components)
            self.name = "HMM (%d states, %d components)" % (N, n_components)
            self.n_components = n_components
            self.chain_size = chain_size
            self.hmms_ = []
            self.N = N
            self.n_states = N

        def getName(self):
            return self.name

        def fit(self, features, parameters):
            self.n_params = len(unique(parameters))



            for n in range(0, self.n_params):
                hmm_ = hmm.GaussianHMM(n_components = self.n_components, covariance_type='full')

                # Get training data for each parameter

                for i in range(0, len(parameters)-self.chain_size):
                    seq = features[i:i+self.chain_size, :]
                    if parameters[i+self.chain_size-1] == n and seq.shape[0] == self.chain_size:
                 #       print "[II] Fitting: %s" % str(seq)

                        hmm_.fit([seq])

                self.hmms_.append(hmm_)

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

            scores = zeros((self.n_states, ))

            for i in range(0, len(labels)):
                testdata = X[i:i+self.chain_size,:]

                for j in range(0, self.n_states):
                    scores[j] = self.hmms_[j].score(testdata)
                labels[i] = argmax(scores)

            return labels

        def cross_validate(self, data, classes):
            print "[II] SVN Classifier Crossvalidating... "
            from copy import deepcopy
            estimator = deepcopy(self)
            estimator_fulldata = deepcopy(self)
            estimator_fulldata.fit(data, classes)

            percents = arange(0.1,0.9,0.1)
            self.percents = percents * 100

            RATIOS = []
            labels = estimator.predict(data)


            print "[II] for full training-testing: %.2f" %  (sum(array(classes==labels).astype(float))/len(labels))

            for p in percents:
                traindata,testdata,trainlabels,testlabels = cross_validation.train_test_split(data,classes,test_size=p,random_state=1)


             #   print "[II] [%.2f-%.2f] Examples used for training: %s" % (p, (1-p), str(traindata))
             #   print "[II] [%.2f-%.2f] Examples used for testing: %s" % (p, (1-p), str(testdata))

                # Training phase



                estimator_ = deepcopy(estimator)
                estimator_.fit(traindata, trainlabels)

                predicted_labels = estimator_.predict(testdata)

                m = sum(array(predicted_labels==testlabels).astype(float))/len(testlabels)

                print "[II] [%.2f-%.2f] 5-fold cross-validation for HMM2 Classifier: %.1f" % (p,  (1-p), m)

                RATIOS.append(m)

            return RATIOS,percents




    class HmmClassifier:
        def __init__(self, N=2, n_components = 1):
            print "[II] HMM Classifier (%d states, %d components)" % (N, n_components)
            self.name = "HMM (%d states, %d components)" % (N, n_components)
            self.n_components = n_components
            self.chain_size = N
            self.hmms_ = []
            self.N = N

        def getName(self):
            return self.name

        def fit(self, X, states):
            self.n_states = len(unique(states))

            for n in range(0,  self.n_states):
                hmm_ = hmm.GaussianHMM(n_components = self.n_components, covariance_type = 'full')

                # Get training data for each class
                vector = X[states == n,:]

                # Fit the HMM
               # print vector
                hmm_.fit([vector])

                # And append to the list
                self.hmms_.append(hmm_)

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

            m = 0

            scores = zeros((1, self.n_states))


            while m*N < X.shape[0]:
                if (m+1)*N > X.shape[0]:
                    testdata = X[m*N:,:]
                else:
                    testdata = X[m*N:(m+1)*N,:]

                  #  print testdata

                for i in range(0, self.n_states):
                    scores[0,i] = self.hmms_[i].score(testdata)

                if (m+1)*N > X.shape[0]:
                    labels[m*N:] = argmax(scores)
                else:
                    labels[m*N:(m+1)*N] = argmax(scores)

                m+= 1

            return labels

#        def cross_validate(self, data, classes, index_vector):
#            print "[II] Crossvalidating... "
#            from copy import deepcopy
#            estimator = deepcopy(self)
#            estimator_fulldata = deepcopy(self)
#            estimator_fulldata.fit(data, classes)
#
#            percents = arange(0.1,0.9,0.1)
#            self.percents = percents * 100
#            RATIOS = []
#            labels = estimator.predict(data)
#
#            print "[II] for full training-testing: %.2f" %  (sum(array(classes==labels).astype(float))/len(labels))
#
#            for p in percents:
#                train,test,trainlabels,testlabels = cross_validation.train_test_split(data,classes,test_size=p,random_state=1)
#                estimator_ = deepcopy(estimator)
#                estimator_.fit(train, trainlabels)
#                labels = estimator_.predict(test)
#                m = sum(array(testlabels==labels).astype(float))/len(labels)
#                RATIOS.append(m)
#                print "[II] for training(%.2f)-testing(%.2f): %.2f" %  ((1-p),p,sum(array(testlabels==labels).astype(float))/len(labels))
#
#            return RATIOS, percents

        def cross_validate(self, data, classes):
            print "[II] SVN Classifier Crossvalidating... "
            from copy import deepcopy
            estimator = deepcopy(self)
            estimator_fulldata = deepcopy(self)
            estimator_fulldata.fit(data, classes)

            percents = arange(0.1,0.9,0.1)
            self.percents = percents * 100

            RATIOS = []
            labels = estimator.predict(data)


            print "[II] for full training-testing: %.2f" %  (sum(array(classes==labels).astype(float))/len(labels))

            for p in percents:
                traindata,testdata,trainlabels,testlabels = cross_validation.train_test_split(data,classes,test_size=p,random_state=1)


       #         print "[II] [%.2f-%.2f] Examples used for training: %s" % (p, (1-p), str(traindata))
        #        print "[II] [%.2f-%.2f] Examples used for testing: %s" % (p, (1-p), str(testdata))
                # Training phase



                estimator_ = deepcopy(estimator)
                estimator_.fit(traindata, trainlabels)

                predicted_labels = estimator_.predict(testdata)
                m = sum(array(predicted_labels==testlabels).astype(float))/len(testlabels)


#                testdata = None
#
#                traindata = None
#                trainlabels = None
#                testlabels = None
#
#                for t in train:
#                    if traindata is None:
#                        traindata = data[index_vector == t, :]
#                        trainlabels = classes[index_vector == t]
#                    else:
#                        traindata = append(traindata, data[index_vector == t, :], 0)
#                        trainlabels = append(trainlabels, classes[index_vector ==t])
#
#                estimator_ = deepcopy(estimator)
#                estimator_.fit(traindata, trainlabels)
#
#
#                for t in test:
#                    if testdata is None:
#                        testdata = data[index_vector == t, :]
#                        testlabels = classes[index_vector == t]
#                    else:
#                        testdata = append(testdata, data[index_vector == t,:], 0)
#                        testlabels = append(testlabels, classes[index_vector == t])



            #    predicted_labels = estimator_.predict(testdata)

                m = sum(array(predicted_labels==testlabels).astype(float))/len(testlabels)

                print "[II] [%.2f-%.2f] 5-fold cross-validation for HMM Classifier: %.1f" % (p,  (1-p), m)

                RATIOS.append(m)

            return RATIOS,percents

    def cross_validate_full_k_fold(data, classes, list_of_classifiers, K=10):
        """ K-Folds Cross-validates all available classifiers, plots and saves barplots. """

        ratios = []
        for l in list_of_classifiers:
            ratio = l.k_fold_cross_validate(data, classes, K)
            ratios.append(ratio)

        L = len(ratios)

        ind = arange(L)

        fig, ax = subplots()
        rects = []



        width = 10/float(len(ratios))

        colors = ['r', 'y', 'g', 'c']

        labels_ = []

        for i in range(0, len(ratios)):
            rects.append(ax.bar(i*width,ratios[i],width,color=colors[i % len(colors)]))
            labels_.append(list_of_classifiers[i].getName())

        return ratios



    def cross_validate_full(seq, classes,  list_of_classifiers):
        """ Cross-validates all available classifiers, plots and saves barplots. """

        ratios = []
        percents = []
        for l in list_of_classifiers:
            ratios_, percents_ = l.cross_validate(seq, classes)
            ratios.append(ratios_)
            percents.append(percents_)


        L = len(ratios[0])

        ind = arange(L)

        fig,ax = subplots()

        rects = []

        W = 0.8
        width = W/float(len(ratios))

        colors = ['r', 'y', 'g', 'c']

        labels_ = []
        for i in range(0, len(ratios)):
            rects.append(ax.bar(ind+i*width,ratios[i],width,color=colors[i % len(colors)]))
            labels_.append(list_of_classifiers[i].getName())

        ax.legend(tuple(labels_))

        ax.set_xticks(ind+W/2)
        ax.set_xticklabels(tuple(percents[0]*100))
        ax.set_xlabel("Percentage % of test data")
        ax.set_ylabel("Success ratio")






      #  rects1 = ax.bar(ind,ratios[0],0.35,color='r')
      #  rects2 = ax.bar(ind+0.351,ratios[1],0.35,color='y')




        return ratios







#    hmmc = HmmClassifier(N = len(unique(parameters_state)), n_components = 1)
#    hmmc.fit(features_vector_upsampled.T, parameters_state)
#
#    hmmc2 = HmmClassifier(N = 12, n_components = 4)
#    hmmc2.fit(features_vector_upsampled.T, parameters_state)
#
    svmc = SVMClassifier(probability=True)
    svmc.fit(features_vector_upsampled.T, parameters_state)
#   # svmc.cross_validate(features_vector_upsampled.T, parameters_state, index_vector)
#
    nbc = NBClassifier()
    nbc.fit(features_vector_upsampled.T, parameters_state)
#
#
    #hchmmc = happy_curve_classification(features_vector_upsampled.T, parameters_state, hmmc, Nd=100)
#  #  cross_validate_classification(features_vector_upsampled.T, parameters_state, hmmc)
#
#    #hmmc.cross_validate(features_vector_upsampled.T, parameters_state)
#

    hmmc3 = HmmClassifier3(N=2, n_components = max(unique(parameters_state)))
    #print(parameters_state)
    plot(parameters_state)


    print(shape(features_vector_upsampled))
    obs = hmmc3.fit(features_vector_upsampled.T, parameters_state)



#    svmhmm = HmmSVMClassifier()
 #   svmhmm.fit(features_vector_upsampled.T, parameters_state)
   # ratios = cross_validate_full(features_vector_upsampled.T, parameters_state,   [nbc, svmc])
  #  ratios = cross_validate_full_k_fold(features_vector_upsampled.T, parameters_state,
   #                                     [hmmc3])

    def approximate_chain_size(data,parameters):
        ratios = []

       # chainsizes = arange(1, 40)
       # for cs in chainsizes:
       #     print "[II] Trying chain size: %d" % cs
       #     hmmc = HmmClassifier3(N=cs, n_components = 2)
       #     ratios.append(hmmc.k_fold_cross_validate(data,parameters,K=10))


        componentssize = arange(1, 12)

        for cs in componentssize:
            print "[II] Trying component size: %d" % cs
            hmmc = HmmClassifier3(N=2, n_components = cs)
            ratios.append(hmmc.k_fold_cross_validate(data, parameters, K=2))




        figure()
        plot(componentssize, ratios)
        xlabel('Chain Size')
        ylabel('Success Ratio')
        title('10-Fold cross validation success ratio vs chain size')


        return ratios

  #  ratios = approximate_chain_size(features_vector_upsampled.T, parameters_state)

    #predictedhmmc2 = hmmc2.predict(features_vector_upsampled.T)

  #  figure()





#    hmms_ = []
#
#    for n in range(0, len(parameter_state_alphabet)):
#        #hmm_ = hmm.GMMHMM(n_components = 1, n_mix = 2)
#        hmm_ = hmm.GaussianHMM(n_components = 1,covariance_type = 'full')
#
#        # Get training data for each class
#        vector = features_vector_upsampled[:,parameters_state == n]
#
#        #if vector.shape[1] < n_clusters:
#        #    hmms_.append(None)
#        #else:
#
#        hmm_.fit([vector.T])
#
#        # Append to the list
#
#        hmms_.append(hmm_)
#
#    labels = zeros((features_vector_upsampled.shape[1],))
#
#    N = 20
#    m = 0
#
#    while m*N < features_vector_upsampled.shape[1]:
#
#        scores = zeros((1, len(parameter_state_alphabet)))
#
#        if (m+1)*N >  features_vector_upsampled.shape[1]:
#            testdata = features_vector_upsampled[:,m*N:]
#        else:
#            testdata = features_vector_upsampled[:,m*N:(m+1)*N]
#
#        for i in range(0, len(parameter_state_alphabet)):
#            if hmms_[i] is not None:
#                scores[0,i] = hmms_[i].score(testdata.T)
#            else:
#                scores[0,i] = -100000 # Very large negative score
#        if (m+1)*N >= features_vector_upsampled.shape[1]:
#            labels[m*N:] = argmax(scores)
#        else:
#            labels[m*N:(m+1)*N] = argmax(scores)
#
#        m += 1


   # figure()
    #plot(labels.T)


#    labels = hmmc.predict(features_vector_upsampled.T)
#    estimated = states_to_vector(labels,parameter_state_alphabet_inv)
#    plot(estimated.T,'r--')
#
#    title('Estimated parameter values')
#    legend(['Naive Bayes Classifier', 'HMM chain size %d (%.1fms)' % (N, float(N)/UpsamplingFactor*23.0)])
#
#    ylabel('value')
#    xlabel('frame #')
#

   # close('all')

   # plot(features_clustering_labels/float(max(features_clustering_labels)))
  #  plot(parameters_vector_upsampled.T/max(ravel(parameters_vector_upsampled)))


    def plot_clusters(x, labels):
        figure()
        symbols = ['>', 'x', '.', '<','v']
        colors =  ['b', 'r', 'g', 'm','c']

        for r in range(0, x.shape[0]):
            scatter(x[r,0], x[r,1], c=colors[int(labels[r]) % len(colors)], marker=symbols[int(labels[r]) % len(symbols)])


  #  plot_clusters(features_vector_upsampled.T, parameters_state)
  #

    # SVM

    class HmmClassifierSVN_Naive:
        def __init__(self, N=2, n_components = 1):
            self.n_components = n_components
            self.chain_size = N
            self.hmms_ = []
            self.N = N

        def fit(self, X, states):
            self.n_states = len(unique(states))

            for n in range(0,  self.n_states):
                hmm_ = hmm.GaussianHMM(n_components = self.n_components, covariance_type = 'full')

                # Get training data for each class
                vector = X[states == n,:]

                # Fit the HMM
               # print vector
                hmm_.fit([vector])

                # And append to the list
                self.hmms_.append(hmm_)

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

            m = 0

            scores = zeros((1, self.n_states))


            while m*N < X.shape[0]:
                if (m+1)*N > X.shape[0]:
                    testdata = X[m*N:,:]
                else:
                    testdata = X[m*N:(m+1)*N,:]

                  #  print testdata

                for i in range(0, self.n_states):
                    scores[0,i] = self.hmms_[i].score(testdata)

                if (m+1)*N > X.shape[0]:
                    labels[m*N:] = argmax(scores)
                else:
                    labels[m*N:(m+1)*N] = argmax(scores)

                m+= 1

            return labels




    def plot_parameters_estimation(list_of_estimators):
        for e in list_of_estimators:
            name_ = e.getName()



            fig = figure()
            fig.suptitle('Parameters estimation using %s' % name_, fontweight='bold')
            subplot(311)
            title ('original parameters')
            param_real = states_to_vector(parameters_state,parameter_state_alphabet_inv).T
            plot(param_real)
            xlabel('sample')
            ylabel('normalized parameter value')
            legend(parameter_captions)
            subplot(312)
            title('estimated parameters')
            pred = e.predict(features_vector_upsampled.T)
            param_est = states_to_vector(pred,parameter_state_alphabet_inv).T
            plot(param_est)
            xlabel('sample')
            ylabel('normalized parameter value')
            legend(parameter_captions)
            subplot(313)
            error_real_est = zeros((len(param_est),))
            for i in range(0, len(error_real_est)):
                error_real_est[i] = mse(param_real[i],param_est[i])

            totalmse = sum(error_real_est)
            title('mean squared error (total: %.3f)' % totalmse)

            plot(error_real_est)
            xlabel('sample')
            ylabel('normalized parameter value')



#    hmm2 = hmm.MultinomialHMM(n_components = n_classes, startprob = prior_p_q, transmat = transmat )



    # TODO: Embed  SVM emissions to HMM. Try to do pick maximum likelihood sequence from
    # classifiers like in HmmClassifier3. For every class, train an SVM, with n_components as
    # states, and compute the log-likelihood for prediction. Should work.



#    class HMMsvmClassifier:
#        def __init__(self, N=6):
#            print "[II] HMM-SVM Classifier (%d time steps)" % N
#            self.name = "HMM-SVM (%d time steps)" % N
#          #  self.n_components = n_components
#            self.chain_size = N
#            self.svm = MyAVAClassifier()
#
#
#        def getName(self):
#            return self.name
#
#        def fit(self, features, parameters):
#
#            # First train emission svm
#
#            self.svm.fit(features, parameters)
#
#         #   print parameters
#            n_classes = max(unique(parameters)) + 1
#
#            print "[II] Number of classes: ", n_classes
#            hmms = [None]*n_classes
#            chain_size = self.chain_size
#            params = [None]*n_classes
#            obs = [None]*n_classes
#
#            for i in range(chain_size, len(parameters)):
#                class_ = parameters[i-1]
#                params_ = parameters[i-chain_size:i]
#                feats_ =features[i-chain_size:i,:]
#
#                if params[class_] is None:
#                    params[class_] = [params_]
#                    obs[class_] = [feats_]
#                else:
#                    params[class_].append(params_)
#                    obs[class_].append(feats_)
#
#
#
#            for i in range(0, len(params)):
#                if params[i] is not None and len(params[i]) != 0:
#                    hmm_ = HMMsvm(self.svm)
#              #      print "[II] %d Fitting params: " % i, params[i]
#                    print "[!!] OBSERVATIONS: "
#                    print obs[i]
#                    hmm_.fit(obs[i], params[i])
#
#                    print "[!!] PARAMETERS: "
#                    print params[i]
#
#                    # TODO: Left here on 06/07 19:56
#
#                  #  hmm_.fit(features,)
#             #       print "obs: ", obs[i]
#            #        print "params: ", params[i]
#             #       hmm_.fit(obs[i], params[i],flr=1e-9)
#                    hmms[i] = hmm_
#
#            self.hmms = hmms
#
#            return obs
#        def predict(self, features, mfilt=20):
#            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:
#                        q,p = hmms[j].logviterbi(seq)
#
#                        scores[j] = logsumexp(p)
#                    else:
#                        scores[j] = -infty
#
#                predicted_classes[i] = argmax(scores)
#
#            # Do a median filter at the end
#
#      #      predicted_classes = median_filter(predicted_classes,mfilt)
#
#            return predicted_classes

    class HmmSvmClassifier3:
        def __init__(self, N=2,n_components = 1):
            print "[II] HMM/SVM Classifier (%d states, %d components)" % (N, n_components)
            self.name = "HMM (%d time steps, %d components)" % (N, n_components)
            self.n_components = n_components
            self.smallname = "hmmsvmc"
            self.chain_size = N

        def getName(self):
            return self.name

        def fit(self, features, parameters):



         #   print parameters
            n_classes = max(unique(parameters)) + 1

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

            print "[II] Number of classes: ", n_classes
            hmms = [None]*n_classes
      #      svms = [None]*n_classes
            chain_size = self.chain_size
            obs = [None]*n_classes
            cls = [None]*n_classes

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

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



            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='full')
                    svm_ = MyAVAClassifier()
                    print obs[i]
                    print cls[i]

                    O = obs[i][0]
                    C = cls[i][0]

                    for j in range(0, len(obs[i])):
                        O = append(O, obs[i][j], 1)
                        C = append(C, cls[i][j], 1)

                    print O
                    print C

                    svm_.fit(O.T,C)
                    hmm_ = HMMsvm(svm_)
                    hmm_.fit(obs[i],cls[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:
                        print hmms[j].score(seq)
                        scores[j] = hmms[j].score(seq)[j]
                    else:
                        scores[j] = -infty

                predicted_classes[i] = argmax(scores)

            # Do a median filter at the end

      #      predicted_classes = median_filter(predicted_classes,mfilt)

            return predicted_classes






             #  newstates = self.logviterbi()

#
#
#            # Alphas and Betas
#            X = features
#            alphas = zeros((T,n_classes))
#            betas = zeros((T,n_classes))
#
#            forward_messages = zeros((T,n_classes))
#            backward_messages = zeros((T, n_classes))
#
#            print "[II] Computing alpha, beta..."
#            for t in range(1, T+1):
#                forward_messages[t-1,:] = self._forward_message(forward_messages[t-2,:],X[0:t,])
#                backward_messages[-t,:] = self._backward_message(backward_messages[-t+1], X[-t:,:])
#
#
#            print "[II] Computing ksu..."
#
#            a_norm = forward_messages[-1,n_classes-1] # alpha_T
#
#            ksu = zeros((T, n_classes, n_classes))
#
#            for i in range(0, n_classes):
#                for j in range(0, n_classes):
#                    for t in range(0, T-1):
#                        ksu[t,i,j] = exp(forward_messages[t, i] + log(transmat[i,j]) + log(self.estimate_emission_probability(X[t,:],j)) + backward_messages[t+1,j] - a_norm)
#
#
#
#            self.alphas = forward_messages
#            self.betas = backward_messages
#            self.ksu = ksu
#
#            transmat_new = zeros((n_classes,n_classes))
#
#            for i in range(0, n_classes):
#
#                for j in range(0, n_classes):
#                    transmat_new[i,j] = sum(ksu[:,i,j])/sum(ksu[:,i,:])
#
#            self.transmat_new = transmat_new
#











#
#        def forward(self, X):
#            # Compute log likelihood using the forward algorithm
#
#            # Number of states
#            N = self.n_classes
#
#            # Number of Observations
#            T = X.shape[0]
#
#
#            transmat = self.transmat
#
#            # Initialization
#         #   a1 = ones((N,))/N
#            a1 = self.prior
#
#            alpha = zeros((N,))
#            for i in range(0, N):
#                alpha[i] = log(a1[i])+log(self.estimate_emission_probability(X[0,:], i))
#          #  print alpha
#            # Recursion
#            for t in range(0, T):
#                alpha_old = alpha.copy()
#                x = X[t, :]
#                for j in range(0, N):
#                    tosum = zeros((N,))
#                    for i in range(0, N):
#                       tosum[i] = alpha_old[i]+log(transmat[i,j])
#
#                  #  Sum = logsum(tosum)
#                    Sum = logsumexp(tosum)
#                    bj = self.estimate_emission_probability(x, j)
#
#                    alpha[j] = Sum+log(bj)
#            #    print alpha
#
#            tosum = zeros((N,))
#            for i in range(0, N):
#                tosum[i] = alpha[i] + log(transmat[i,N-1])
#
#            return tosum
#
#        def backward(self, X):
#            # Number of states
#            N = self.n_classes
#
#            # Number of Observations
#            T = X.shape[0]
#
#            transmat = self.transmat
#
#            # Initialization
#
#            b1 = ones((N,))/N
#
#            beta = zeros((N,))
#            for i in range(0, N):
#                beta[i] = log(b1[i])+log(self.estimate_emission_probability(X[-1,:], i))
#
#            for t in range(0, T):
#                beta_old = beta.copy()
#                x = X[-t, :]
#                for j in range(0, N):
#                    tosum = zeros((N, ))
#                    for i in range(0, N):
#                        tosum[i] = beta_old[i]+log(transmat[i,j])
#
#                    Sum = logsumexp(tosum)
#                    bj = self.estimate_emission_probability(x, j)
#                    beta[j] = Sum+log(bj)
#
#            tosum = zeros((N,))
#
#            for i in range(0, N):
#                tosum[i] = beta[i] + log(transmat[i,0])
#
#        #3    p = logsumexp(tosum)
#
#            return tosum


        def _log_likelihood(self, X):

            return logsumexp(self.forward(X))





        def _likelihood(self, X):
            # Number of states
            N = self.n_classes

            # Number of Observations
            T = X.shape[0]


            transmat = self.transmat

            # Initialization
           # a1 = ones((N,))/N
            a1 = self.prior

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

            # Recursion
            print alpha
            for t in range(1, T):
                alpha_old = alpha.copy()
                x = X[t, :]
                for j in range(0, N):
                    Sum = 0
                    for i in range(0, N):
                        Sum += alpha_old[i]*transmat[i,j]

                    alpha[j] = Sum*self.estimate_emission_probability(x, j)
                print alpha


            # Termination

            Sum = 0
            for i in range(0, N):
                Sum += alpha[i]*transmat[i, N-1]

            p = Sum

            return p







#        def fit(self, X, states):
#         #   print parameters
#            n_classes = max(unique(states)) + 1
#
#      #      svms = [None]*n_classes
#            obs = [None]*n_classes
#            sta = [None]*n_classes
#            chain_size = self.chain_size
#
#         #22   obs = None
#        #    sta = None
#
#            print "[II] Number of classes: ", n_classes
#
#            # For each class:
#                # Concatenate examples
#                # Fit SVM
#
#            for i in range(chain_size, len(states)):
#                class_ = states[i-1]
#                seq = X[i-chain_size:i, :]
#                states_ = states[i-chain_size:i]
#
#
#                if obs[class_] is None:
#                    obs[class_] = [seq]
#                    sta[class_] = [states_]
#                    self.svms.append(MyAVAClassifier())
#                else:
#                    obs[class_].append(seq)
#                    sta[class_].append(states_)
#
#
#            for i in range(0, len(obs)):
#
#                o = None
#                s = None
#
#                for j in range(0, len(obs[i])):
#                    if o is None:
#                        o = obs[i][j]
#                        s = sta[i][j]
#
#                    else:
#                        o = append(o, obs[i][j],0)
#                        s = append(s, sta[i][j])
#
#
#                self.svms[i].fit(o, s)

#        def predict(self, features):
#            chain_size = self.chain_size
#            svms = self.svms
#
#            predicted_classes = zeros((features.shape[0],))
#            for i in range(chain_size, features.shape[0]):
#                scores = zeros((len(svms)))
#
#                seq = features[i-chain_size:i, :]
#
#                for j in range(0, len(svms)):
#                    if svms[j] is not None:
#                        scores[j] = svms[j].compute_log_likelihood(seq)
#                    else:
#                        scores[k] = -infty
#                predicted_classes[i] = argmax(scores)
#
#            return predicted_classes





            # Marginalize over the latent variable
#            for i in range(0, X.shape[0]):
#                P = zeros((self.n_classes,))
#                x = X[i,:]
#                for j in range(0, len(P)):
#                    P[j] = self.estimate_emission_probability(j, x)
#
#                S[i] = log(sum(P[j]))
#
#            return sum(S)




    # Continue from here
    X = features_vector_upsampled.T
    y = parameters_state

   # clf = svm.SVC()
#    clf.fit(X,y)


 #   parameters_state_y = clf.predict(X)


    predhmmc3 = hmmc3.predict(features_vector_upsampled.T)

    myava = MyAVAClassifier()
    myava.fit(features_vector_upsampled.T, parameters_state)

    predmyava = myava.predict(features_vector_upsampled.T)

 #   hmmsvmc = HMMsvmClassifier(N=2)
#    hmmsvmc.fit(features_vector_upsampled.T, parameters_state)
    #predhmmsvmc = hmmsvmc.predict(features_vector_upsampled.T)

    obsc = MyAVAClassifier()
    obsc.fit(features_vector_upsampled.T, parameters_state)
 #   hmm2 = hmm.GaussianHMM(n_components=6)

  #  hmm2.fit([features_vector_upsampled.T])
  #  hmm3 = HMMsvm(obsc)
   # hmm3.fit([features_vector_upsampled.T],[parameters_state])
   # predhmmsvm = hmm3.predict(features_vector_upsampled.T)

#    hmmsvmc = HMMsvmClassifier()

    hmmsvmc = HMMsvmClassifier()
  #  hmmsvmc = HmmSvmClassifier3()
    hmmsvmc.fit(features_vector_upsampled.T, parameters_state)
    predhmmsvm = hmmsvmc.predict(features_vector_upsampled.T)
    plot_parameters_estimation([nbc, myava, hmmc3, hmmsvmc])


    #plot(parameters)
    figure()
    plot(parameters_state,'b')

  #  plot(parameters_state_y,'r--')
    plot(predhmmc3,'g+-')

    plot(predhmmsvm, 'co-')

 #   plot(predhmmsvmc, 'bo-')

    legend(['Real', 'SVM', 'HMM','HMM-SVM'])

    xlabel('sample')
    ylabel('state')
    # TODO, HMM with SVN, Cross-validation

    # Using hidden markov svm

    svmhmm = ""

   # print "[II] Success ratio for SVM-RBF Classifier: %.2f" % (sum(parameters_state==parameters_state_y).astype(float)/float(len(parameters_state)))

    print "[II] All-vs-All custom Support Vector Classifier Success Ratio: %.2f " %(sum(predmyava==parameters_state).astype(float)/len(parameters_state))
    print "[II] Success ratio for HMM: %.2f" % (sum(parameters_state==predhmmc3).astype(float)/float(len(predhmmc3)))
    print "[II] Success ratio for HMM-SVM: %.2f" % (sum(parameters_state==predhmmsvm).astype(float)/float(len(predhmmsvm)))





    model_ghmm = ReverbModel("GaussianHMM Classifier Model", kernel, q, features_to_keep, hmmc3, parameter_dictionary=parameter_state_alphabet_inv, moments_vector=moments_vector)
    model_gnb = ReverbModel("Naive Bayes Classifier Model", kernel, q, features_to_keep, nbc, parameter_dictionary=parameter_state_alphabet_inv, moments_vector=moments_vector  )
    model_svm = ReverbModel("AVA LinearSVM Classifier Model", kernel, q, features_to_keep, myava, parameter_dictionary=parameter_state_alphabet_inv, moments_vector=moments_vector  )
    model_hmmsvm = ReverbModel("HMM0SVM Classifier Model", kernel, q, features_to_keep, hmmsvmc, parameter_dictionary=parameter_state_alphabet_inv, moments_vector=moments_vector  )

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



    def split_data_to_k_samples(data, params, K):
        L = shape(data)[0]
        M = int(round(L/float(K)))

        samples_data = []
        samples_parameters = []
        for k in range(K):
            if (k+1)*M > L:
                samples_data.append(data[k*M:,:])
                samples_parameters.append(params[k*M:])
            else:
                samples_data.append(data[k*M:(k+1)*M,:])
                samples_parameters.append(params[k*M:(k+1)*M])

        return samples_data, samples_parameters

    def join_data(data, params):
        L = 0
        for p in params:
            L += len(p)

        new_data = zeros((L, data[0].shape[1]))
        new_params = zeros((L,))

        idx = 0
        for n in range(0,len(data)):
            new_data[idx:idx+data[n].shape[0],:] = data[n]
            new_params[idx:idx+len(params[n])] = params[n]
            idx += len(params[n])

        return new_data, new_params

    cross_validation_ratios = {}
    cross_validation_mses = {}


    features_samples, parameters_samples = split_data_to_k_samples(features_vector_upsampled.T, parameters_state, 5)

    from sklearn.cross_validation import KFold

    from sklearn import metrics



    KF = 10

    s1,s2 = split_data_to_k_samples(features_vector_upsampled.T, parameters_state, KF)
    r1,r2 = join_data(s1,s2)

    kf = KFold(KF, n_folds=KF) # Leave one out

    trial = 1.0

    rhmmc = 0
    rgnbc = 0
    ravac = 0
    rhmmsvmc = 0


    msegbnc = 0
    msehmmc = 0
    mseavac = 0
    msehmmsvmc = 0

    fullratios = []
    fullmses = []


    def calculate_mse(pred, real):
        realparams = states_to_vector(real, parameter_state_alphabet_inv)
        estparams = states_to_vector(pred, parameter_state_alphabet_inv)
        return mse(realparams, estparams)
    ratios_gnb = []
    ratios_hmm = []
    ratios_svm = []
    ratios_hmmsvm = []
    ratios_sinkhole = []

    mse_gnb = []
    mse_hmm = []
    mse_svm = []
    mse_hmmsvm = []
    mse_sinkhole = []


    for train, test in kf:

        train_idx = 100

        print train,test
        f1 = []
        p1 = []

        tf1 = []
        tp1 = []

        for i in train:
            f1.append(s1[i-1])
            p1.append(s2[i-1])

        for i in test:
            tf1.append(s1[i-1])
            tp1.append(s2[i-1])

        new_trainset_data, new_trainset_params = join_data(f1, p1)
        new_trainset_params = new_trainset_params.astype(int)


        # Detect and remove outliers

#        ee = EllipticEnvelope()

#        ee.fit(new_trainset_data, new_trainset_params)


#        inliers = ee.predict(new_trainset_data) == 1

#        new_trainset_data = new_trainset_data[inliers, :]
#        new_trainset_params = new_trainset_params[inliers]


        new_testset_data, new_testset_params = join_data(tf1, tp1)
        new_testset_params = new_testset_params.astype(int)



        param_est = states_to_vector(new_testset_params,parameter_state_alphabet_inv).T







        X = arange(50, len(new_trainset_params), 50)
        print X
    #    X = append(X, [len(new_trainset_params)], 1)

        
#
        for M in X:
#                    print CHAINSIZE
 #                   print NCOMPONENTS
                    hmmc = HmmClassifier3(N=CHAINSIZE, n_components = NCOMPONENTS)
                    gnbc = GaussianNB()
                    avac = MyAVAClassifier()
                    hmmsvmc = HMMsvmClassifier()

                    gnbc.fit(new_trainset_data[0:M,:], new_trainset_params[0:M])
                    pgnbc = gnbc.predict(new_testset_data)
                    rgnbc = metrics.accuracy_score(pgnbc, new_testset_params)
                    mgnbc = calculate_mse(pgnbc,new_testset_params)

#
                    hmmc.fit(new_trainset_data[0:M,:], new_trainset_params[0:M])
                    phmmc = hmmc.predict(new_testset_data)
                    rhmmc = metrics.accuracy_score(phmmc, new_testset_params)
                    mhmmc = calculate_mse(phmmc,new_testset_params)

#
                    try:
                        avac.fit(new_trainset_data[0:M,:], new_trainset_params[0:M])
                        pavac = avac.predict(new_testset_data)
                        ravac = metrics.accuracy_score(pavac, new_testset_params)
                        mavac = calculate_mse(pavac,new_testset_params)
                    except:
                        ravac = 0
                        mavac = inf
#
#
#
#
#
                    try:
                        hmmsvmc.fit(new_trainset_data[0:M,:], new_trainset_params[0:M])
                        phmmsvm = hmmsvmc.predict(new_testset_data)
                        rhmmsvm = metrics.accuracy_score(phmmsvm, new_testset_params)
                        mhmmsvm = calculate_mse(phmmsvm,new_testset_params)
#
                    except:
                        rhmmsvm = 0
                        mavac = inf
#
#                    ratios_gnb.append(rgnbc)
#                    ratios_hmm.append(rhmmc)
#                    ratios_svm.append(ravac)
#                    ratios_hmmsvm.append(rhmmsvm)
#
#                    mse_gnb.append(mgnbc)
#                    mse_hmm.append(mhmmc)
#                    mse_svm.append(mavac)
#                    mse_hmmsvm.append(mhmmsvm)
#
  #      fullratios.append((trial,X,ratios_gnb,ratios_hmm,ratios_svm,ratios_hmmsvm))
  #      fullmses.append((trial,X,mse_gnb,mse_hmm,mse_svm,mse_hmmsvm))
        
#
#
#
#






  #      hmmc = HmmClassifier3(N=3, n_components = max(unique(parameters_state))+1)
        hmmc = HmmClassifier3(N=CHAINSIZE, n_components = NCOMPONENTS)
        gnbc = GaussianNB()
        avac = MyAVAClassifier()
        hmmsvmc = HMMsvmClassifier()
        sinkholec = SinkHoleClassifier()

        hmmc.fit(new_trainset_data, new_trainset_params)
        gnbc.fit(new_trainset_data, new_trainset_params)
#        try:
        avac.fit(new_trainset_data, new_trainset_params)
        sinkholec.fit(new_trainset_data, new_trainset_params)
        pavac = avac.predict(new_testset_data)


#        except:
#            pavac = ones((new_testset_data.shape[0],))*new_trainset_params[0]


 #       try:
        hmmsvmc.fit(new_trainset_data, new_trainset_params)
        phmmsvm = hmmsvmc.predict(new_testset_data)
#        except:
#
#        #    phmmsvm = hmmsvmc.predict(new_testset_data)
#            phmmsvm = ones((new_testset_data.shape[0],))*new_trainset_params[0]


        phmmc = hmmc.predict(new_testset_data)
        pgnbc = gnbc.predict(new_testset_data)
        psinkhole = sinkholec.predict(new_testset_data)

        print "[II] (Trial %d) Classification ratio for GNB: %.2f" % (trial, metrics.accuracy_score(pgnbc, new_testset_params))
        print "[II] (Trial %d) Classification ratio for SVM: %.2f" % (trial, metrics.accuracy_score(pavac, new_testset_params))
        print "[II] (Trial %d) Classification ratio for HMM: %.2f" % (trial, metrics.accuracy_score(phmmc, new_testset_params))
        print "[II] (Trial %d) Classification ratio for HMM/SVM: %.2f" % (trial, metrics.accuracy_score(phmmsvm, new_testset_params))
        print "[II] (Trial %d) Classification ratio for Sinkhole Approach: %.2f" % (trial, metrics.accuracy_score(psinkhole, new_testset_params))


        mse

      #  rgnbc = rhmmc*(trial-1)/trial +  metrics.accuracy_score(pgnbc, new_testset_params)/trial
      #  ravac = rhmmc*(trial-1)/trial +  metrics.accuracy_score(pavac, new_testset_params)/trial
      #  rhmmc = rhmmc*(trial-1)/trial +  metrics.accuracy_score(phmmc, new_testset_params)/trial
      #  rhmmsvmc = rhmmc*(trial-1)/trial +  metrics.accuracy_score(phmmsvm, new_testset_params)/trial

        rgnbc = metrics.accuracy_score(pgnbc, new_testset_params)
        ravac = metrics.accuracy_score(pavac, new_testset_params)
        rhmmc = metrics.accuracy_score(phmmc, new_testset_params)
        rhmmsvmc = metrics.accuracy_score(phmmsvm, new_testset_params)
        rsinkhole = metrics.accuracy_score(psinkhole, new_testset_params)

        ratios_gnb.append(rgnbc)
        ratios_svm.append(ravac)
        ratios_hmm.append(rhmmc)
        ratios_hmmsvm.append(rhmmsvmc)
        ratios_sinkhole.append(rsinkhole)





     #   msegbnc = msegbnc*(trial-1)/trial +  calculate_mse(pgnbc,new_testset_params)/trial
      #  msehmmc = msehmmc*(trial-1)/trial +  calculate_mse(phmmc,new_testset_params)/trial
      #  mseavac = mseavac*(trial-1)/trial +  calculate_mse(pavac,new_testset_params)/trial
      #  msehmmsvmc = msehmmsvmc*(trial-1)/trial +  calculate_mse(phmmsvm,new_testset_params)/trial

        msegnb = calculate_mse(pgnbc,new_testset_params)
        msesvm = calculate_mse(pavac,new_testset_params)
        msehmm = calculate_mse(phmmc,new_testset_params)
        msehmmsvm = calculate_mse(phmmsvm,new_testset_params)
        msesinkhole = calculate_mse(psinkhole,new_testset_params)

        mse_gnb.append(msegnb)
        mse_svm.append(msesvm)
        mse_hmm.append(msehmm)
        mse_hmmsvm.append(msehmmsvm)
        mse_sinkhole.append(msesinkhole)






        trial += 1

    print ratios_gnb
    

    sucrats = [mean(ratios_gnb),mean(ratios_svm), mean(ratios_hmm),mean(ratios_hmmsvm), mean(ratios_sinkhole)]
    sucratsstd = [std(ratios_gnb,ddof=1),std(ratios_svm,ddof=1), std(ratios_hmm,ddof=1),std(ratios_hmmsvm,ddof=1), std(ratios_sinkhole,ddof=1)]
    mserats = [mean(mse_gnb),mean(mse_svm), mean(mse_hmm),mean(mse_hmmsvm), mean(mse_sinkhole)]
    #mserats = [msegbnc, mseavac, msehmmc,msehmmsvmc]

    print tuple(sucrats)
#    print tuple(sucratsstd)


#    scores = zeros((len(sucrats), len(sucrats)))
#    for i in range(0, len(sucrats)):
#        for j in range(0, len(sucrats)):
#            scores[i,j] = (sucrats[i] - sucrats[j])/sqrt(sucratsstd[i]**2 + sucratsstd[j]**2)
    print tuple(mserats)
#    print scores



#    print (msegbnc, mseavac, msehmmc,msehmmsvmc )

    sample_sizes = X/float(max(X))
    gnbc_ratios_curves = zeros((len(X),len(fullratios)))
    gnbc_mses_curves = zeros((len(X),len(fullratios)))
#
    avac_ratios_curves = zeros((len(X),len(fullratios)))
    avac_mses_curves = zeros((len(X),len(fullratios)))
#
    hmmc_ratios_curves = zeros((len(X),len(fullratios)))
    hmmc_mses_curves = zeros((len(X),len(fullratios)))
#
    hmmsvmc_ratios_curves = zeros((len(X),len(fullratios)))
    hmmsvmc_mses_curves = zeros((len(X),len(fullratios)))
#
#
    for i in range(0, len(fullratios)):
        for j in range(0, len(X)):
            gnbc_ratios_curves[j,i] = fullratios[i][2][j]
            gnbc_mses_curves[j,i] = fullmses[i][2][j]

            avac_ratios_curves[j,i] = fullratios[i][3][j]
            avac_mses_curves[j,i] = fullmses[i][3][j]

            hmmc_ratios_curves[j,i] = fullratios[i][4][j]
            hmmc_mses_curves[j,i] = fullmses[i][4][j]

            hmmsvmc_ratios_curves[j,i] = fullratios[i][5][j]
            hmmsvmc_mses_curves[j,i] = fullmses[i][5][j]
#
#
#
#
#
    gnbc_mses_curve = mean(gnbc_mses_curves,1)
    avac_mses_curve = mean(avac_mses_curves,1)
    hmmc_mses_curve = mean(hmmc_mses_curves,1)
    hmmsvmc_mses_curve = mean(hmmsvmc_mses_curves,1)
#
    gnbc_ratios_curve = mean(gnbc_ratios_curves,1)
    avac_ratios_curve = mean(avac_ratios_curves,1)
    hmmc_ratios_curve = mean(hmmc_ratios_curves,1)
    hmmsvmc_ratios_curve = mean(hmmsvmc_ratios_curves,1)
#
    figure()
    subplot(211)
    plot(sample_sizes,gnbc_ratios_curve)
    plot(sample_sizes,avac_ratios_curve)
    plot(sample_sizes,hmmc_ratios_curve)
    plot(sample_sizes,hmmsvmc_ratios_curve)
    xlabel('Part of training set used')
    ylabel('Accuracy')
    legend(['GNB','SVM','HMM','HMM/SVM'])
#
    subplot(212)
    plot(sample_sizes,gnbc_mses_curve)
    plot(sample_sizes,avac_mses_curve)
    plot(sample_sizes,hmmc_mses_curve)
    plot(sample_sizes,hmmsvmc_mses_curve)
    xlabel('Part of training set used')
    ylabel('Mean Squared Error')
#
    legend(['GNB','SVM','HMM','HMM/SVM'])
#
    show()
#   # tight_layout()
    savefig('ratios-2.png',dpi=500)
#
#
#