<commit_before>try:
    import psyco
    psyco.full()
except ImportError: pass


__doc__ = """
The present module Generative Linear Gaussian Models (GLGM) implements in Python programming language
the following models based mainly on references nr.1 (see References below):

        x(t+1) = Ax(t) + w(t) = Ax(t) + w0,     w0 = N(0, Q)
        y(t)   = Cx(t) + v(t) = Cx(t) + v0,     v0 = N(0, R)
     where:
        x is a (k,1) vector of latent factors or state or causes (hidden) variables
        y is a (p,1) vector of observations
        A is a (k,k) matrix of states transition
        C is a (p,k) matrix of observation measurement or generative matrix
        w is a (k,1) vector of state evolution White Gaussian Noise (WGN)
        v is a (k,1) vector of observation evolution WGN
        w and v are statistically independent of each other and of x and y
        N stands for for Gaussian or Normal probability density function (pdf)

In particular the module aims to implement the following:
    - fa    (Factor Analysis)
    - ppca  (Probabilistic Principal Component Analysis)
    - pca   (Principal Component Analysis)
    - mog   (Mixture of Gaussians)
    - vq    (Vector Quantization)
    - k-means clustering
    - mofa  (Mixture of Factor Analyzers)
    - ica   (Independent Component Analysis)

References:
    1  Unifying Review of Linear Gaussian Models
        Roweis, Ghahramani
        Neural Computation 1999 11:2, 305-345.
    2  The EM algorithm for mixtures of factor analyzers
        Ghahramani, Hinton
        1997, Technical Report CRG-TR-96-1
        Dept. of Computer Science, University of Toronto, Toronto, Canada, MSS 1A4
    3  Max Welling's tutorials 
        (available @ http://www.ics.uci.edu/~welling/classnotes/classnotes.html)
    4  Book of Johnson and Wichern
    5  Book of Joreskog
    6  Book of Duda, Hart
    7  Book of MacKay
    8  Book of Golub, <NAME>
    9  Maximum Likelihood and Covariant Algorithms for Independent Component Analysis
        MacKay
    10 Solving inverse problems using an EM approach to density estimation
        Ghahramani
    11 Maximum likelihood and minimum classification error Factor Analysis for automatic speech recognition
        Saul, Rahim
    12 Finite mixture models
        <NAME>, <NAME> [book of] 
    13 Numerical recipes: the art of scientific computing
        Press et al [book of]
    14 Unsupervised Classification with Non-Gaussian Mixture Models Using ICA
        Lee, Lewicki, Sejnowski
    15 EM Algorithms for PCA and SPCA
        Roweis
    16 A tutorial on Hidden Markov Models and selected applications in speech recognition
        Rabiner
    17 Maximum Likelihood Estimation for Multivariate Mixture Observations of Markov Chains
        Juang, Levinson, Sondhi   

NB Matrix x and y have both shape given by the tuple ('variables nr','samples nr').
"""


from numpy import (array, arange, dot, inner, outer, vdot, cov,
                   diag, ones, eye, zeros, argmax, nanargmax, 
                   mean, std, multiply, sum, product, sqrt,
                   log, abs, exp, power, hstack, vstack, append,
                   concatenate, pi, inf, amin, amax, empty,
                   tanh, any, isnan)
from numpy.linalg import (norm, inv, det, svd, solve, cholesky, linalg)
from numpy.random import (normal, randn, rand, multivariate_normal, uniform)


class lm(object):
    """See details explained above"""

    y, n, p = _default_values_ = None, 0, None
    cumulate = True
    
    def __init__(self, y, *args, **kwargs):

        if not(hasattr(y, '__iter__')):
            raise AttributeError('Input has to an iterable')
        self.y = y
        self.p, self.n = y.shape
             
    def __call__(self, **kw): return self.InferandLearn(**kw)

    def InferandLearn(self, max_iter_nr = 30, **kwargs):
        """
        This method implement and run the EM algorithm, in order to:
        - learn model's parameter A, C, Q, R (learning or system identification)
        - estimate hidden states from observations (inference or filtering or smoothing)

        Inference and Learning, respectively in methods E and M, are overrided by subclasses.
        """
        
        if not isinstance(max_iter_nr, int):
            raise TypeError('The maximum number of iterations of the EM procedure must be an integer')
        if max_iter_nr <= 0:
            raise ValueError('The maximum number of iterations of the EM procedure must be positive')
        for kw, val in kwargs.iteritems(): self.kw = val
        E_EM, M_EM  = self.Inference, self.Learning 
        logLH, Break = self.logLikelihood, self.break_condition
        self.logLH_delta = kwargs.get('logLH_delta', None)
        
        for iter_nr in xrange(max_iter_nr):
            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            # E step
            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            E_EM()
            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            # log-LH
            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            logLH()
            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            # M step
            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            M_EM()
            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            # Break condition of the for loop
            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            if Break(): break
        
        self.iter_nr = iter_nr
        self.trained = True

    def next(self, **kwargs):
        """Calling InferandLearn method to run only one iteration"""

        self.InferandLearn(max_iter_nr = 1, **kwargs)
    
    def Inference(self):
        """The E(xpectation) step of EM algorithm"""
        
        pass

    def Learning(self):
        """The M(aximization) step of EM algorithm"""

        pass
    
    def break_condition(self): 
        """A method verifying an ad-hoc condition to exit from EM iter."""

        pass
    
    def logLikelihood(self):
        """The logLikelihood method for the given model"""

        pass  
    
    def mu_y(self):
        
        return mean(self.y, axis = 1).reshape(self.p, 1)
    
    def cov_obs(self, cov_bias = 1):
        
        return cov(self.y, bias = cov_bias)
    
    def centered_input(self): return self.y - self.mu_y()

    def erase(self): self.y, self.n, self.p = lm._default_values_

    

class fa(lm):
    """
    Factor Analysis model of static data y.
    Model:
        A = 0     (because data are static in time)
        x = x0 = w0,        w0 = N(0, Q)
        y = y0 = Cx + v0,   v0 = N(0, R)
    then
        y ~ N(0, CQC.T + R) 
    and in order to solve any model degeneracy
        Q = I
        R is diagonal
    finally
        y   ~ N(0, C*C.T + R)
        x_y ~ N(beta*y, I-beta*C)   useful for the Inference task
        beta = C.T(C*C.T + R)**-1
    
    C is also called the factor loadings matrix,
    R's diagonal elements are the uniquenesses,
    v the sensor noise.
    
    Hint: apply Template Design Pattern to scale-down from fa
          to spca (or ppca), pca and whitening sub-classes.
    
    Based on ref.1
    """

    k = None
    
    def __init__(self, y, k):
        
        super(fa, self).__init__(y)
        if not isinstance(k, int):
            raise TypeError('k (the number of latent factors) must be an integer')
        if k <= 0:
            raise ValueError('k (the number of latent factors) must be positive')
        if k > self.p:
            raise ValueError('k (the number of latent factors) must not be greater than p (the number of observables)')
        self.k = k
        print self.n, self.k, self.p

        self.initialize()

    def initialize_Q(self):

        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        # Covariance matrix Q of hidden factors = I matrix
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        self.Q = eye(self.k)
    
    def initialize_C(self):
        
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        # Following MDP init settings
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        scale = product(self.yyTdiag) ** (1./self.p)
        self.C = normal(0, sqrt(scale / self.k), size = (self.p, self.k))
    
    def initialize_R(self, with_WN = False):
        """There are other (here commented) ways to init R"""
        
        self.R = self.yyTdiag
        if with_WN: self.R += randn(self.p)
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        # From pag.531 of ref.4 
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        ## self.R = (1-.5*self.p/self.k)*self.yyTdiag

    def initializelogLH(self):
        
        self.logLH = -inf
        self.deltalogLH = inf
        self.logLH_const = -.5 * self.p * log(2. * pi)
        self.logLH_tracks = []
        self.logLH__break = False

    def initialize(self):
        """Initialization step: C and other vars, get observed data covariance"""

        self.arangek, self.arangep = arange(self.k), arange(self.p)
        
        self.yyT = self.cov_obs()
        self.yyTdiag = self.yyT[self.arangep, self.arangep]
        self.initialize_C()
        self.initialize_R()
        self.initialize_Q()
        self.initializelogLH()

    def InferandLearn(self, max_iter_nr = 20, logLH_delta = 1e-3,
                      inferwithLemma = True, **kwargs):

        self.betaInferenceMethod = self.betaInferenceLemma \
                            if inferwithLemma else self.betaInference

        lm.InferandLearn(self, max_iter_nr = max_iter_nr, logLH_delta = logLH_delta)
        
    def break_condition(self):
        
        if -self.logLH_delta < self.deltalogLH < self.logLH_delta:
            self.logLH__break = True
            return True
        return False
    
    def betaInferenceLemma(self):

        C, CT, Rinv = self.C, self.C.T, self.R ** -1
        
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        # Applying the matrix inversion lemma
        # See: - pag.334 in ref.1
        #      - aka Sherman-Morrison-Woodbury formula at pag.50 
        #        in ref.8 (formula (2.1.4)
        #      - aka binomial inverse theorem at 
        #        http://en.wikipedia.org/wiki/Binomialinverse_theorem
        #      - or derived from matrix blockwise inversion as in 
        #        http://en.wikipedia.org/wiki/Invertible_matrix#Blockwiseinversion
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        beta_temp = multiply(CT, Rinv)
        beta = dot(beta_temp, C)
        beta[self.arangek, self.arangek] += 1.
        beta = -dot(C, dot(inv(beta), beta_temp))
        beta[self.arangep, self.arangep] += 1.
        self.logLH_temp = multiply(Rinv.reshape(self.p, 1), beta) ## Rinv * beta #multiply(Rinv, beta)
        self.beta = dot(beta_temp, beta)

    def betaInference(self):

        C, R = self.C, self.R
        CT = C.T
        
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        # Applying the classical method to invert beta
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        beta = dot(C, CT)
        beta[self.arangep, self.arangep] += R
        self.logLH_temp = beta         # = inv(beta)
        self.beta = dot(CT, inv(beta)) #beta)

    def Inference(self):
        """
        E step of EM algorithm
        Inference of sufficient statistic E(x|y) (here x_y)
        
        NB I tried to compute beta via the matrix inversion lemma
           but performances doesn't seem to be better than applying the ordinary formula.
        NB beta.ravel()[arange(0, k**2, k) + arange(k)] += 1.   #<--- code at left is inefficient!
           beta.ravel()[arange(0, p**2, p) + arange(p)] += 1.   #<--- code at left is inefficient!
           self.V.ravel()[arange(0, k**2, k) + arange(k)] += 1. #<--- code at left is inefficient!

        TODO Test if betaInferenceLemma gives performances advantages for (very) high p.
        """
        
        self.betaInferenceMethod()
        self.V = - dot(self.beta, self.C)
        self.V[self.arangek, self.arangek] += 1.

    def Learning_R(self):
        """
        Learning R in M step of EM algorithm
        It is necessary to learn R separately in order to give
        sub-classes able to override the present method
        """
        
        self.R = self.yyTdiag - sum(multiply(self.C, self.delta), axis=1) / self.n

    def Learning(self):
        """
        M step of EM algorithm
        Computing delta and gamma (ref.1 at pag.335).
        Learning and updating model's parameters C and R
        """
        
        delta = dot(self.yyT, self.beta.T)
        self.gamma = self.n * (dot(self.beta, delta) + self.V)
        delta *= self.n
        self.delta = delta
        self.C = dot(self.delta, inv(self.gamma))
        self.Learning_R()

    def logLikelihood(self, mu = None):
        """
        Log-likelihood (LogLH)
        This the LogLH per sample, apart the constant self.logLH_const and
        factor .5 both affecting very sligthly the logLH values convergence
        to a (local, hopefully global) maximum.
        """
        
        logLH_old = self.logLH
        _yyT =  self.yyT - dot(mu, mu.T) if mu != None else self.yyT 
        self.logLH = self.logLH_const - .5 * (-log(det(self.logLH_temp)) + \
                     vdot(_yyT, self.logLH_temp))
        self.deltalogLH = self.logLH - logLH_old
        self.logLH_tracks.append(self.logLH)

    def infer(self):
        
        if self.trained: return self.get_expected_latent(), self.infer_observed()
        print 'Not even trained!'
        return

    def get_expected_latent(self):
        
        if not(self.trained): self.InferandLearn()<commit_msg>Modified method fa.infer<commit_after>try:
    import psyco
    psyco.full()
except ImportError: pass


__doc__ = """
The present module Generative Linear Gaussian Models (GLGM) implements in Python programming language
the following models based mainly on references nr.1 (see References below):

        x(t+1) = Ax(t) + w(t) = Ax(t) + w0,     w0 = N(0, Q)
        y(t)   = Cx(t) + v(t) = Cx(t) + v0,     v0 = N(0, R)
     where:
        x is a (k,1) vector of latent factors or state or causes (hidden) variables
        y is a (p,1) vector of observations
        A is a (k,k) matrix of states transition
        C is a (p,k) matrix of observation measurement or generative matrix
        w is a (k,1) vector of state evolution White Gaussian Noise (WGN)
        v is a (k,1) vector of observation evolution WGN
        w and v are statistically independent of each other and of x and y
        N stands for for Gaussian or Normal probability density function (pdf)

In particular the module aims to implement the following:
    - fa    (Factor Analysis)
    - ppca  (Probabilistic Principal Component Analysis)
    - pca   (Principal Component Analysis)
    - mog   (Mixture of Gaussians)
    - vq    (Vector Quantization)
    - k-means clustering
    - mofa  (Mixture of Factor Analyzers)
    - ica   (Independent Component Analysis)

References:
    1  Unifying Review of Linear Gaussian Models
        Roweis, Ghahramani
        Neural Computation 1999 11:2, 305-345.
    2  The EM algorithm for mixtures of factor analyzers
        Ghahramani, Hinton
        1997, Technical Report CRG-TR-96-1
        Dept. of Computer Science, University of Toronto, Toronto, Canada, MSS 1A4
    3  <NAME>'s tutorials 
        (available @ http://www.ics.uci.edu/~welling/classnotes/classnotes.html)
    4  Book of Johnson and Wichern
    5  Book of Joreskog
    6  Book of <NAME>
    7  Book of MacKay
    8  Book of Golub, Van Loan
    9  Maximum Likelihood and Covariant Algorithms for Independent Component Analysis
        MacKay
    10 Solving inverse problems using an EM approach to density estimation
        Ghahramani
    11 Maximum likelihood and minimum classification error Factor Analysis for automatic speech recognition
        Saul, Rahim
    12 Finite mixture models
        <NAME>, <NAME> [book of] 
    13 Numerical recipes: the art of scientific computing
        Press et al [book of]
    14 Unsupervised Classification with Non-Gaussian Mixture Models Using ICA
        Lee, Lewicki, Sejnowski
    15 EM Algorithms for PCA and SPCA
        Roweis
    16 A tutorial on Hidden Markov Models and selected applications in speech recognition
        Rabiner
    17 Maximum Likelihood Estimation for Multivariate Mixture Observations of Markov Chains
        Juang, Levinson, Sondhi   

NB Matrix x and y have both shape given by the tuple ('variables nr','samples nr').
"""


from numpy import (array, arange, dot, inner, outer, vdot, cov,
                   diag, ones, eye, zeros, argmax, nanargmax, 
                   mean, std, multiply, sum, product, sqrt,
                   log, abs, exp, power, hstack, vstack, append,
                   concatenate, pi, inf, amin, amax, empty,
                   tanh, any, isnan)
from numpy.linalg import (norm, inv, det, svd, solve, cholesky, linalg)
from numpy.random import (normal, randn, rand, multivariate_normal, uniform)


class lm(object):
    """See details explained above"""

    y, n, p = _default_values_ = None, 0, None
    cumulate = True
    
    def __init__(self, y, *args, **kwargs):

        if not(hasattr(y, '__iter__')):
            raise AttributeError('Input has to an iterable')
        self.y = y
        self.p, self.n = y.shape
             
    def __call__(self, **kw): return self.InferandLearn(**kw)

    def InferandLearn(self, max_iter_nr = 30, **kwargs):
        """
        This method implement and run the EM algorithm, in order to:
        - learn model's parameter A, C, Q, R (learning or system identification)
        - estimate hidden states from observations (inference or filtering or smoothing)

        Inference and Learning, respectively in methods E and M, are overrided by subclasses.
        """
        
        if not isinstance(max_iter_nr, int):
            raise TypeError('The maximum number of iterations of the EM procedure must be an integer')
        if max_iter_nr <= 0:
            raise ValueError('The maximum number of iterations of the EM procedure must be positive')
        for kw, val in kwargs.iteritems(): self.kw = val
        E_EM, M_EM  = self.Inference, self.Learning 
        logLH, Break = self.logLikelihood, self.break_condition
        self.logLH_delta = kwargs.get('logLH_delta', None)
        
        for iter_nr in xrange(max_iter_nr):
            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            # E step
            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            E_EM()
            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            # log-LH
            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            logLH()
            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            # M step
            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            M_EM()
            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            # Break condition of the for loop
            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            if Break(): break
        
        self.iter_nr = iter_nr
        self.trained = True

    def next(self, **kwargs):
        """Calling InferandLearn method to run only one iteration"""

        self.InferandLearn(max_iter_nr = 1, **kwargs)
    
    def Inference(self):
        """The E(xpectation) step of EM algorithm"""
        
        pass

    def Learning(self):
        """The M(aximization) step of EM algorithm"""

        pass
    
    def break_condition(self): 
        """A method verifying an ad-hoc condition to exit from EM iter."""

        pass
    
    def logLikelihood(self):
        """The logLikelihood method for the given model"""

        pass  
    
    def mu_y(self):
        
        return mean(self.y, axis = 1).reshape(self.p, 1)
    
    def cov_obs(self, cov_bias = 1):
        
        return cov(self.y, bias = cov_bias)
    
    def centered_input(self): return self.y - self.mu_y()

    def erase(self): self.y, self.n, self.p = lm._default_values_

    

class fa(lm):
    """
    Factor Analysis model of static data y.
    Model:
        A = 0     (because data are static in time)
        x = x0 = w0,        w0 = N(0, Q)
        y = y0 = Cx + v0,   v0 = N(0, R)
    then
        y ~ N(0, CQC.T + R) 
    and in order to solve any model degeneracy
        Q = I
        R is diagonal
    finally
        y   ~ N(0, C*C.T + R)
        x_y ~ N(beta*y, I-beta*C)   useful for the Inference task
        beta = C.T(C*C.T + R)**-1
    
    C is also called the factor loadings matrix,
    R's diagonal elements are the uniquenesses,
    v the sensor noise.
    
    Hint: apply Template Design Pattern to scale-down from fa
          to spca (or ppca), pca and whitening sub-classes.
    
    Based on ref.1
    """

    k = None
    
    def __init__(self, y, k):
        
        super(fa, self).__init__(y)
        if not isinstance(k, int):
            raise TypeError('k (the number of latent factors) must be an integer')
        if k <= 0:
            raise ValueError('k (the number of latent factors) must be positive')
        if k > self.p:
            raise ValueError('k (the number of latent factors) must not be greater than p (the number of observables)')
        self.k = k

        self.initialize()

    def initialize_Q(self):

        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        # Covariance matrix Q of hidden factors = I matrix
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        self.Q = eye(self.k)
    
    def initialize_C(self):
        
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        # Following MDP init settings
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        scale = product(self.yyTdiag) ** (1./self.p)
        self.C = normal(0, sqrt(scale / self.k), size = (self.p, self.k))
    
    def initialize_R(self, with_WN = False):
        """There are other (here commented) ways to init R"""
        
        self.R = self.yyTdiag
        if with_WN: self.R += randn(self.p)
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        # From pag.531 of ref.4 
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        ## self.R = (1-.5*self.p/self.k)*self.yyTdiag

    def initializelogLH(self):
        
        self.logLH = -inf
        self.deltalogLH = inf
        self.logLH_const = -.5 * self.p * log(2. * pi)
        self.logLH_tracks = []
        self.logLH__break = False

    def initialize(self):
        """Initialization step: C and other vars, get observed data covariance"""

        self.arangek, self.arangep = arange(self.k), arange(self.p)
        
        self.yyT = self.cov_obs()
        self.yyTdiag = self.yyT[self.arangep, self.arangep]
        self.initialize_C()
        self.initialize_R()
        self.initialize_Q()
        self.initializelogLH()

    def InferandLearn(self, max_iter_nr = 20, logLH_delta = 1e-3,
                      inferwithLemma = True, **kwargs):

        self.betaInferenceMethod = self.betaInferenceLemma \
                            if inferwithLemma else self.betaInference

        lm.InferandLearn(self, max_iter_nr = max_iter_nr, logLH_delta = logLH_delta)
        
    def break_condition(self):
        
        if -self.logLH_delta < self.deltalogLH < self.logLH_delta:
            self.logLH__break = True
            return True
        return False
    
    def betaInferenceLemma(self):

        C, CT, Rinv = self.C, self.C.T, self.R ** -1
        
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        # Applying the matrix inversion lemma
        # See: - pag.334 in ref.1
        #      - aka Sherman-Morrison-Woodbury formula at pag.50 
        #        in ref.8 (formula (2.1.4)
        #      - aka binomial inverse theorem at 
        #        http://en.wikipedia.org/wiki/Binomialinverse_theorem
        #      - or derived from matrix blockwise inversion as in 
        #        http://en.wikipedia.org/wiki/Invertible_matrix#Blockwiseinversion
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        beta_temp = multiply(CT, Rinv)
        beta = dot(beta_temp, C)
        beta[self.arangek, self.arangek] += 1.
        beta = -dot(C, dot(inv(beta), beta_temp))
        beta[self.arangep, self.arangep] += 1.
        self.logLH_temp = multiply(Rinv.reshape(self.p, 1), beta) ## Rinv * beta #multiply(Rinv, beta)
        self.beta = dot(beta_temp, beta)

    def betaInference(self):

        C, R = self.C, self.R
        CT = C.T
        
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        # Applying the classical method to invert beta
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        beta = dot(C, CT)
        beta[self.arangep, self.arangep] += R
        self.logLH_temp = beta         # = inv(beta)
        self.beta = dot(CT, inv(beta)) #beta)

    def Inference(self):
        """
        E step of EM algorithm
        Inference of sufficient statistic E(x|y) (here x_y)
        
        NB I tried to compute beta via the matrix inversion lemma
           but performances doesn't seem to be better than applying the ordinary formula.
        NB beta.ravel()[arange(0, k**2, k) + arange(k)] += 1.   #<--- code at left is inefficient!
           beta.ravel()[arange(0, p**2, p) + arange(p)] += 1.   #<--- code at left is inefficient!
           self.V.ravel()[arange(0, k**2, k) + arange(k)] += 1. #<--- code at left is inefficient!

        TODO Test if betaInferenceLemma gives performances advantages for (very) high p.
        """
        
        self.betaInferenceMethod()
        self.V = - dot(self.beta, self.C)
        self.V[self.arangek, self.arangek] += 1.

    def Learning_R(self):
        """
        Learning R in M step of EM algorithm
        It is necessary to learn R separately in order to give
        sub-classes able to override the present method
        """
        
        self.R = self.yyTdiag - sum(multiply(self.C, self.delta), axis=1) / self.n

    def Learning(self):
        """
        M step of EM algorithm
        Computing delta and gamma (ref.1 at pag.335).
        Learning and updating model's parameters C and R
        """
        
        delta = dot(self.yyT, self.beta.T)
        self.gamma = self.n * (dot(self.beta, delta) + self.V)
        delta *= self.n
        self.delta = delta
        self.C = dot(self.delta, inv(self.gamma))
        self.Learning_R()

    def logLikelihood(self, mu = None):
        """
        Log-likelihood (LogLH)
        This the LogLH per sample, apart the constant self.logLH_const and
        factor .5 both affecting very sligthly the logLH values convergence
        to a (local, hopefully global) maximum.
        """
        
        logLH_old = self.logLH
        _yyT =  self.yyT - dot(mu, mu.T) if mu != None else self.yyT 
        self.logLH = self.logLH_const - .5 * (-log(det(self.logLH_temp)) + \
                     vdot(_yyT, self.logLH_temp))
        self.deltalogLH = self.logLH - logLH_old
        self.logLH_tracks.append(self.logLH)

    def infer(self):
        
        if self.trained: return self.get_expected_latent(), self.infer_observed()

    def get_expected_latent(self):
        
        if not(self.trained): self.InferandLearn()
{'fraction_non_alphanumeric': 0.16318472252135832, 'fraction_numerical': 0.009045875511522722, 'mean_word_length': 2.2010800873262095, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2661140', 'n_tokens_mistral': 8559, 'n_tokens_neox': 7849, 'n_words': 2563}