def rand(self):

        m, n = self.__m, self.__n

        s = linalg.cholesky(self.__prod).transpose()
        w = self.__weight

        # Compute the parameters of the posterior distribution.
        mu = linalg.solve(s[:m, :m], s[:m, m:])
        omega = np.dot(s[:m, :m].transpose(), s[:m, :m])
        sigma = np.dot(s[m:, m:].transpose(), s[m:, m:]) / w
        eta = w

        # Simulate the marginal Wishart distribution.
        f = linalg.solve(np.diag(np.sqrt(2.0*random.gamma(
            (eta - np.arange(n))/2.0))) + np.tril(random.randn(n, n), -1),
                         np.sqrt(eta)*linalg.cholesky(sigma).transpose())
        b = np.dot(f.transpose(), f)

        # Simulate the conditional Gauss distribution.
        a = mu + linalg.solve(linalg.cholesky(omega).transpose(),
                              np.dot(random.randn(m, n),
                                     linalg.cholesky(b).transpose()))

        return a, b

def compute_trajectory(which_compute = False):
    
    b = 10E-4
    g = 9.8
    t = .1 #this is the delta t
    sx,sy,vx,vy = 0,0,300,600
    Q = .1*np.eye(4)
    R = 500*np.eye(2)
    F = np.array(([1.,0.,t,0.],[0.,1.,0.,t],[0.,0.,1-b,0.],[0.,0.,0.,1-b]))
    u = np.array([0,0,0,-g*t])
    x0 = np.array([sx,sy,vx,vy])
    rand_vec = np.random.rand(4)
    w = np.dot((la.cholesky(Q)).T, rand_vec)
    N = 1200
    x_list = []

    if which_compute == True:
        y = []
        H = np.eye(4)[:2]
        R = np.eye(2)*500

    for i in xrange(N):
        if which_compute == True and i+1 >= 400 and i+1 <= 600:
            y.append(np.dot(H,x0) + np.dot(la.cholesky(R),np.random.randn(2)))
        x_new = np.dot(F,x0) + u + np.dot(la.cholesky(Q),np.random.randn(4))
        x0 = x_new
        x_list.append(x_new)


    if which_compute == True:
        return np.array(x_list), np.array(y) #.shape
    else:
        return np.array(x_list)

def decomposeSubsetKM(K_r, K_rr):
    """decomposes r*m kernel matrix, where r is the number of basis vectors and m the
    number of training examples
    
    @param K_r: r*m kernel matrix, where only the lines corresponding to basis vectors are present
    @type K_r: numpy matrix
    @param basis_vectors: the indices of the basis vectors
    @type basis_vectors: list of integers
    @return svals, evecs, U, C_T_inv
    @rtype tuple of numpy matrices"""
    #K_rr = K_r[:, basis_vectors]
    #C = la.cholesky(K_rr)
    try:
        C = la.cholesky(K_rr)
    except LinAlgError:
        print "Warning: chosen basis vectors not linearly independent"
        print "Shifting the diagonal of kernel matrix"
        #__shiftKmatrix(K_r, basis_vectors)
        #K_rr = K_r[:, basis_vectors]
        #C = la.cholesky(K_rr)
        C = la.cholesky(K_rr+0.000000001 * np.eye(K_rr.shape[0]))
    C_T_inv = la.inv(C.T)
    H = np.dot(K_r.T, C_T_inv)
    svals, evecs, U = decomposeDataMatrix(H.T)
    return svals, evecs, U, C_T_inv

def CCA(X, Y, eps=1.e-15):
    """
    Canonical corelation analysis of two matrices
    
    Parameters
    ----------
    X array of shape (nbitem,p) 
    Y array of shape (nbitem,q) 
    eps=1.e-15, float is a small biasing constant
                to grant invertibility of the matrices
    
    Returns
    -------
    ccs, array of shape(min(n,p,q) the canonical correlations
        
    Note
    ----
    It is expected that nbitem>>max(p,q)
    """
    from numpy.linalg import cholesky, inv, svd
    if Y.shape[0]!=X.shape[0]:
        raise ValueError,"Incompatible dimensions for X and Y"
    p = X.shape[1]
    q = Y.shape[1]
    sqX = np.dot(X.T,X)
    sqY = np.dot(Y.T,Y)
    sqX += np.trace(sqX)*eps*np.eye(p)
    sqY += np.trace(sqY)*eps*np.eye(q)
    rsqX = cholesky(sqX)
    rsqY = cholesky(sqY)
    iX = inv(rsqX).T
    iY = inv(rsqY).T
    Cxy = np.dot(np.dot(X,iX).T,np.dot(Y,iY))
    uv, ccs, vv = svd(Cxy)
    return ccs

def gp_pred(logtheta, covfunc, X, y, Xstar, R=None, w=None, Rstar=None):
    
    #else:
    #    print '        xgp_pred()'
        
    # compute training set covariance matrix (K) and
    # (marginal) test predictions (Kss = self-cov; Kstar = corss-cov)
    if R==None:
        K = feval(covfunc, logtheta, X)                     # training covariances
        [Kss, Kstar] = feval(covfunc, logtheta, X, Xstar)   # test covariances (Kss = self covariances, Kstar = cov between train and test cases)
    else:
        K = feval(covfunc, logtheta, X, R, w)               # training covariances
        [Kss, Kstar] = feval(covfunc, logtheta, X, R, w, Xstar, Rstar)   # test covariances
   # K += sn2*eye(X.shape[0]
    try:
        n = X.shape[0]
        K = K + identity(n)*sn2 #numerical stability shit
    except TypeError:
        raise Exception(str(K) + " " + str(X.shape) + " " + str(identity(sn2).shape) + " " + str(np.array(K).shape))
    try:
        L = linalg.cholesky(K,lower=True) # lower triangular matrix
    except linalg.LinAlgError:
        L = linalg.cholesky(nearPD(K),lower=True)
    #L = linalg.cholesky(K,lower=True) # lower triangular matrix    
    alpha = solve_chol(L.transpose(),y)         # compute inv(K)*y
      
    out1 = dot(Kstar.transpose(),alpha)         # predicted means
    v = linalg.solve(L, Kstar)                  
    tmp=v*v                
    out2 = Kss - array([tmp.sum(axis=0)]).transpose()  # predicted variances  

    return [out1, out2]

 def __init__(self, mu=array([0, 0]), Sigma=eye(2), is_cholesky=False, ell=None):
     Distribution.__init__(self, len(Sigma))
     
     assert(len(shape(mu)) == 1)
     assert(max(shape(Sigma)) == len(mu))
     self.mu = mu
     self.ell = ell
     if is_cholesky: 
         self.L = Sigma
         if ell == None:
             assert(shape(Sigma)[0] == shape(Sigma)[1])
         else:
             assert(shape(Sigma)[1] == ell)
     else: 
         assert(shape(Sigma)[0] == shape(Sigma)[1])
         if ell is not None:
             self.L, _, _ = MatrixTools.low_rank_approx(Sigma, ell)
             self.L = self.L.T
             assert(shape(self.L)[1] == ell)
         else:
             try:
                 self.L = cholesky(Sigma)
             except LinAlgError:
                 # some really crude check for PSD (which only corrects for orunding errors
                 self.L = cholesky(Sigma+eye(len(Sigma))*1e-5)

    def rand(self):

        dim=self.__dim__

        mu=self.__param__.mu
        omega=self.__param__.omega
        sigma=self.__param__.sigma
        eta=self.__param__.eta

        if numpy.isfinite(eta):

            # Simulate the marginal Wishart distribution.
            diag=2.0*random.gamma((eta-numpy.arange(dim))/2.0)
            fact=numpy.diag(numpy.sqrt(diag))+numpy.tril(random.randn(dim,dim),-1)
            fact=linalg.solve(fact,math.sqrt(eta)*linalg.cholesky(sigma).transpose())
            disp=numpy.dot(fact.transpose(),fact)

        else:

            # Account for the special case where the
            # marginal distribution is singular.
            disp=numpy.copy(sigma)

        if numpy.isfinite(omega):

            # Simulate the conditional Gauss distribution.
            loc=mu+numpy.dot(linalg.cholesky(disp),random.randn(dim))/math.sqrt(omega)

        else:

            # Account for the special case where the
            # conditional distribution is singular.
            loc=numpy.copy(mu)

        return loc,disp

    def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084, -2.1825367], [-2.1825367, 9.85910595]], dtype=">f8")
        b = array(a, dtype="<f8")

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp)

 def is_positive(matrix):
     """ All main main minors are positive
     """
     try:
         cholesky(matrix)
         return True
     except LinAlgError:
         return False

    def factor(self, X, rho):
        n, p = X.shape
        if n >= p:
            L = li.cholesky(np.dot(X.T, X) + rho * np.eye(p))
        else:
            L = li.cholesky(np.eye(n) + 1.0 / rho * np.dot(X, X.T))

        return L, L.T  # L, U

def _is_sympd(M):
    'Check that the matrix M is symmetric positive definite'
    try:
        linalg.cholesky(M)  # check that is a symmetric pd
        is_sympd = True 
    except linalg.LinAlgError as err:
        is_sympd = False
    return is_sympd

    def __init__(self,initmean,initvar,transgain,transnoise,measgain,measnoise):

        # Check the initial mean.
        try:
            numstate,=numpy.shape(initmean)
        except ValueError:
            raise Exception('Initial mean must be a vector.')

        # Check the initial variance.
        if numpy.shape(initvar)!=(numstate,numstate):
            raise Exception('Initial variance must be a {}-by-{} matrix.'.format(numstate,numstate))
        if not numpy.allclose(numpy.transpose(initvar),initvar):
            raise Exception('Initial variance matrix must be symmetric.')
        try:
            cholfact=linalg.cholesky(initvar)
        except linalg.LinAlgError:
            raise Exception('Initial variance matrix must be positive-definite.')

        # Check the transition gain.
        if numpy.ndim(transgain)!=2 or numpy.shape(transgain)!=(numstate,numstate):
            raise Exception('Transition gain must be a {}-by-{} matrix.'.format(numstate,numstate))

        # Check the transition noise.
        if numpy.ndim(transnoise)!=2 or numpy.shape(transnoise)!=(numstate,numstate):
            raise Exception('Transition noise must be a {}-by-{} matrix.'.format(numstate,numstate))
        if not numpy.allclose(numpy.transpose(transnoise),transnoise):
            raise Exception('Transition noise matrix must be symmetric.')
        if numpy.any(linalg.eigvalsh(transnoise)<0.0):
            raise Exception('Transition noise matrix must be positive-semi-definite.')

        # Check the measurement gain.
        try:
            numobs,numcol=numpy.shape(measgain)
        except ValueError:
            raise Exception('Measurement gain must be a matrix.')
        if numcol!=numstate:
            raise Exception('Measurement gain matrix must have {} columns.'.format(numstate))

        # Check the measurement noise.
        if numpy.ndim(measnoise)!=2 or numpy.shape(measnoise)!=(numobs,numobs):
            raise Exception('Measurement noise must be a {}-by-{} matrix.'.format(numobs,numobs))
        if not numpy.allclose(numpy.transpose(measnoise),measnoise):
            raise Exception('Measurement noise matrix must be symmetric.')
        try:
            cholfact=linalg.cholesky(measnoise)
        except linalg.LinAlgError:
            raise Exception('Measurement noise matrix must be positive-definite.')

        # Set the model.
        self.initmean=numpy.asarray(initmean)
        self.initvar=numpy.asarray(initvar)
        self.transgain=numpy.asarray(transgain)
        self.transnoise=numpy.asarray(transnoise)
        self.measgain=numpy.asarray(measgain)
        self.measnoise=numpy.asarray(measnoise)

        self.__size__=numstate,numobs

def factor(X,rho):
    m,n = X.shape
    if m>=n:
       L = cholesky(X.T.dot(X)+rho*sparse.eye(n))
    else:
       L = cholesky(sparse.eye(m)+1./rho*(X.dot(X.T)))
    L = sparse.csc_matrix(L)
    U = sparse.csc_matrix(L.T)
    return L,U

def induceRankCorr(R, Cstar):
    """Induces rank correlation Cstar onto a sample R [N x k].
    Note that it is easy to specify correlations that are not possible to generate.
    Results generated with a given Cstar should be checked.

    Iman, R. L., and W. J. Conover. 1982. A Distribution-free Approach to Inducing Rank
    Correlation Among Input Variables. Communications in Statistics: Simulation and
    Computations 11:311-334.
    
    Parameters
    ----------
    R : ndarray [N x k]
        Matrix of random samples (with no pre-existing correlation)
    Cstar : ndarray [k x k]
        Desired positive, symetric correlation matrix with ones along the diagonal.
    
    Returns
    -------
    corrR : ndarray [N x k]
        A correlated matrix of samples."""

    """Define inverse complimentary error function (erfcinv in matlab)
    x is on interval [0,2]
    its also defined in scipy.special"""
    #erfcinv = lambda x: -stats.norm.ppf(x/2)/sqrt(2)

    C = Cstar
    N, k = R.shape
    """Calculate the sample correlation matrix T"""
    T = np.corrcoef(R.T)

    """Calculate lower triangular cholesky
        decomposition of Cstar (i.e. P*P' = C)"""
    P = cholesky(C).T

    """Calculate lower triangular cholesky decomposition of T, i.e. Q*Q' = T"""
    Q = cholesky(T).T

    """S*T*S' = C"""
    S = P.dot(inv(Q))

    """Replace values in samples with corresponding
    rank-indices and convert to van der Waerden scores"""

    RvdW = -np.sqrt(2) * special.erfcinv(2*((_columnRanks(R)+1)/(N+1)))

    """Matrix RBstar has a correlation matrix exactly equal to C"""
    RBstar = RvdW.dot(S.T)
    
    """Match up the rank pairing in R according to RBstar"""
    ranks = _columnRanks(RBstar)
    sortedR = np.sort(R, axis=0)
    corrR = np.zeros(R.shape)
    for j in np.arange(k):
        corrR[:, j] = sortedR[ranks[:, j], j]

    return corrR

def factor(A, rho):
    m, n = A.shape
    if m >= n:  # if skinny
        AA = np.dot(A.T, A) + rho * np.identity(n)
        L = npl.cholesky(AA)
    else:  # if fat
        AA = np.indentity(m) + 1.0 / rho * np.dot(A, A.T)
        L = npl.cholesky(AA)
    U = L.T
    return L, U

def factor(A, rho):
    m, n = A.shape;
    At = np.transpose(A)

    if m >= n:    # if skinny
        L = npl.cholesky(np.dot(At, A) + rho * np.eye(n));
    else:         # if fat
        L = npl.cholesky(np.eye(m) + 1 / rho * np.dot(A, At));

    U = np.transpose(L)

    return (L, U)

	def test_mogsm(self):
		mcgsm = MCGSM(
			dim_in=0,
			dim_out=3,
			num_components=2,
			num_scales=2,
			num_features=0)

		p0 = 0.3
		p1 = 0.7
		N = 20000
		m0 = array([[2], [0], [0]])
		m1 = array([[0], [2], [1]])
		C0 = cov(randn(mcgsm.dim_out, mcgsm.dim_out**2))
		C1 = cov(randn(mcgsm.dim_out, mcgsm.dim_out**2))
		input = zeros([0, N])
		output = hstack([
			dot(cholesky(C0), randn(mcgsm.dim_out, round(p0 * N))) + m0,
			dot(cholesky(C1), randn(mcgsm.dim_out, round(p1 * N))) + m1]) * (rand(1, N) + 0.5)

		mcgsm.train(input, output, parameters={
			'verbosity': 0,
			'max_iter': 10,
			'train_means': True})

		mogsm = MoGSM(3, 2, 2)

		# translate parameters from MCGSM to MoGSM
		mogsm.priors = sum(exp(mcgsm.priors), 1) / sum(exp(mcgsm.priors))

		for k in range(mogsm.num_components):
			mogsm[k].mean = mcgsm.means[:, k]
			mogsm[k].covariance = inv(dot(mcgsm.cholesky_factors[k], mcgsm.cholesky_factors[k].T))
			mogsm[k].scales = exp(mcgsm.scales[k, :])
			mogsm[k].priors = exp(mcgsm.priors[k, :]) / sum(exp(mcgsm.priors[k, :]))

		self.assertAlmostEqual(mcgsm.evaluate(input, output), mogsm.evaluate(output), 5)

		mogsm_samples = mogsm.sample(N)
		mcgsm_samples = mcgsm.sample(input)

		# generated samples should have the same distribution
		for i in range(mogsm.dim):
			self.assertTrue(ks_2samp(mogsm_samples[i], mcgsm_samples[0]) > 0.0001)
			self.assertTrue(ks_2samp(mogsm_samples[i], mcgsm_samples[1]) > 0.0001)
			self.assertTrue(ks_2samp(mogsm_samples[i], mcgsm_samples[2]) > 0.0001)

		posterior = mcgsm.posterior(input, mcgsm_samples)

		# average posterior should correspond to prior
		for k in range(mogsm.num_components):
			self.assertLess(abs(1 - mean(posterior[k]) / mogsm.priors[k]), 0.1)

    def compute_cholesky(self, x, train=False):
        """
        Computes the cholesky decomposition

        """
        Kxx = self.kernel(x, x)
        # import pylab as pl
        # pl.imshow(Kxx)
        # pl.show()
        if not train:
            self.cholesky = cholesky(Kxx)
        else:
            self.train_cholesky = cholesky(Kxx)

 def _is_pos_def(self, A):
     """
     Check if matrix A is positive definite. Code from
     https://stackoverflow.com/questions/16266720/find-out-if-matrix-is-positive-definite-with-numpy
     """
     if allclose(A, A.T, rtol=1e-10, atol=1e-12):
         try:
             cholesky(A)
             return True
         except LinAlgError:
             return False
     else:
         return False

def factor(A, rho):
    """ Returns the factorisation of rho*I + At*A"""
    m, n = A.shape;
    At = np.transpose(A)

    if m >= n:    # if skinny
        L = npl.cholesky(np.dot(At, A) + rho * np.eye(n));
    else:         # if fat
        L = npl.cholesky(np.eye(m) + (1 / float(rho)) * np.dot(A, At));

    U = np.transpose(L)

    return (L, U)