diff --git a/cooltools/api/saddle.py b/cooltools/api/saddle.py index 1c9499e1..f146fc7b 100644 --- a/cooltools/api/saddle.py +++ b/cooltools/api/saddle.py @@ -784,6 +784,18 @@ def saddle_strength(S, C): raise ValueError("`saddledata` should be square.") ratios = np.zeros(n) + + # extent=0 reduces to a single pair of corner bins on each side + intra0_sum = np.nansum(S[0, 0]) + np.nansum(S[-1, -1]) + intra0_count = np.nansum(C[0, 0]) + np.nansum(C[-1, -1]) + intra0 = intra0_sum / intra0_count + + inter0_sum = np.nansum(S[0, -1]) + np.nansum(S[-1, 0]) + inter0_count = np.nansum(C[0, -1]) + np.nansum(C[-1, 0]) + inter0 = inter0_sum / inter0_count + + ratios[0] = intra0 / inter0 + for k in range(1, n): intra_sum = np.nansum(S[0:k, 0:k]) + np.nansum(S[n - k : n, n - k : n]) intra_count = np.nansum(C[0:k, 0:k]) + np.nansum(C[n - k : n, n - k : n]) diff --git a/tests/test_compartments_saddle.py b/tests/test_compartments_saddle.py index d125aa62..85375f28 100644 --- a/tests/test_compartments_saddle.py +++ b/tests/test_compartments_saddle.py @@ -419,3 +419,34 @@ def test_saddle(request, tmpdir): # TODO: tests after adding input agreement, e.g. # asserting saddle.saddle(clr, cis-type-expected, track, "trans") # throws an error + + +def test_saddle_strength(): + # asymmetric bin counts on the diagonal (C[0, 0] != C[-1, -1]) so that + # a naive extent=0 default of 0 would be distinguishable from the + # correct value + S = np.array( + [ + [10.0, 2.0, 1.0], + [2.0, 5.0, 2.0], + [1.0, 2.0, 8.0], + ] + ) + C = np.array( + [ + [4.0, 2.0, 1.0], + [2.0, 5.0, 2.0], + [1.0, 2.0, 2.0], + ] + ) + + ratios = saddle.saddle_strength(S, C) + + # extent=0: ((S/C)[0, 0] + (S/C)[-1, -1]) / (2 * (S/C)[-1, 0]) + # hand-computed from S, C above + expected = np.array([3.0, 3.0, 1.5]) + + assert np.allclose(ratios, expected) + # regression check: extent=0 must not be left at its zero-initialized + # default (see https://github.com/open2c/cooltools/issues/565) + assert ratios[0] != 0