python类diagonal()的实例源码

def test_diagonal(self):
        a = np.arange(12).reshape((3, 4))
        assert_equal(a.diagonal(), [0, 5, 10])
        assert_equal(a.diagonal(0), [0, 5, 10])
        assert_equal(a.diagonal(1), [1, 6, 11])
        assert_equal(a.diagonal(-1), [4, 9])

        b = np.arange(8).reshape((2, 2, 2))
        assert_equal(b.diagonal(), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(0), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(1), [[2], [3]])
        assert_equal(b.diagonal(-1), [[4], [5]])
        assert_raises(ValueError, b.diagonal, axis1=0, axis2=0)
        assert_equal(b.diagonal(0, 1, 2), [[0, 3], [4, 7]])
        assert_equal(b.diagonal(0, 0, 1), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(offset=1, axis1=0, axis2=2), [[1], [3]])
        # Order of axis argument doesn't matter:
        assert_equal(b.diagonal(0, 2, 1), [[0, 3], [4, 7]])

def zz(matrix, nb):
    r"""Zig-zag traversal of the input matrix
    :param matrix: input matrix
    :param nb: number of coefficients to keep
    :return: an array of nb coefficients
    """
    flipped = np.fliplr(matrix)
    rows, cols = flipped.shape  # nb of columns

    coefficient_list = []

    for loop, i in enumerate(range(cols - 1, -rows, -1)):
        anti_diagonal = np.diagonal(flipped, i)

        # reversing even diagonals prioritizes the X resolution
        # reversing odd diagonals prioritizes the Y resolution
        # for square matrices, the information content is the same only when nb covers half of the matrix
        #  e.g. [ nb = n*(n+1)/2 ]
        if loop % 2 == 0:
            anti_diagonal = anti_diagonal[::-1]  # reverse anti_diagonal

        coefficient_list.extend([x for x in anti_diagonal])

    # flattened = [val for sublist in coefficient_list for val in sublist]
    return coefficient_list[:nb]

def predict_y(self, inputs):
        """Summary

        Args:
            inputs (TYPE): Description

        Returns:
            TYPE: Description
        """
        mf, vf = self.dyn_layer.forward_prop_thru_post(inputs)
        if self.gp_emi:
            mg, vg = self.emi_layer.forward_prop_thru_post(mf, vf)
            my, vy = self.lik_layer.output_probabilistic(mg, vg)
        else:
            my, _, vy = self.emi_layer.output_probabilistic(mf, vf)
            vy = np.diagonal(vy, axis1=1, axis2=2)
        return my, vy

def get_posterior_y(self):
        """Summary

        Returns:
            TYPE: Description
        """
        mx, vx = self.get_posterior_x()
        if self.Dcon_emi > 0:
            mx = np.hstack((mx, self.x_control))
            vx = np.hstack((vx, np.zeros((self.N, self.Dcon_emi))))
        if self.gp_emi:
            mf, vf = self.emi_layer.forward_prop_thru_post(mx, vx)
            my, vyn = self.lik_layer.output_probabilistic(mf, vf)
        else:
            my, vy, vyn = self.emi_layer.output_probabilistic(mx, vx)
            vf = np.diagonal(vy, axis1=1, axis2=2)
            vyn = np.diagonal(vyn, axis1=1, axis2=2)
        return my, vf, vyn

def test_diagonal(self):
        a = np.arange(12).reshape((3, 4))
        assert_equal(a.diagonal(), [0, 5, 10])
        assert_equal(a.diagonal(0), [0, 5, 10])
        assert_equal(a.diagonal(1), [1, 6, 11])
        assert_equal(a.diagonal(-1), [4, 9])

        b = np.arange(8).reshape((2, 2, 2))
        assert_equal(b.diagonal(), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(0), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(1), [[2], [3]])
        assert_equal(b.diagonal(-1), [[4], [5]])
        assert_raises(ValueError, b.diagonal, axis1=0, axis2=0)
        assert_equal(b.diagonal(0, 1, 2), [[0, 3], [4, 7]])
        assert_equal(b.diagonal(0, 0, 1), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(offset=1, axis1=0, axis2=2), [[1], [3]])
        # Order of axis argument doesn't matter:
        assert_equal(b.diagonal(0, 2, 1), [[0, 3], [4, 7]])

def lsInversion(self):
        """
        LS Inversion from Hwang et al (2002)
        """

        At=np.transpose(self.A)
        St=np.transpose(self.S)
        N=At.dot(self.P).dot(self.A)
        #solution:
        self.X=np.linalg.inv(N+self.S.dot(St)).dot(At).dot(self.P).dot(self.Obs)

        self.r=self.A.dot(self.X)-self.Obs
        rt=np.transpose(self.r)
        self.VtPV=rt.dot(self.P).dot(self.r)
        var_post_norm=self.VtPV/self.dof
        self.SDaposteriori=np.sqrt(var_post_norm)

        cov_post=np.linalg.inv(N)*var_post_norm
        self.var=np.diagonal(cov_post)

def _fisher_vector(self, img_descriptors):
        """
        :param img_descriptors: X
        :return: fisher vector
        :rtype: np.array
        """
        means, covariances, weights = self.gmm.means, self.gmm.covariances, self.gmm.weights
        s0, s1, s2 = self._likelihood_statistics(img_descriptors)
        T = img_descriptors.shape[0]
        diagonal_covariances = np.float32([np.diagonal(covariances[k]) for k in range(0, covariances.shape[0])])
        """ Refer page 4, first column of reference [1] """
        g_weights = self._fisher_vector_weights(s0, s1, s2, means, diagonal_covariances, weights, T)
        g_means = self._fisher_vector_means(s0, s1, s2, means, diagonal_covariances, weights, T)
        g_sigma = self._fisher_vector_sigma(s0, s1, s2, means, diagonal_covariances, weights, T)
        fv = np.concatenate([np.concatenate(g_weights), np.concatenate(g_means), np.concatenate(g_sigma)])
        fv = self.normalize(fv)
        return fv

def compute_mar_likelihood(X_train, X_test, y_train, sigma, l):
    """
    compute log marginal likelihood for tuning parameters using Bayesian optimization
    :param X_train: training data
    :param X_test: test data
    :param y_train: training targets
    :param sigma: output variance
    :param l: lengthscalar
    :return: log marginal likelihood
    """
    s = 0.0005  # noise variance and zero mean for noise
    n = len(X_train)

    # choose RBF kernel in this regression case
    K_train = RBF_kernel(X_train, X_train, sigma, l)
    L = np.linalg.cholesky(K_train + s * np.eye(n))
    m = np.linalg.solve(L, y_train)
    alpha = np.linalg.solve(L.T, m)

    # compute log marginal likelihood
    log_marg_likelihood = -.5 * np.dot(y_train.T, alpha) - np.log(np.diagonal(L)).sum(0) - n / 2.0 * np.log(2 * np.pi)
    return log_marg_likelihood

def test_diagonal(self):
        a = np.arange(12).reshape((3, 4))
        assert_equal(a.diagonal(), [0, 5, 10])
        assert_equal(a.diagonal(0), [0, 5, 10])
        assert_equal(a.diagonal(1), [1, 6, 11])
        assert_equal(a.diagonal(-1), [4, 9])

        b = np.arange(8).reshape((2, 2, 2))
        assert_equal(b.diagonal(), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(0), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(1), [[2], [3]])
        assert_equal(b.diagonal(-1), [[4], [5]])
        assert_raises(ValueError, b.diagonal, axis1=0, axis2=0)
        assert_equal(b.diagonal(0, 1, 2), [[0, 3], [4, 7]])
        assert_equal(b.diagonal(0, 0, 1), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(offset=1, axis1=0, axis2=2), [[1], [3]])
        # Order of axis argument doesn't matter:
        assert_equal(b.diagonal(0, 2, 1), [[0, 3], [4, 7]])

def test_diagonal(self):
        a = np.arange(12).reshape((3, 4))
        assert_equal(a.diagonal(), [0, 5, 10])
        assert_equal(a.diagonal(0), [0, 5, 10])
        assert_equal(a.diagonal(1), [1, 6, 11])
        assert_equal(a.diagonal(-1), [4, 9])

        b = np.arange(8).reshape((2, 2, 2))
        assert_equal(b.diagonal(), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(0), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(1), [[2], [3]])
        assert_equal(b.diagonal(-1), [[4], [5]])
        assert_raises(ValueError, b.diagonal, axis1=0, axis2=0)
        assert_equal(b.diagonal(0, 1, 2), [[0, 3], [4, 7]])
        assert_equal(b.diagonal(0, 0, 1), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(offset=1, axis1=0, axis2=2), [[1], [3]])
        # Order of axis argument doesn't matter:
        assert_equal(b.diagonal(0, 2, 1), [[0, 3], [4, 7]])

def print_performance(cm):
    tp = np.diagonal(cm).astype(np.float)
    tpfp = np.sum(cm, axis=0).astype(np.float) # sum of each col
    tpfn = np.sum(cm, axis=1).astype(np.float) # sum of each row
    acc = np.sum(tp) / np.sum(cm)
    precision = tp / tpfp
    recall = tp / tpfn
    f1 = (2 * precision * recall) / (precision + recall)
    mf1 = np.mean(f1)

    print "Sample: {}".format(np.sum(cm))
    print "W: {}".format(tpfn[W])
    print "N1: {}".format(tpfn[N1])
    print "N2: {}".format(tpfn[N2])
    print "N3: {}".format(tpfn[N3])
    print "REM: {}".format(tpfn[REM])
    print "Confusion matrix:"
    print cm
    print "Precision: {}".format(precision)
    print "Recall: {}".format(recall)
    print "F1: {}".format(f1)
    print "Overall accuracy: {}".format(acc)
    print "Macro-F1 accuracy: {}".format(mf1)

def test_diagonal(self):
        a = np.arange(12).reshape((3, 4))
        assert_equal(a.diagonal(), [0, 5, 10])
        assert_equal(a.diagonal(0), [0, 5, 10])
        assert_equal(a.diagonal(1), [1, 6, 11])
        assert_equal(a.diagonal(-1), [4, 9])

        b = np.arange(8).reshape((2, 2, 2))
        assert_equal(b.diagonal(), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(0), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(1), [[2], [3]])
        assert_equal(b.diagonal(-1), [[4], [5]])
        assert_raises(ValueError, b.diagonal, axis1=0, axis2=0)
        assert_equal(b.diagonal(0, 1, 2), [[0, 3], [4, 7]])
        assert_equal(b.diagonal(0, 0, 1), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(offset=1, axis1=0, axis2=2), [[1], [3]])
        # Order of axis argument doesn't matter:
        assert_equal(b.diagonal(0, 2, 1), [[0, 3], [4, 7]])

def compute_precision_and_recall(confusion_matrix):
        correct_predictions = np.diagonal(confusion_matrix)
        samples_per_class = np.sum(confusion_matrix, axis=0)
        false_positives = np.sum(confusion_matrix, axis=1) - correct_predictions
        false_negatives = samples_per_class - correct_predictions

        prectmp = correct_predictions / (correct_predictions + false_positives)
        prectmp[np.where(correct_predictions == 0)[0]] = 0
        prectmp[np.where(samples_per_class == 0)[0]] = float('nan')
        precision = np.nanmean(prectmp)

        rectmp = correct_predictions / (correct_predictions + false_negatives)
        rectmp[np.where(correct_predictions == 0)[0]] = 0
        rectmp[np.where(samples_per_class == 0)[0]] = float('nan')
        recall = np.nanmean(rectmp)
        return precision, recall

def GetClassMetric(gtLabal, testLabel, numClass=-1, labelSet=npy.array([])):
    if numClass>0:
        labelSet=npy.arange(numClass)
    else:
        if labelSet.size()==0 or npy.min(labelSet)<0:
            return
        numClass=npy.max(labelSet)+1   

    confMat=npy.zeros((numClass,numClass),dtype=npy.float32)
    vecOnes=npy.ones(len(gtLabal))
    for ii in labelSet:
        for jj in labelSet:
            confMat[ii,jj]=npy.sum(vecOnes[npy.logical_and(testLabel==ii, gtLabal==jj)])

    ccn=npy.diagonal(confMat)
    oa=npy.sum(ccn)/npy.sum(confMat) 
    pa=ccn/npy.sum(confMat, axis=0) 
    ua=ccn/npy.sum(confMat, axis=1) 
    temp1=npy.sum(confMat)*npy.sum(ccn)-npy.sum(npy.sum(confMat,axis=1)*npy.sum(confMat,axis=0));
    temp2=npy.power(npy.sum(confMat),2)-npy.sum(npy.sum(confMat,axis=1)*npy.sum(confMat,axis=0));
    kappa=temp1/temp2
    confMat=confMat.astype(npy.int32)
    accMetric={"confMat":confMat, "oa":oa, "pa":pa, "ua":ua, "kappa": kappa}
    return accMetric

def GetClassMetric(gtLabal, testLabel, numClass=-1, labelSet=npy.array([])):
    if numClass>0:
        labelSet=npy.arange(numClass)
    else:
        if labelSet.size()==0 or npy.min(labelSet)<0:
            return
        numClass=npy.max(labelSet)+1   

    confMat=npy.zeros((numClass,numClass),dtype=npy.float32)
    vecOnes=npy.ones(len(gtLabal))
    for ii in labelSet:
        for jj in labelSet:
            confMat[ii,jj]=npy.sum(vecOnes[npy.logical_and(testLabel==ii, gtLabal==jj)])

    ccn=npy.diagonal(confMat)
    oa=npy.sum(ccn)/npy.sum(confMat) 
    pa=ccn/npy.sum(confMat, axis=0) 
    ua=ccn/npy.sum(confMat, axis=1) 
    temp1=npy.sum(confMat)*npy.sum(ccn)-npy.sum(npy.sum(confMat,axis=1)*npy.sum(confMat,axis=0));
    temp2=npy.power(npy.sum(confMat),2)-npy.sum(npy.sum(confMat,axis=1)*npy.sum(confMat,axis=0));
    kappa=temp1/temp2
    confMat=confMat.astype(npy.int32)
    accMetric={"confMat":confMat, "oa":oa, "pa":pa, "ua":ua, "kappa": kappa}
    return accMetric

def write_overlap_densities(
        path_hdf5: str, paths_fragment_overlaps: List, swaps: Matrix, dt: int=1):
    """
    Write the diagonal of the overlap matrices
    """
    logger.info("writing densities in human readable format")

    # Track the crossing between MOs
    for paths_overlaps in paths_fragment_overlaps:
        overlaps = np.stack(retrieve_hdf5_data(path_hdf5, paths_overlaps))
        for k, mtx in enumerate(np.rollaxis(overlaps, 0)):
            overlaps[k] = mtx[:, swaps[k]][swaps[k]]

    # Print to file the densities for each fragment on a given MO
    for ifrag, paths_overlaps in enumerate(paths_fragment_overlaps):
        # time frame
        frames = overlaps.shape[0]
        ts = np.arange(1, frames + 1).reshape(frames, 1) * dt
        # Diagonal of the 3D-tensor
        densities = np.diagonal(overlaps, axis1=1, axis2=2)
        data = np.hstack((ts, densities))
        # Save data in human readable format
        file_name = 'densities_fragment_{}.txt'.format(ifrag)
        np.savetxt(file_name, data, fmt='{:^3}'.format('%e'))

def test_diagonal(self):
        a = np.arange(12).reshape((3, 4))
        assert_equal(a.diagonal(), [0, 5, 10])
        assert_equal(a.diagonal(0), [0, 5, 10])
        assert_equal(a.diagonal(1), [1, 6, 11])
        assert_equal(a.diagonal(-1), [4, 9])

        b = np.arange(8).reshape((2, 2, 2))
        assert_equal(b.diagonal(), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(0), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(1), [[2], [3]])
        assert_equal(b.diagonal(-1), [[4], [5]])
        assert_raises(ValueError, b.diagonal, axis1=0, axis2=0)
        assert_equal(b.diagonal(0, 1, 2), [[0, 3], [4, 7]])
        assert_equal(b.diagonal(0, 0, 1), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(offset=1, axis1=0, axis2=2), [[1], [3]])
        # Order of axis argument doesn't matter:
        assert_equal(b.diagonal(0, 2, 1), [[0, 3], [4, 7]])

def tril(m, k=0):
    """
    Lower triangle of an array.

    Return a copy of an array with elements above the `k`-th diagonal zeroed.

    Parameters
    ----------
    m : array_like, shape (M, N)
        Input array.
    k : int, optional
        Diagonal above which to zero elements.  `k = 0` (the default) is the
        main diagonal, `k < 0` is below it and `k > 0` is above.

    Returns
    -------
    array, shape (M, N)
        Lower triangle of `m`, of same shape and data-type as `m`.

    See Also
    --------
    triu : Same thing, only for the upper triangle.

    """
    return m * tri(m.shape[0], m.shape[1], k=k, dtype=m.dtype)

def test_log_normal_matrix_full():
    n_points, n_components, n_features = 10,5,2

    points = np.random.randn(n_points,n_features)
    means = np.random.randn(n_components,n_features)
    cov = generate_covariance_matrices_full(n_components,n_features)

    # Beginnig of the test
    log_det_cov = np.log(np.linalg.det(cov))
    precisions = np.linalg.inv(cov)
    log_prob = np.empty((n_points,n_components))
    for i in range(n_components):
        diff = points - means[i]
        y = np.dot(diff,np.dot(precisions[i],diff.T))
        log_prob[:,i] = np.diagonal(y)

    expected_log_normal_matrix = -0.5 * (n_features * np.log(2*np.pi) +
                                         log_prob + log_det_cov)

    predected_log_normal_matrix = _log_normal_matrix(points,means,cov,'full')

    assert_almost_equal(expected_log_normal_matrix,predected_log_normal_matrix)

def test_log_normal_matrix_full():
    n_points, n_components, n_features = 10,5,2

    points = np.random.randn(n_points,n_features)
    means = np.random.randn(n_components,n_features)
    cov = generate.generate_covariance_matrices_full(n_components,n_features)
    cov_chol = np.empty((n_components,n_features,n_features))
    for i in range(n_components):
        cov_chol[i] = linalg.cholesky(cov[i],lower=True)

    # Beginnig of the test
    log_det_cov = np.log(np.linalg.det(cov))
    precisions = np.linalg.inv(cov)
    log_prob = np.empty((n_points,n_components))
    for i in range(n_components):
        diff = points - means[i]
        y = np.dot(diff,np.dot(precisions[i],diff.T))
        log_prob[:,i] = np.diagonal(y)

    expected_log_normal_matrix = -0.5 * (n_features * np.log(2*np.pi) +
                                         log_prob + log_det_cov)

    predected_log_normal_matrix = _log_normal_matrix(points,means,cov_chol,'full')

    assert_almost_equal(expected_log_normal_matrix,predected_log_normal_matrix)

def _log_multivariate_normal_density_full(X, means, covars, min_covar=1.e-7):
    """Log probability for full covariance matrices."""
    n_samples, n_dim = X.shape
    nmix = len(means)
    log_prob = np.empty((n_samples, nmix))
    for c, (mu, cv) in enumerate(zip(means, covars)):
        try:
            cv_chol = linalg.cholesky(cv, lower=True)
        except linalg.LinAlgError:
            # The model is most probably stuck in a component with too
            # few observations, we need to reinitialize this components
            try:
                cv_chol = linalg.cholesky(cv + min_covar * np.eye(n_dim),
                                          lower=True)
            except linalg.LinAlgError:
                raise ValueError("'covars' must be symmetric, "
                                 "positive-definite")

        cv_log_det = 2 * np.sum(np.log(np.diagonal(cv_chol)))
        cv_sol = linalg.solve_triangular(cv_chol, (X - mu).T, lower=True).T
        log_prob[:, c] = - .5 * (np.sum(cv_sol ** 2, axis=1) +
                                 n_dim * np.log(2 * np.pi) + cv_log_det)

    return log_prob

def test_diagonal(self):
        a = np.arange(12).reshape((3, 4))
        assert_equal(a.diagonal(), [0, 5, 10])
        assert_equal(a.diagonal(0), [0, 5, 10])
        assert_equal(a.diagonal(1), [1, 6, 11])
        assert_equal(a.diagonal(-1), [4, 9])

        b = np.arange(8).reshape((2, 2, 2))
        assert_equal(b.diagonal(), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(0), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(1), [[2], [3]])
        assert_equal(b.diagonal(-1), [[4], [5]])
        assert_raises(ValueError, b.diagonal, axis1=0, axis2=0)
        assert_equal(b.diagonal(0, 1, 2), [[0, 3], [4, 7]])
        assert_equal(b.diagonal(0, 0, 1), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(offset=1, axis1=0, axis2=2), [[1], [3]])
        # Order of axis argument doesn't matter:
        assert_equal(b.diagonal(0, 2, 1), [[0, 3], [4, 7]])

def test_item_is_self_similar(self):
        sim_matrix, _ = self.out
        diagonal = np.diagonal(sim_matrix)
        self.assertEqual(diagonal.tolist(), [2.0, 1.0, 2.0, 2.0, 1.0])

def diagonal_mpa(entries, sites):
    """Returns an MPA with ``entries`` on the diagonal and zeros otherwise.

    :param numpy.ndarray entries: one-dimensional array
    :returns: :class:`~mpnum.mparray.MPArray` with rank ``len(entries)``.

    """
    assert sites > 0

    if entries.ndim != 1:
        raise NotImplementedError("Currently only supports diagonal MPA with "
                                  "one leg per site.")

    if sites < 2:
        return mp.MPArray.from_array(entries)

    ldim = len(entries)
    leftmost_ltens = np.eye(ldim).reshape((1, ldim, ldim))
    rightmost_ltens = np.diag(entries).reshape((ldim, ldim, 1))
    center_ltens = np.zeros((ldim,) * 3)
    np.fill_diagonal(center_ltens, 1)
    ltens = it.chain((leftmost_ltens,), it.repeat(center_ltens, sites - 2),
                     (rightmost_ltens,))

    return mp.MPArray(LocalTensors(ltens, cform=(sites - 1, sites)))


#########################
#  More physical stuff  #
#########################

def _unitary_haar(dim, randstate=None):
    """Returns a sample from the Haar measure of the unitary group of given
    dimension.

    :param int dim: Dimension
    :param randn: Function to create real N(0,1) distributed random variables.
        It should take the shape of the output as numpy.random.randn does
        (default: numpy.random.randn)
    """
    z = _zrandn((dim, dim), randstate) / np.sqrt(2.0)
    q, r = qr(z)
    d = np.diagonal(r)
    ph = d / np.abs(d)
    return q * ph

def test_diagonal(self):
        a = [[0, 1, 2, 3],
             [4, 5, 6, 7],
             [8, 9, 10, 11]]
        out = np.diagonal(a)
        tgt = [0, 5, 10]

        assert_equal(out, tgt)

def test_diagonal_view_notwriteable(self):
        # this test is only for 1.9, the diagonal view will be
        # writeable in 1.10.
        a = np.eye(3).diagonal()
        assert_(not a.flags.writeable)
        assert_(not a.flags.owndata)

        a = np.diagonal(np.eye(3))
        assert_(not a.flags.writeable)
        assert_(not a.flags.owndata)

        a = np.diag(np.eye(3))
        assert_(not a.flags.writeable)
        assert_(not a.flags.owndata)

def test_diagonal_memleak(self):
        # Regression test for a bug that crept in at one point
        a = np.zeros((100, 100))
        assert_(sys.getrefcount(a) < 50)
        for i in range(100):
            a.diagonal()
        assert_(sys.getrefcount(a) < 50)

def _half_log_det(self, M):
        """ Return log(|M|)*0.5. For positive definite matrix M
            of more than 2 dimensions, calculate this for the
            last two dimension and return a value corresponding
            to each element in the first few dimensions.
        """
        chol = np.linalg.cholesky(M)
        if M.ndim == 2:
            return np.sum(np.log(np.abs(np.diag(chol))))
        else:
            return np.sum(np.log(np.abs(np.diagonal(
                chol, axis1=-2, axis2=-1))), axis=-1)

def predict_forward_mm(self, T, x_control):
        """Summary

        Args:
            T (TYPE): Description
            x_control (None, optional): Description

        Returns:
            TYPE: Description
        """
        mx = np.zeros((T, self.Din))
        vx = np.zeros((T, self.Din))
        my = np.zeros((T, self.Dout))
        vy_noiseless = np.zeros((T, self.Dout))
        vy = np.zeros((T, self.Dout))
        post_m, post_v = self.get_posterior_x()
        mtm1 = post_m[[-1], :]
        vtm1 = post_v[[-1], :]
        for t in range(T):
            if self.Dcon_dyn > 0:
                mtm1 = np.hstack((mtm1, x_control[[t], :]))
                vtm1 = np.hstack((vtm1, np.zeros((1, self.Dcon_dyn))))
            mt, vt = self.dyn_layer.forward_prop_thru_post(mtm1, vtm1)
            if self.Dcon_emi > 0:
                mtc = np.hstack((mt, x_control[[t], :]))
                vtc = np.hstack((vt, np.zeros((1, self.Dcon_emi))))
            else:
                mtc, vtc = mt, vt
            if self.gp_emi:
                mft, vft = self.emi_layer.forward_prop_thru_post(mtc, vtc)
                myt, vyt_n = self.lik_layer.output_probabilistic(mft, vft)
            else:
                myt, vyt, vyt_n = self.emi_layer.output_probabilistic(mt, vt)
                vft = np.diagonal(vyt, axis1=1, axis2=2)
                vyt_n = np.diagonal(vyt_n, axis1=1, axis2=2)
            mx[t, :], vx[t, :] = mt, vt
            my[t, :], vy_noiseless[t, :], vy[t, :] = myt, vft, vyt_n
            mtm1 = mt
            vtm1 = vt
        return mx, vx, my, vy_noiseless, vy