def LDAnalysis_manual(X, y):
    n_features = X.shape[1]
    n_classes = len(np.unique(y))

    print("Mean vectors...")
    mean_vectors = []
    for cl in range(n_classes):
        mean_vectors.append(np.mean(X[y == cl], axis=0))
        # print("Mean vector class {}: {}".format(cl, mean_vectors[cl - 1]))

    print("In-class scatter matrix...")
    S_W = np.zeros((n_features, n_features))
    for cl, mv in zip(range(1, n_classes), mean_vectors):
        class_sc_mat = np.zeros((n_features, n_features))  # each class' scatter matrix
        for row in X[y == cl]:
            row, mv = row.reshape(n_features, 1), mv.reshape(n_features, 1)  # column vectors
            class_sc_mat += (row - mv).dot((row - mv).T)
        S_W += class_sc_mat  # sum class scatter matrices

    overall_mean = np.mean(X, axis=0)
    print("Between-class scatter matrix...")
    S_B = np.zeros((n_features, n_features))

    for i, mean_vec in enumerate(mean_vectors):
        n = X[y == i + 1].shape[0]
        mean_vec = mean_vec.reshape(n_features, 1)  # make column vector
        overall_mean = overall_mean.reshape(n_features, 1)
        S_B += n * (mean_vec - overall_mean).dot((mean_vec - overall_mean).T)

    eig_vals, eig_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_B))

    print("Eigenvector test")
    for i in range(len(eig_vals)):
        print("r{:3}".format(i), end=" ")
        sys.stdout.flush()

        eigv = eig_vecs[:, i].reshape(n_features, 1)
        np.testing.assert_array_almost_equal(np.linalg.inv(S_W).dot(S_B).dot(eigv).real,
                                             (eig_vals[i] * eigv).real,
                                             decimal=6, err_msg='', verbose=True)
    __log.debug("nAll values ok.")

    eig_pairs = [(np.abs(eig_vals[i]), eig_vecs[:, i]) for i in
                 range(len(eig_vals))]  # make list of value & vector tuples
    eig_pairs = sorted(eig_pairs, key=lambda k: k[0], reverse=True)  # Sort tuple-list from high to low

    __log.info("nEigenvalues (decending):")
    for i in eig_pairs:
        __log.info(i[0])

    tot = sum(eig_vals)
    var_exp = [(i / tot) for i in sorted(eig_vals, reverse=True)]
    cum_var_exp = np.cumsum(var_exp)
    cum_var_exp = cum_var_exp.real
    plot(len(var_exp), var_exp, cum_var_exp)

    idx_98 = next(idx for idx, val in enumerate(cum_var_exp) if val > .98)
    return idx_98 + 1