Write tests for core functions

[j-i-l]

[j-i-l]

We have some issues with the type declarations for indices and indptr:

https://github.com/alexbovet/flow_stability/actions/runs/10513029536/job/29127746857#step:5:116

It also seems not completely consistent in the defs:

E.g. for indptr some functions require int:

https://github.com/alexbovet/flow_stability/blob/8648a8bd715e3b086ed13cb7453d41754c176398/src/flowstab/_cython_sparse_stoch.pyx#L859-L863

some long long: https://github.com/alexbovet/flow_stability/blob/8648a8bd715e3b086ed13cb7453d41754c176398/src/flowstab/_cython_sparse_stoch.pyx#L613-L618

I'd put everything indices and indptr to long and declare with np.int64 on the python site. @alexbovet is that OK for you or is there actually a use-case for long long?

[alexbovet]

I know that I had to switch to long long to allow for more than 2**31 -1 elements for a project. I also had to set MKL_INTERFACE_LAYER=ILP64 for sparse_dot_mkl... But this is not necessary for all projects and probably wastes memory, so an option somewhere would be nice.

[github-actions[bot]]

Coverage report

Click to see where and how coverage changed

File	Statements	Missing	Coverage	Coverage (new stmts)	Lines missing
src/flowstab
FlowStability.py
SparseStochMat.py					1009
SynthTempNetwork.py
_cython_sparse_stoch_subst.py					37-40
Project Total

This report was generated by python-coverage-comment-action

[j-i-l]

Starting to add some tests for the usage of mkl:

For dot_product_mkl it seems not obvious that there is a speed benefit as compared to scipy's @.

Using very sparse random matrices: https://github.com/alexbovet/flow_stability/blob/3e15f1f640cab87cb5fe39087a23f15e12e5b659/tests/conftest.py#L31-L41

and comparing timing and memory usage with:

https://github.com/alexbovet/flow_stability/blob/3e15f1f640cab87cb5fe39087a23f15e12e5b659/tests/test_sparse_stoch_matrix.py#L204-L217

we get:

================================ MEMRAY REPORT =================================
Allocation results for tests/test_sparse_stoch_matrix.py::test_sparse_matmul_memory at the high watermark

	 📦 Total memory allocated: 152.6MiB
	 📏 Total allocations: 6
	 📊 Histogram of allocation sizes: |█   ▄|
	 🥇 Biggest allocating functions:
		- _matmul_sparse:/home/runner/work/flow_stability/flow_stability/venv/lib/python3.9/site-packages/scipy/sparse/_compressed.py:535 -> 114.4MiB
		- _matmul_sparse:/home/runner/work/flow_stability/flow_stability/venv/lib/python3.9/site-packages/scipy/sparse/_compressed.py:530 -> 38.1MiB
		- __init__:/home/runner/work/flow_stability/flow_stability/venv/lib/python3.9/site-packages/scipy/sparse/_compressed.py:27 -> 992.0B
		- check_shape:/home/runner/work/flow_stability/flow_stability/venv/lib/python3.9/site-packages/scipy/sparse/_sputils.py:296 -> 864.0B
		- check_format:/home/runner/work/flow_stability/flow_stability/venv/lib/python3.9/site-packages/scipy/sparse/_compressed.py:132 -> 608.0B


Allocation results for tests/test_sparse_stoch_matrix.py::test_sparse_matmul_mkl_memory at the high watermark

	 📦 Total memory allocated: 114.5MiB
	 📏 Total allocations: 16
	 📊 Histogram of allocation sizes: |▃█▁  |
	 🥇 Biggest allocating functions:
		- _matmul_mkl:/home/runner/work/flow_stability/flow_stability/venv/lib/python3.9/site-packages/sparse_dot_mkl/_sparse_sparse.py:35 -> 114.5MiB
		- _pass_mkl_handle_csr_csc:/home/runner/work/flow_stability/flow_stability/venv/lib/python3.9/site-packages/sparse_dot_mkl/_mkl_interface/_common.py:307 -> 16.6KiB
		- _pass_mkl_handle_csr_csc:/home/runner/work/flow_stability/flow_stability/venv/lib/python3.9/site-packages/sparse_dot_mkl/_mkl_interface/_common.py:307 -> 16.6KiB
		- _ctype_ndarray:/home/runner/work/flow_stability/flow_stability/venv/lib/python3.9/site-packages/numpy/ctypeslib.py:354 -> 936.0B
		- as_array:/home/runner/work/flow_stability/flow_stability/venv/lib/python3.9/site-packages/numpy/ctypeslib.py:521 -> 936.0B

...

============================== slowest durations ===============================
48.14s call     tests/test_sparse_stoch_matrix.py::test_sparse_matmul_mkl_memory
46.04s call     tests/test_sparse_stoch_matrix.py::test_sparse_matmul_memory
...

So, the mkl based approach is even slightly slower.

One thing to note here, though, is that the github runner has probably only few threads (if even more than just 1).

[j-i-l]

In a more multi-threading friendly environment (i.e. my laptop) things look indeed different:

==================================== MEMRAY REPORT =====================================
Allocation results for tests/test_sparse_stoch_matrix.py::test_sparse_matmul_mkl_memory at the high watermark

	 📦 Total memory allocated: 190.8MiB
	 📏 Total allocations: 12
	 📊 Histogram of allocation sizes: |█▁   |
	 🥇 Biggest allocating functions:
		- _matmul_mkl:/home/jonas/Projects/FlowStability/flow_stability/venv/lib/python3.9/site-packages/sparse_dot_mkl/_sparse_sparse.py:35 -> 190.8MiB
		- _pass_mkl_handle_csr_csc:/home/jonas/Projects/FlowStability/flow_stability/venv/lib/python3.9/site-packages/sparse_dot_mkl/_mkl_interface/_common.py:307 -> 16.7KiB
		- _pass_mkl_handle_csr_csc:/home/jonas/Projects/FlowStability/flow_stability/venv/lib/python3.9/site-packages/sparse_dot_mkl/_mkl_interface/_common.py:307 -> 16.7KiB
		- _ctype_ndarray:/home/jonas/Projects/FlowStability/flow_stability/venv/lib/python3.9/site-packages/numpy/ctypeslib.py:354 -> 936.0B
		- _ctype_ndarray:/home/jonas/Projects/FlowStability/flow_stability/venv/lib/python3.9/site-packages/numpy/ctypeslib.py:354 -> 936.0B


Allocation results for tests/test_sparse_stoch_matrix.py::test_sparse_matmul_memory at the high watermark

	 📦 Total memory allocated: 152.6MiB
	 📏 Total allocations: 7
	 📊 Histogram of allocation sizes: |█▂  ▂|
	 🥇 Biggest allocating functions:
		- _matmul_sparse:/home/jonas/Projects/FlowStability/flow_stability/venv/lib/python3.9/site-packages/scipy/sparse/_compressed.py:535 -> 114.4MiB
		- _matmul_sparse:/home/jonas/Projects/FlowStability/flow_stability/venv/lib/python3.9/site-packages/scipy/sparse/_compressed.py:530 -> 38.1MiB
		- _matmul_sparse:/home/jonas/Projects/FlowStability/flow_stability/venv/lib/python3.9/site-packages/scipy/sparse/_compressed.py:532 -> 8.1KiB
		- _matmul_sparse:/home/jonas/Projects/FlowStability/flow_stability/venv/lib/python3.9/site-packages/scipy/sparse/_compressed.py:531 -> 4.0KiB
		- __init__:/home/jonas/Projects/FlowStability/flow_stability/venv/lib/python3.9/site-packages/scipy/sparse/_compressed.py:27 -> 992.0B

================================== slowest durations ===================================
145.27s call     tests/test_sparse_stoch_matrix.py::test_sparse_matmul_memory
108.48s call     tests/test_sparse_stoch_matrix.py::test_sparse_matmul_mkl_memory
...

Which looks more like expected.

One thing that bugs me with the current approach of using MKL is that MKL also provides BLAS and LAPACK and thus could be linked against when building scipy. In doing so we should be able to get a scipy version that relies on MKL for matrix operations and we would not need to write any extra code and checks if some package is installed or not.

[j-i-l]

In a more multi-threading friendly environment (i.e. my laptop) things look indeed different:

==================================== MEMRAY REPORT =====================================
Allocation results for tests/test_sparse_stoch_matrix.py::test_sparse_matmul_mkl_memory at the high watermark

	 📦 Total memory allocated: 190.8MiB
	 📏 Total allocations: 12
	 📊 Histogram of allocation sizes: |█▁   |
	 🥇 Biggest allocating functions:
		- _matmul_mkl:/home/jonas/Projects/FlowStability/flow_stability/venv/lib/python3.9/site-packages/sparse_dot_mkl/_sparse_sparse.py:35 -> 190.8MiB
		- _pass_mkl_handle_csr_csc:/home/jonas/Projects/FlowStability/flow_stability/venv/lib/python3.9/site-packages/sparse_dot_mkl/_mkl_interface/_common.py:307 -> 16.7KiB
		- _pass_mkl_handle_csr_csc:/home/jonas/Projects/FlowStability/flow_stability/venv/lib/python3.9/site-packages/sparse_dot_mkl/_mkl_interface/_common.py:307 -> 16.7KiB
		- _ctype_ndarray:/home/jonas/Projects/FlowStability/flow_stability/venv/lib/python3.9/site-packages/numpy/ctypeslib.py:354 -> 936.0B
		- _ctype_ndarray:/home/jonas/Projects/FlowStability/flow_stability/venv/lib/python3.9/site-packages/numpy/ctypeslib.py:354 -> 936.0B


Allocation results for tests/test_sparse_stoch_matrix.py::test_sparse_matmul_memory at the high watermark

	 📦 Total memory allocated: 152.6MiB
	 📏 Total allocations: 7
	 📊 Histogram of allocation sizes: |█▂  ▂|
	 🥇 Biggest allocating functions:
		- _matmul_sparse:/home/jonas/Projects/FlowStability/flow_stability/venv/lib/python3.9/site-packages/scipy/sparse/_compressed.py:535 -> 114.4MiB
		- _matmul_sparse:/home/jonas/Projects/FlowStability/flow_stability/venv/lib/python3.9/site-packages/scipy/sparse/_compressed.py:530 -> 38.1MiB
		- _matmul_sparse:/home/jonas/Projects/FlowStability/flow_stability/venv/lib/python3.9/site-packages/scipy/sparse/_compressed.py:532 -> 8.1KiB
		- _matmul_sparse:/home/jonas/Projects/FlowStability/flow_stability/venv/lib/python3.9/site-packages/scipy/sparse/_compressed.py:531 -> 4.0KiB
		- __init__:/home/jonas/Projects/FlowStability/flow_stability/venv/lib/python3.9/site-packages/scipy/sparse/_compressed.py:27 -> 992.0B

================================== slowest durations ===================================
145.27s call     tests/test_sparse_stoch_matrix.py::test_sparse_matmul_memory
108.48s call     tests/test_sparse_stoch_matrix.py::test_sparse_matmul_mkl_memory
...

Which looks more like expected.

One thing that bugs me with the current approach of using MKL is that MKL also provides BLAS and LAPACK and thus could be linked against when building scipy. In doing so we should be able to get a scipy version that relies on MKL for matrix operations and we would not need to write any extra code and checks if some package is installed or not.

One approach might be to remove the sparse_dot_mkl package and add instructions (or simply links) to how one can build scipy with MKL for BLAS and LAPAK. For now (basically as long as dependency groups are not really a thing) it would be rather cumbersome to provide a "[mkl]" installation option but we could even add a simple 'how-to' section if using the intel based scipy from the anaconda mirror, or the inter-scipy python package provide functional, pre-compiled versions of scipy with mkl.

[j-i-l]

Comparing the output of the csr_ methods with native implementations we get mismatches on .data and .indices (see e.g. https://github.com/alexbovet/flow_stability/actions/runs/10865864561/job/30152834847#step:6:102).

A smaller example illustrates this better (also comparing csr_add to a native addition):

# data
[0.16911086 0.26525225 0.04051438 0.09085281 0.55254878 0.82202173
 0.51055564 0.50439578 0.85289212 0.74480397]  # native
[0.16911086 0.26525225 0.04051438 0.09085281 0.55254878 0.82202173
 0.51055564 0.85289212 0.50439578 0.74480397]  # csr_add
# indicies
[0 1 1 2 3 3 4 5 7 9]  # native
[0 1 1 2 3 3 4 7 5 9]  # csr_add
# indptr
[ 0  0  1  2  2  4  5  5  7 10 10]  # native
[ 0  0  1  2  2  4  5  5  7 10 10]  # csr_add

Basically, the indices are not ordered within a indptr slice but data match the indices (at least in this example). Not sure if that leads to any trouble down the road.

[alexbovet]

Comparing the output of the csr_ methods with native implementations we get mismatches on .data and .indices (see e.g. https://github.com/alexbovet/flow_stability/actions/runs/10865864561/job/30152834847#step:6:102).

A smaller example illustrates this better (also comparing csr_add to a native addition):

# data
[0.16911086 0.26525225 0.04051438 0.09085281 0.55254878 0.82202173
 0.51055564 0.50439578 0.85289212 0.74480397]  # native
[0.16911086 0.26525225 0.04051438 0.09085281 0.55254878 0.82202173
 0.51055564 0.85289212 0.50439578 0.74480397]  # csr_add
# indicies
[0 1 1 2 3 3 4 5 7 9]  # native
[0 1 1 2 3 3 4 7 5 9]  # csr_add
# indptr
[ 0  0  1  2  2  4  5  5  7 10 10]  # native
[ 0  0  1  2  2  4  5  5  7 10 10]  # csr_add

Basically, the indices are not ordered within a indptr slice but data match the indices (at least in this example). Not sure if that leads to any trouble down the road.

It is expected that the indices may not be sorted after some operations. csr matrices have a has_sorted_indices attribute. If you do a sort_indices() and eliminate_zeros() before the np.testing.assert_allclose(), it should hopefully fix the problem.

[j-i-l]

Follow up from the last tests: https://github.com/alexbovet/flow_stability/actions/runs/10902255803/job/30253907061#step:6:789:

Generating sparse_autocov_mat with from_T_forward generates an object which stores in PTT^{inv} in the attribute PTwhich might be a bit confusing:

https://github.com/alexbovet/flow_stability/blob/8648a8bd715e3b086ed13cb7453d41754c176398/src/flowstab/SparseStochMat.py#L1596-L1607

We could adapt the class to keep only T, p1 and p2 as first-class citizens, populates _PT and/or _PTTinfas attributes (on demand) and defines the properties PT and PTTinf (that compute _PT and _PTTinf on the first call and return them). Something along these lines. If we keep it the way it is now, then one might add a comment that the attribute PT does not necessarily hold always diag(p1) @ T.

---

One other thing to note for sparse_autocov_mat:

Handling the option of p1 and p2 being scalars leads to a lot of noise of the form:

def some_method(self,...):
    if self.p_scalar:
        # do something with self.p1p2
        pass
    else:
        # do something with self.p1 and self.p2
        pass

It would make this class easier to read if self.p1p2 were a custom object that handles the scalar check internally and exposes just the necessary methods, like a toarray method:

def toarray(self,):  # a method of P1P2 class
    if self.scalars:
        return np.ones(self.shape)*self.p1*self.p2
    else:
        return np.outer(self.p1, self.p2)

Then we do not need to bother in each method of sparse_autocov_mat if self.p1p2 actually contains just scalars or vectors. Also we could expose p1 and p2 and keep them as properties in sparse_autocov_mat (like return self.p1p2.p1) which would always return a np.array.

Or, we simply always expand self.p1 and self.p2 (self here being an object of sparse_autocov_mat) to np.arrays which will increase the memory footprint but could be faster as we avoid object initiations we might have when using self.p1p2 (like e.g. in the aggregate method).

[alexbovet]

Generating sparse_autocov_mat with from_T_forward generates an object which stores in PTT^{inv} in the attribute PTwhich might be a bit confusing:

The standard use case is to initialize a sparse_autocov_mat as sparse_autocov_mat(PT, p1, p2), and here, the PT in general will be computed as a PTTinv. The from_T and from_T_forward methods are helpers that will probably not really be used that much.

It would be better to keep it like this but clarify the class description to explain that it can be used differently.

[j-i-l]

Howdie @alexbovet :wave:

This PR implements quite a few test and test setups. It also implements a monitoring of the code base that tracks test coverage of the existing code and of Pull Requests.

The test coverage: https://github.com/alexbovet/flow_stability/pull/35#issuecomment-2324894373 is still quite low, but for most sub-modules there are now corresponding test scripts to be found under tests/ that contain for some classes and functions some unit-test and for other only empty test functions that need to be filled eventually. However, writing unit tests for all functions is quite time consuming (writing them from scratch or using something like llama or yi-coder takes about the same amount of time - correcting the errors and inconsistencies of the llms is also quite time consuming).

Thus, I'd suggest we merge this PR and will try to slowly increase test-coverage as we work on specific issues and features.

Other than the unit tests, there are quite a few type hinting and documentation edits in this PR. The PR does not contain functional changes, so it should be save to merge it.

[j-i-l]

will need adaptations after #49 and #50