michalk8 commented 4 years ago

I've ported most of the tests from the matlab repo, some don't pass because the checks are not implemented in some of the functions specifically (but are present in the functions which call them). I've also tried to preserve the test names and variables as in the Matlab repo.

TODOs:

[ ] whether sparse/dense implemetation give similar results
[ ] test brandts/krylov/scipy
[ ] test missing high lvl function (like gpcca_coarse_grain).

closes #17 closes #4

michalk8 commented 4 years ago

The only test that fails and is worrying me is this one: https://github.com/msmdev/msmtools/pull/18/files#diff-a77e704d8e9642584dac7cea5a986e5cf755403a1056410c6718fc75b012d381R60 I've checked whether I've made a typo when copying the data, but haven't found anything there.

Marius1311 commented 4 years ago

to do's for Marius

add an example with absorbing states
add assertion for invariant membership vectors
add new test with transition matrix from Supplementary Fig.

Marius1311 commented 4 years ago

Thanks for adding these new test @michalk8 !

Marius1311 commented 4 years ago

Hi Mike, this notebook in our reproduciblity repo contains the following very small toy example matrix you could use:

:
p = np.array([
    # 0.   1.   2.   3.   4.   5.   6.   7.   8.   9.   10.  11.   
    [0.0, 0.8, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], #0
    [0.2, 0.0, 0.6, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], #1
    [0.6, 0.2, 0.0, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], #2

    [0.0, 0.05, 0.05, 0.0, 0.6, 0.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], #3
    [0.0, 0.0, 0.0, 0.25, 0.0, 0.25, 0.4, 0.0, 0.0, 0.1, 0.0, 0.0], #4
    [0.0, 0.0, 0.0, 0.25, 0.25, 0.0, 0.1, 0.0, 0.0, 0.4, 0.0, 0.0], #5

    [0.0, 0.0, 0.0, 0.0, 0.05, 0.05, 0.0, 0.7, 0.2, 0.0, 0.0, 0.0], #6
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.2, 0.0, 0.8, 0.0, 0.0, 0.0], #7
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.2, 0.0, 0.0, 0.0, 0.0], #8

    [0.0, 0.0, 0.0, 0.0, 0.05, 0.05, 0.0, 0.0, 0.0, 0.0, 0.7, 0.2], #9
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.2, 0.0, 0.8], #10
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.2, 0.0], #11
            ])

I wanted to do this myself but things are a bit hectic on my end atm, so I though I will just give you the matrix for now in case you find the time to already set up some tests.

Marius1311 commented 4 years ago

Besides some tests on this irreducible example, you could introduce an absorbing state (e.g. make state 11 absorbing by giving it 0 outgoing probability and 1 self-loop prob) and apply GPCCA. As far as I'm aware, this fails atm, something with D-orthogonality, and we need to look into this.

Marius1311 commented 3 years ago

Hi Mike, any updates here?

Marius1311 commented 3 years ago

@michalk8, I recall we had one failing tests where membership vectors differed between brands and krylov or something similar. I went through old PRs and realised that I discussed this with Bernhard before, see here: https://github.com/msmdev/msmtools/pull/8

Indeed, the Schur decomposition is not unique, so we would not expect to get exactly the same schur vectors back. However, we expect them to span the same invariant subspaces, which you can check by using the subspace_angles, method, see the utility function I posted in the linked PR.

BTW, this is a nice way to check whether something is a correct schur form: np.max(subspace_angles(T.A @ X_b, X_b @ sf_b)), see my utility function.

So, this should fix another failing test. Could you please describe in bullet points what's still failing and why it is? Then I can give you more targeted feedback.

Marius1311 commented 3 years ago

Hi Mike, would be great to have some test cases with the toy matrix I posted above ready for our chat tomorrow. In general, we should have everything ready by Friday to discuss with Bernhard.

michalk8 commented 3 years ago

I've added 1 test, however 1 check there fails: https://github.com/msmdev/msmtools/pull/18/commits/6a4e9e5fca957d086bc981e8d5a0136a9bd54f3b#diff-a77e704d8e9642584dac7cea5a986e5cf755403a1056410c6718fc75b012d381R623

I've regenerated the ground truth for sparse matrices, so the only test failing apart from the above mentioned is: https://github.com/msmdev/msmtools/pull/18/commits/6a4e9e5fca957d086bc981e8d5a0136a9bd54f3b#diff-a77e704d8e9642584dac7cea5a986e5cf755403a1056410c6718fc75b012d381R98

Lastly, any idea why this happens when mu=0: https://github.com/msmdev/msmtools/pull/18/commits/6a4e9e5fca957d086bc981e8d5a0136a9bd54f3b#diff-a77e704d8e9642584dac7cea5a986e5cf755403a1056410c6718fc75b012d381R595

michalk8 commented 3 years ago

There are all the tests that fail for me (2 from GPCCA, the rest from other modules):

FAILED tests/analysis/test_mfpt.py::TestMfptSparse::test_mfpt_between_sets - TypeError: Cannot use scipy.linalg.eig for sparse A with k >= N - 1. Use scipy.linalg.eig(A.toarray()) or reduce k.
FAILED tests/analysis/impl/dense/gpcca_test.py::TestGPCCAMatlabRegression::test_normal_case - AssertionError: 
FAILED tests/analysis/impl/dense/gpcca_test.py::TestCustom::test_P2 - AssertionError: array([[7.77156117e-16, 1.27631890e-03],
FAILED tests/analysis/impl/sparse/mean_first_passage_time_test.py::TestMfpt::test_mfpt_between_sets - TypeError: Cannot use scipy.linalg.eig for sparse A with k >= N - 1. Use scipy.linalg.eig(A.toarray()) or reduce k.
FAILED tests/estimation/tests/test_effective_count_matrix.py::TestEffectiveCountMatrix::test_njobs_speedup - AssertionError: 1.3874175548553467 not less than 1.3511504530906677 : does not scale for njobs=4
FAILED tests/estimation/tests/test_transition_matrix.py::TestTransitionRevPiSym::test_sparse - TypeError: Cannot use scipy.linalg.eig for sparse A with k >= N - 1. Use scipy.linalg.eig(A.toarray()) or reduce k.

I've tested it with scipy version: 1.5.3.

Marius1311 commented 3 years ago

I've added 1 test, however 1 check there fails: 6a4e9e5#diff-a77e704d8e9642584dac7cea5a986e5cf755403a1056410c6718fc75b012d381R623

For this one, we don't know what's going on. Need to discuss

I've regenerated the ground truth for sparse matrices, so the only test failing apart from the above mentioned is: 6a4e9e5#diff-a77e704d8e9642584dac7cea5a986e5cf755403a1056410c6718fc75b012d381R98

This one looks fine, was a problem with membership vectors being not exactly the same but they do actually span the same subspace.

Lastly, any idea why this happens when mu=0: 6a4e9e5#diff-a77e704d8e9642584dac7cea5a986e5cf755403a1056410c6718fc75b012d381R595

Need to discuss.

Marius1311 commented 3 years ago

Reducible matrix

I looked into applying GPCCA to reducible matrices. First, I tried a case which had one transient and one single recurrent class: all good here, GPCCA correctly recovers the macrostates. Then, I moved onto an example with one transient and two recurrent classes. This is where the trouble starts.

The matrix

This is the matrix:

p = np.array([
    # 0.   1.   2.   3.   4.   5.   6.   7.   8.   9.   10.  11.   
    [0.0, 0.8, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], #0
    [0.2, 0.0, 0.6, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], #1
    [0.6, 0.2, 0.0, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], #2

    [0.0, 0.05, 0.05, 0.0, 0.6, 0.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], #3
    [0.0, 0.0, 0.0, 0.25, 0.0, 0.25, 0.4, 0.0, 0.0, 0.1, 0.0, 0.0], #4
    [0.0, 0.0, 0.0, 0.25, 0.25, 0.0, 0.1, 0.0, 0.0, 0.4, 0.0, 0.0], #5

    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.7, 0.3, 0.0, 0.0, 0.0], #6
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.2, 0.0, 0.8, 0.0, 0.0, 0.0], #7
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.2, 0.0, 0.0, 0.0, 0.0], #8

    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.7, 0.3], #9
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.2, 0.0, 0.8], #10
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.2, 0.0], #11
            ])

Maybe a bit easier to see when I plot the sate graph:

Essentially, (6, 7, 8) form one recurrent class while (9, 10, 11) form the other.

The problem

When I run g.optimize, I get

~/Projects/msmtools/msmtools/analysis/dense/gpcca.py in _do_schur(P, eta, m, z, method, tol_krylov)
    256     if not np.allclose(X.T.dot(X*eta[:, None]), np.eye(X.shape[1]), atol=1e-6, rtol=1e-5):
    257         print(X.T.dot(X*eta[:, None]))
--> 258         raise ValueError("Schur vectors appear to not be D-orthogonal.")
    259 
    260     # Raise, if X doesn't fullfill the invariant subspace condition!

ValueError: Schur vectors appear to not be D-orthogonal.

It also prints the matrix X^T D X, which is

[[ 1.00000000e+00 -5.64545831e-01  2.22044605e-16 -1.10328413e-15]
 [-5.64545831e-01  1.00000000e+00  2.47749896e-16 -1.88096579e-16]
 [ 2.19434836e-16  2.72066758e-16  1.00000000e+00  1.38261193e-16]
 [-1.10249066e-15 -2.10303071e-16  8.14451571e-17  1.00000000e+00]]

which means that the first two Schur vectors are not D-orthogonal, the others are fine. I check the Schur matrix to find out what these first two correspond to:

array([[ 1.     , -0.     , -0.06543,  0.11101],
       [ 0.     ,  1.     , -0.03596,  0.07668],
       [ 0.     ,  0.     ,  0.8871 , -0.23621],
       [ 0.     ,  0.     ,  0.     ,  0.59929]])

So they correspond to the 1-eigenvalues. There are two 1-eigenvalues here, one for each recurrent class.

Marius1311 commented 3 years ago

Okay, I'm going through the code now and I'm wondering why we're checking orthogonality wrt standart inner product and not wrt. D in 228:

if not np.allclose(Q.T.dot(Q), np.eye(Q.shape[1]), rtol=1e6*eps, atol=1e6*eps)

We should be checking this wrt. D, no? @msmdev

Marius1311 commented 3 years ago

Okay, I'm going through the code now and I'm wondering why we're checking orthogonality wrt standart inner product and not wrt. D in 228:
if not np.allclose(Q.T.dot(Q), np.eye(Q.shape[1]), rtol=1e6*eps, atol=1e6*eps)
We should be checking this wrt. D, no? @msmdev

Changing this fixed the above error with reducible matrices.

msmdev commented 3 years ago

Okay, I'm going through the code now and I'm wondering why we're checking orthogonality wrt standart inner product and not wrt. D in 228:
if not np.allclose(Q.T.dot(Q), np.eye(Q.shape[1]), rtol=1e6*eps, atol=1e6*eps)
We should be checking this wrt. D, no? @msmdev
Changing this fixed the above error with reducible matrices.

Makes perfect sense... thx!

Marius1311 commented 3 years ago

Oh, someone's back from holidays!

Marius1311 commented 3 years ago

@michalk8 , I would suggest to have two functions to check schur decompositions:

Full decomposition

If we have computed the full decomposition, we can do what you are currently doing in the tests:

def _assert_schur(
    P: np.ndarray, X: np.ndarray, RR: np.ndarray, N: Optional[int] = None
):
    if N is not None:
        np.testing.assert_array_equal(P.shape, [N, N])
        np.testing.assert_array_equal(X.shape, [N, N])
        np.testing.assert_array_equal(RR.shape, [N, N])

    assert np.all(np.abs(X @ RR - P @ X) < eps), np.abs(X @ RR - P @ X)
    assert np.all(np.abs(X[:, 0] - 1) < eps), np.abs(X[:, 0])

Partial Schur decomposition'

If we have however only computed a partial schur decomposition, the above will fail. In that case, we can use the checks we already have in do_schur:

    # Raise, if X doesn't fullfill the invariant subspace condition!
    dp = np.dot(P, sp.csr_matrix(X) if issparse(P) else X)
    dummy = subspace_angles(dp.toarray() if issparse(dp) else dp, np.dot(X, R))

    test = np.allclose(dummy, 0.0, atol=1e-6, rtol=1e-5)
    test1 = (dummy.shape[0] == m)

Marius1311 commented 3 years ago

@msmdev , is this todo something we should be worried about?

    # Search for the constant (Schur) vector, if explicitly present.
    max_i = 0
    for i in range(m):
        vsum = np.sum(X[:,i])
        dummy = ( np.ones(X[:,i].shape) * (vsum / n) )
        if np.allclose(X[:,i], dummy, rtol=1e-6, atol=1e-5 ):  
            max_i = i #TODO: check, if more than one vec fulfills this

msmdev commented 3 years ago

@msmdev , is this todo something we should be worried about?

    # Search for the constant (Schur) vector, if explicitly present.
    max_i = 0
    for i in range(m):
        vsum = np.sum(X[:,i])
        dummy = ( np.ones(X[:,i].shape) * (vsum / n) )
        if np.allclose(X[:,i], dummy, rtol=1e-6, atol=1e-5 ):  
            max_i = i #TODO: check, if more than one vec fulfills this

Not sure. This never happened to me, but I think we can't rule this case out with 100% certainty. Would need to check the Perron Frobenius theorem for non-negative matrices to find out if we need to worry about this.

Marius1311 commented 3 years ago

@msmdev , is this todo something we should be worried about?
    # Search for the constant (Schur) vector, if explicitly present.
    max_i = 0
    for i in range(m):
        vsum = np.sum(X[:,i])
        dummy = ( np.ones(X[:,i].shape) * (vsum / n) )
        if np.allclose(X[:,i], dummy, rtol=1e-6, atol=1e-5 ):  
            max_i = i #TODO: check, if more than one vec fulfills this
Not sure. This never happened to me, but I think we can't rule this case out with 100% certainty. Would need to check the Perron Frobenius theorem for non-negative matrices to find out if we need to worry about this.

Would be great if you could look into this, this is one of the last open questions here.

Marius1311 commented 3 years ago

@michalk8, I removed this check from TestCustom:

            # AssertionError: 0.0012763188973086148
            assert np.all(np.abs(X @ RR - P_2 @ X) <= eps), np.abs(X @ RR - P_2 @ X).max()

This is expected to fail for partial Schur decomposition, so it does not make sense to check this for partial Schur decompositions.

Marius1311 commented 3 years ago

@michalk8 , I added matrices P_3 - P_7 to conftest.py. Could you add some standard tests for all of these to TestCustom please? Could actually be as simple as what we currently have for test_P2. Maybe just add one test on the coarse grained transition matrix that checks whether rows sum to one.

Marius1311 commented 3 years ago

Looking forward to discussing this tomorrow! I'm optimistic we can open a PR for this soon.

Marius1311 commented 3 years ago

Thanks for adding these Mike!

Marius1311 commented 3 years ago

Hi @msmdev , did you have the chance to go through our tests yet? Would be great if we could merge this this week still!

msmdev commented 3 years ago

Hi @msmdev , did you have the chance to go through our tests yet? Would be great if we could merge this this week still!

Sorry @Marius1311 , I couldn't find the time so far, but I will do my very best to do so as fast as possible!

Marius1311 commented 3 years ago

thanks!

Marius1311 commented 3 years ago

Hi @michalk8 , please have a look through @msmdev s comments and in particular, could you implement the row-matching scheme we discussed?

michalk8 commented 3 years ago

As of now, the only failing tests are here: https://github.com/msmdev/msmtools/pull/18#discussion_r529078100

Marius1311 commented 3 years ago

@msmdev , Mike went through your comments and fixed basically everything, I think this is ready for another round of your review.

Marius1311 commented 3 years ago

Todo for @michalk8 , add test for membership vectors in the dense/sparse setting.

michalk8 commented 3 years ago

Done in 2 places: https://github.com/msmdev/msmtools/pull/18/files#diff-a77e704d8e9642584dac7cea5a986e5cf755403a1056410c6718fc75b012d381R660 and https://github.com/msmdev/msmtools/pull/18/files#diff-a77e704d8e9642584dac7cea5a986e5cf755403a1056410c6718fc75b012d381R656 Small note: I've added the test requirement (pytest-mock) because I monkeypatch 2 methods in https://github.com/msmdev/msmtools/pull/18/files#diff-a77e704d8e9642584dac7cea5a986e5cf755403a1056410c6718fc75b012d381R338

msmdev commented 3 years ago

Currently, I have an important deadline to keep, so I will try to do the review on the weekend... Please, excuse the delay!

Marius1311 commented 3 years ago

No worries, keep us updated.

Marius1311 commented 3 years ago

Hi @msmdev, what's the state of the art here? Did you manage to look into this over the weekend? I think our aim should be to open a PR to upstream this week - it does not need to be perfect for that yet. We just really have to push this so we get done before we get revisions in, otherwise, this will delay the process for very long and hinder maintainability of cellrank.

msmdev commented 3 years ago

@Marius1311 @michalk8 I tried to review the changes, but a had troubles with many of the supplied links to the changes, so I had trouble finding the related changes. Also I had very few time to do so - please excuse, if I might have misunderstood or overlooked some things. In those cases, please try to elaborate the things in more detail to help me getting to the point very fast. Thank you!

michalk8 commented 3 years ago

Hi @msmdev , sorry for the sloppy linking, I must've link the previous commit - I will link now not the changes in the PR, but rather the branch. Regarding the errors in https://github.com/msmdev/msmtools/pull/18#discussion_r533187270, they only occur in test_sort_real_schur for 3 matrices (out of 5):

# M_2
E       AssertionError: 
E       Not equal to tolerance rtol=1e-05, atol=5
E       
E       Mismatched elements: 4 / 4 (100%)
E       Max absolute difference: 296.85321743
E       Max relative difference: 296.85321743
E        x: array([297.853217, 297.853217, 297.667196, 297.667196])
E        y: array(1.)

# M_3
E       AssertionError: 
E       Not equal to tolerance rtol=1e-05, atol=5
E       
E       Mismatched elements: 4 / 4 (100%)
E       Max absolute difference: 317.45258363
E       Max relative difference: 317.45258363
E        x: array([318.452584, 318.452584, 318.293397, 318.293397])
E        y: array(1.)

# M_5
E       AssertionError: 
E       Not equal to tolerance rtol=1e-05, atol=5
E       
E       Mismatched elements: 4 / 4 (100%)
E       Max absolute difference: 26991.94731741
E       Max relative difference: 26991.94731741
E        x: array([26992.947317, 26992.947317, 26022.44675 , 26022.44675 ])
E        y: array(1.)

ping @Marius1311 https://github.com/msmdev/msmtools/pull/18#discussion_r533223339 - is this fine for us?

Marius1311 commented 3 years ago

Hi @msmdev , sorry for the sloppy linking, I must've link the previous commit - I will link now not the changes in the PR, but rather the branch. Regarding the errors in #18 (comment), they only occur in test_sort_real_schur for 3 matrices (out of 5):

# M_2
E       AssertionError: 
E       Not equal to tolerance rtol=1e-05, atol=5
E       
E       Mismatched elements: 4 / 4 (100%)
E       Max absolute difference: 296.85321743
E       Max relative difference: 296.85321743
E        x: array([297.853217, 297.853217, 297.667196, 297.667196])
E        y: array(1.)

# M_3
E       AssertionError: 
E       Not equal to tolerance rtol=1e-05, atol=5
E       
E       Mismatched elements: 4 / 4 (100%)
E       Max absolute difference: 317.45258363
E       Max relative difference: 317.45258363
E        x: array([318.452584, 318.452584, 318.293397, 318.293397])
E        y: array(1.)

# M_5
E       AssertionError: 
E       Not equal to tolerance rtol=1e-05, atol=5
E       
E       Mismatched elements: 4 / 4 (100%)
E       Max absolute difference: 26991.94731741
E       Max relative difference: 26991.94731741
E        x: array([26992.947317, 26992.947317, 26022.44675 , 26022.44675 ])
E        y: array(1.)

ping @Marius1311 #18 (comment) - is this fine for us?

Hi @michalk8 , I still have no idea what this test is testing for, so I can't say whether it's okay for us whether it fails. @msmdev , could you please explain the logic of this test? To facilitate this, here's the code:

    def test_sort_real_schur(self, R_i: np.ndarray):
        # test_SRSchur_num_t
        Q = np.eye(4)
        QQ, RR, ap = sort_real_schur(Q, R_i, z="LM", b=0)

        assert np.all(np.array(ap) <= 1), ap

        EQ = np.true_divide(np.linalg.norm(Q - QQ.T @ QQ, ord=1), eps)
        assert_allclose(EQ, 1.0, atol=5)

        EA = np.true_divide(np.linalg.norm(R_i - QQ @ RR @ QQ.T, ord=1), eps * np.linalg.norm(R_i, ord=1))
        assert_allclose(EA, 1.0, atol=5)

        l = eigs(R_i, k=4, which="SM", return_eigenvectors=False)
        ll = eigs(RR, k=4, which="SM", return_eigenvectors=False)
        EL = np.true_divide(np.abs(l - ll), eps * np.abs(l))
        assert_allclose(EL, 1.0, atol=5)

According to @michalk8 , we're failing at the last assertion, which checks whether the smallest eigenvalues of these matrices are similar, if I understand this correctly.

michalk8 commented 3 years ago

@msmdev regarding the discrepancy between sparse/dense The rotation matrix is similar in absolute difference (10e-2), just the signs are sometimes flipped and when comparing the memberships from Matlab's ground truth, the max. absolute difference is less than 1-e4. https://github.com/michalk8/msmtools/commit/212cb240f97333039ecd0b938c0ef8579f80b42b#diff-a77e704d8e9642584dac7cea5a986e5cf755403a1056410c6718fc75b012d381R619

msmdev commented 3 years ago

Hi @msmdev , sorry for the sloppy linking, I must've link the previous commit - I will link now not the changes in the PR, but rather the branch. Regarding the errors in #18 (comment), they only occur in test_sort_real_schur for 3 matrices (out of 5):

# M_2
E       AssertionError: 
E       Not equal to tolerance rtol=1e-05, atol=5
E       
E       Mismatched elements: 4 / 4 (100%)
E       Max absolute difference: 296.85321743
E       Max relative difference: 296.85321743
E        x: array([297.853217, 297.853217, 297.667196, 297.667196])
E        y: array(1.)

# M_3
E       AssertionError: 
E       Not equal to tolerance rtol=1e-05, atol=5
E       
E       Mismatched elements: 4 / 4 (100%)
E       Max absolute difference: 317.45258363
E       Max relative difference: 317.45258363
E        x: array([318.452584, 318.452584, 318.293397, 318.293397])
E        y: array(1.)

# M_5
E       AssertionError: 
E       Not equal to tolerance rtol=1e-05, atol=5
E       
E       Mismatched elements: 4 / 4 (100%)
E       Max absolute difference: 26991.94731741
E       Max relative difference: 26991.94731741
E        x: array([26992.947317, 26992.947317, 26022.44675 , 26022.44675 ])
E        y: array(1.)

ping @Marius1311 #18 (comment) - is this fine for us?

Hi @michalk8 , I still have no idea what this test is testing for, so I can't say whether it's okay for us whether it fails. @msmdev , could you please explain the logic of this test? To facilitate this, here's the code:

    def test_sort_real_schur(self, R_i: np.ndarray):
        # test_SRSchur_num_t
        Q = np.eye(4)
        QQ, RR, ap = sort_real_schur(Q, R_i, z="LM", b=0)

        assert np.all(np.array(ap) <= 1), ap

        EQ = np.true_divide(np.linalg.norm(Q - QQ.T @ QQ, ord=1), eps)
        assert_allclose(EQ, 1.0, atol=5)

        EA = np.true_divide(np.linalg.norm(R_i - QQ @ RR @ QQ.T, ord=1), eps * np.linalg.norm(R_i, ord=1))
        assert_allclose(EA, 1.0, atol=5)

        l = eigs(R_i, k=4, which="SM", return_eigenvectors=False)
        ll = eigs(RR, k=4, which="SM", return_eigenvectors=False)
        EL = np.true_divide(np.abs(l - ll), eps * np.abs(l))
        assert_allclose(EL, 1.0, atol=5)

According to @michalk8 , we're failing at the last assertion, which checks whether the smallest eigenvalues of these matrices are similar, if I understand this correctly.

@Marius1311 @michalk8 Okay, this is not so trivial: those tests are taken from Brandts original publication of the SRSchur sorting algorithm: They test normal execution with test cases 1-4 from Table I in (1) Brandts, J. H. Matlab Code for Sorting Real Schur Forms. Numer. Linear Algebr. with Appl. 2002, 9 (3), 249-261 and a fifth test case designed from the first test matrix M1 by swapping of the first block A11 with the second block A22. The formula for EA had a typo in (1), the correct formula was taken from (2) Bai, Z.; Demmel, J. W. On Swapping Diagonal Blocks in Real Schur Form. Linear Algebra Appl. 1993, 186 (C), 75-95.

Please take a look there to get an impression of what we are testing here... Anyway, I would consider this test important and it shouldn't fail so badly... Just one idea: you are using scipy.sparse.linalg.eigs to get the eigenvalues, right? This might again (I already mentioned this some time ago...) not be such a good idea .

michalk8 commented 3 years ago

Thanks a lot @msmdev , I completely forgot about that, switching to eigvals solved it! I'm happy to confirm that all the tests now pass. @Marius1311 @msmdev would be great if both of you could confirm this by running pytest tests/analysis/impl/dense/gpcca_test.py You might need to pip install -e.[test], since I've added pytest-mock requirement (needed for this 1 test): https://github.com/michalk8/msmtools/blob/gpcca_tests/tests/analysis/impl/dense/gpcca_test.py#L339

Marius1311 commented 3 years ago

@msmdev regarding the discrepancy between sparse/dense The rotation matrix is similar in absolute difference (10e-2), just the signs are sometimes flipped and when comparing the memberships from Matlab's ground truth, the max. absolute difference is less than 1-e4. michalk8@212cb24#diff-a77e704d8e9642584dac7cea5a986e5cf755403a1056410c6718fc75b012d381R619

Amazing, thanks for looking into this @michalk8 . This is as I expected: the rotation matrix differs, but the final outcome in terms of membership vectors is the same. There are more things we could check for here - however, I would prioritise closing this PR and opening the upstream PR. This is more important from my point of view right now as the current situation is really a bit messy, and this PR has gotten very large already.

Marius1311 commented 3 years ago

Thanks a lot @msmdev , I completely forgot about that, switching to eigvals solved it! I'm happy to confirm that all the tests now pass. @Marius1311 @msmdev would be great if both of you could confirm this by running pytest tests/analysis/impl/dense/gpcca_test.py You might need to pip install -e.[test], since I've added pytest-mock requirement (needed for this 1 test): https://github.com/michalk8/msmtools/blob/gpcca_tests/tests/analysis/impl/dense/gpcca_test.py#L339

Thanks @msmdev for pointing this out, and thanks to @michalk8 for solving this! I will attempt to run the tests locally on my machine now.

Marius1311 commented 3 years ago

pip install -e.[test] on mac gives zsh: no matches found: -e.[test]. Next, tried pip install -e '.[test]'. This works, but says WARNING: msmtools 1.2.5+273.gf6e93f8 does not provide the extra 'test'. I pulled before doing this, so should be up to date.

Marius1311 commented 3 years ago

Here's the output I get upon running pytest:

==================================================================================== short test summary info =====================================================================================
FAILED tests/analysis/test_mfpt.py::TestMfptSparse::test_mfpt_between_sets - TypeError: Cannot use scipy.linalg.eig for sparse A with k >= N - 1. Use scipy.linalg.eig(A.toarray()) or reduce k.
FAILED tests/analysis/impl/dense/gpcca_test.py::TestPETScSLEPc::test_gpcca_krylov_sparse_eq_dense[0] - AssertionError: 
FAILED tests/analysis/impl/dense/gpcca_test.py::TestPETScSLEPc::test_gpcca_krylov_sparse_eq_dense[10] - AssertionError: 
FAILED tests/analysis/impl/dense/gpcca_test.py::TestPETScSLEPc::test_gpcca_krylov_sparse_eq_dense[200] - AssertionError: 
FAILED tests/analysis/impl/dense/gpcca_test.py::TestPETScSLEPc::test_gpcca_krylov_sparse_eq_dense[500] - AssertionError: 
FAILED tests/analysis/impl/dense/gpcca_test.py::TestPETScSLEPc::test_gpcca_krylov_sparse_eq_dense[1000] - AssertionError: 
FAILED tests/analysis/impl/sparse/mean_first_passage_time_test.py::TestMfpt::test_mfpt_between_sets - TypeError: Cannot use scipy.linalg.eig for sparse A with k >= N - 1. Use scipy.linalg.eig...
FAILED tests/estimation/tests/test_effective_count_matrix.py::TestEffectiveCountMatrix::test_njobs_speedup - AssertionError: 1.5465750694274902 not less than 1.368762493133545 : does not scal...
FAILED tests/estimation/tests/test_transition_matrix.py::TestTransitionRevPiSym::test_sparse - TypeError: Cannot use scipy.linalg.eig for sparse A with k >= N - 1. Use scipy.linalg.eig(A.toar...
ERROR tests/analysis/impl/dense/gpcca_test.py::TestGPCCAMatlabUnit::test_objective_1st_col

The ones relevant for us, I guess, are

(I) sparse/dense mismatches

FAILED tests/analysis/impl/dense/gpcca_test.py::TestPETScSLEPc::test_gpcca_krylov_sparse_eq_dense[0] - AssertionError: 
FAILED tests/analysis/impl/dense/gpcca_test.py::TestPETScSLEPc::test_gpcca_krylov_sparse_eq_dense[10] - AssertionError: 
FAILED tests/analysis/impl/dense/gpcca_test.py::TestPETScSLEPc::test_gpcca_krylov_sparse_eq_dense[200] - AssertionError: 
FAILED tests/analysis/impl/dense/gpcca_test.py::TestPETScSLEPc::test_gpcca_krylov_sparse_eq_dense[500] - AssertionError: 
FAILED tests/analysis/impl/dense/gpcca_test.py::TestPETScSLEPc::test_gpcca_krylov_sparse_eq_dense[1000] - AssertionError:

(II) something with the 1st column

ERROR tests/analysis/impl/dense/gpcca_test.py::TestGPCCAMatlabUnit::test_objective_1st_col

Marius1311 commented 3 years ago

(I) sparse/dense mismatches

For the first class (I) of failing tests, here's the code:

    def test_gpcca_krylov_sparse_eq_dense(self, example_matrix_mu: np.ndarray):
        P, sd = get_known_input(example_matrix_mu)

        # for 3 it's fine
        g_s = GPCCA(csr_matrix(P), eta=sd, method="krylov").optimize(4)
        g_d = GPCCA(P, eta=sd, method="krylov").optimize(4)

        assert issparse(g_s.P)
        assert not issparse(g_d.P)

        assert_allclose(g_s.memberships.sum(1), 1.0)
        assert_allclose(g_d.memberships.sum(1), 1.0)

        X_k, X_b = g_s.schur_vectors, g_d.schur_vectors
        RR_k, RR_b = g_s.schur_matrix, g_d.schur_matrix

        # check if it's a correct Schur form
        _assert_schur(P, X_b, RR_b, N=None)
        _assert_schur(P, X_k, RR_k, N=None)
        # check if they span the same subspace
        assert np.max(subspace_angles(X_b, X_k)) < eps

        ms, md = g_s.memberships, g_d.memberships
        cs, cd = g_s.coarse_grained_transition_matrix, g_d.coarse_grained_transition_matrix
        perm = _find_permutation(md, ms)

        cs = cs[:, perm]
>      assert_allclose(cs, g_d.coarse_grained_transition_matrix)

        ms = ms[:, perm]
        assert_allclose(ms, md)

We're failing at the indicated assertion, which checks for equality of coarse grained transition probs. Looking into the first of these failing tests in class (I) for \mu = 0, I get the output

E       AssertionError: 
E       Not equal to tolerance rtol=1e-05, atol=1e-08
E       
E       Mismatched elements: 16 / 16 (100%)
E       Max absolute difference: 0.99858631
E       Max relative difference: 348.37650063
E        x: array([[ 7.603764e-01,  8.171464e-02,  7.559078e-02,  8.231817e-02],
E              [ 8.063700e-03, -2.684025e-03,  9.948072e-01, -1.868863e-04],
E              [ 8.651307e-03,  9.944040e-01, -2.670513e-03, -3.847842e-04],
E              [ 8.805835e-03, -4.138702e-04, -1.865480e-04,  9.917946e-01]])
E        y: array([[ 7.613608e-01,  7.934900e-02,  9.091987e-02,  6.837028e-02],
E              [ 1.013324e-02,  9.939923e-01, -3.779101e-03, -3.464779e-04],
E              [ 8.836629e-03, -2.862612e-03,  9.940104e-01,  1.561961e-05],
E              [ 9.662125e-03, -6.172231e-04, -1.063530e-03,  9.920186e-01]])

The problem is that the vectors are not permuted correctly. Permuting 2 and 3 (start counting from 1) in the first matrix resolves the large absolute deviation.

Marius1311 commented 3 years ago

@michalk8, the problem for class (I) failures is the permutation scheme. Could you look into it please?

michalk8 commented 3 years ago

@Marius1311 I forgot to do the permutation in 1 place + also permute the rows of coarse-grained TM, could you please check again?

Marius1311 commented 3 years ago

(II) something with the first column

Here's the code of the testing function

 def test_objective_1st_col(self, mocker):
        # check_in_matlab: _objective
        P, _ = get_known_input(mu(0))
        N, M = P.shape[0], 4

        _, L, _ = lu(P[:, :M])
        mocker.patch(
            "msmtools.util.sorted_schur.sorted_schur",
            return_value=(np.eye(M), L, np.array([np.nan] * M)),
        )
        mocker.patch("msmtools.analysis.dense.gpcca._gram_schmidt_mod", return_value=L)

        with pytest.raises(
            ValueError,
            match=r"The first column X\[:, 0\] of the Schur "
            r"vector matrix isn't constantly equal 1.",
        ):
            _do_schur(P, eta=np.true_divide(np.ones((N,), dtype=np.float64), N), m=M)

This was a problem with me missing the pytest-mock dependency. When I manually installed it and run the tests again using pytest, I got

==================================================================================== short test summary info =====================================================================================
FAILED tests/analysis/test_mfpt.py::TestMfptSparse::test_mfpt_between_sets - TypeError: Cannot use scipy.linalg.eig for sparse A with k >= N - 1. Use scipy.linalg.eig(A.toarray()) or reduce k.
FAILED tests/analysis/impl/dense/gpcca_test.py::TestPETScSLEPc::test_gpcca_krylov_sparse_eq_dense[0] - AssertionError: 
FAILED tests/analysis/impl/dense/gpcca_test.py::TestPETScSLEPc::test_gpcca_krylov_sparse_eq_dense[10] - AssertionError: 
FAILED tests/analysis/impl/dense/gpcca_test.py::TestPETScSLEPc::test_gpcca_krylov_sparse_eq_dense[200] - AssertionError: 
FAILED tests/analysis/impl/dense/gpcca_test.py::TestPETScSLEPc::test_gpcca_krylov_sparse_eq_dense[500] - AssertionError: 
FAILED tests/analysis/impl/dense/gpcca_test.py::TestPETScSLEPc::test_gpcca_krylov_sparse_eq_dense[1000] - AssertionError: 
FAILED tests/analysis/impl/sparse/mean_first_passage_time_test.py::TestMfpt::test_mfpt_between_sets - TypeError: Cannot use scipy.linalg.eig for sparse A with k >= N - 1. Use scipy.linalg.eig...
FAILED tests/estimation/tests/test_effective_count_matrix.py::TestEffectiveCountMatrix::test_njobs_speedup - AssertionError: 1.9440791606903076 not less than 1.553575038909912 : does not scal...
FAILED tests/estimation/tests/test_transition_matrix.py::TestTransitionRevPiSym::test_sparse - TypeError: Cannot use scipy.linalg.eig for sparse A with k >= N - 1. Use scipy.linalg.eig(A.toar...
========================================================================== 9 failed, 387 passed, 74 warnings in 45.66s ===========================================================================

i.e. the Error was gone, problem solved.

Marius1311 commented 3 years ago

@Marius1311 I forgot to do the permutation in 1 place + also permute the rows of coarse-grained TM, could you please check again?

I @michalk8 , just pulled and re-run pytest, the problem persits, i.e. all class (I) test still fail. For the first matrix, the output is

E       AssertionError: 
E       Not equal to tolerance rtol=1e-05, atol=1e-08
E       
E       Mismatched elements: 16 / 16 (100%)
E       Max absolute difference: 0.01532909
E       Max relative difference: 12.9648493
E        x: array([[ 7.603764e-01,  8.171464e-02,  7.559078e-02,  8.231817e-02],
E              [ 8.651307e-03,  9.944040e-01, -2.670513e-03, -3.847842e-04],
E              [ 8.063700e-03, -2.684025e-03,  9.948072e-01, -1.868863e-04],
E              [ 8.805835e-03, -4.138702e-04, -1.865480e-04,  9.917946e-01]])
E        y: array([[ 7.613608e-01,  7.934900e-02,  9.091987e-02,  6.837028e-02],
E              [ 1.013324e-02,  9.939923e-01, -3.779101e-03, -3.464779e-04],
E              [ 8.836629e-03, -2.862612e-03,  9.940104e-01,  1.561961e-05],
E              [ 9.662125e-03, -6.172231e-04, -1.063530e-03,  9.920186e-01]])

msmdev / msmtools

Gpcca tests #18

Reducible matrix

The matrix

The problem

Full decomposition

Partial Schur decomposition'

(I) sparse/dense mismatches

(II) something with the 1st column

(I) sparse/dense mismatches

(II) something with the first column