diff options
Diffstat (limited to 'worldalign/field_anatomy.py')
| -rw-r--r-- | worldalign/field_anatomy.py | 29 |
1 files changed, 22 insertions, 7 deletions
diff --git a/worldalign/field_anatomy.py b/worldalign/field_anatomy.py index 6156018..9e04ec6 100644 --- a/worldalign/field_anatomy.py +++ b/worldalign/field_anatomy.py @@ -49,23 +49,38 @@ def correlate(first: np.ndarray, second: np.ndarray) -> float: return float(np.corrcoef(offdiagonal(first), offdiagonal(second))[0, 1]) -def degree_part(matrix: np.ndarray) -> tuple[np.ndarray, np.ndarray]: +def degree_part( + matrix: np.ndarray, sweeps: int = 200 +) -> tuple[np.ndarray, np.ndarray]: """Split into the additive row/column-mean model and its residual. The two-way additive fit m + r_i + c_j is what a pure busyness effect produces: every entry explained by how prominent each of its two scenes is, with nothing said about the pair. + + The fit is swept to convergence rather than taken in one pass. With the + diagonal excluded the design is unbalanced, so a single pass of row and + column means leaves an O(1/n) piece of a pure degree effect behind -- + small, but it is exactly the quantity this function exists to remove. + Alternating the two sweeps converges to the correct fit. """ size = len(matrix) mask = ~np.eye(size, dtype=bool) - work = matrix.copy().astype(np.float64) - np.fill_diagonal(work, np.nan) + work = np.where(mask, matrix.astype(np.float64), np.nan) grand = np.nanmean(work) - rows = np.nanmean(work, axis=1, keepdims=True) - grand - cols = np.nanmean(work, axis=0, keepdims=True) - grand + rows = np.zeros((size, 1)) + cols = np.zeros((1, size)) + for _ in range(sweeps): + residual = work - (grand + rows + cols) + step = np.nanmean(residual, axis=1, keepdims=True) + rows = rows + step + residual = work - (grand + rows + cols) + step_columns = np.nanmean(residual, axis=0, keepdims=True) + cols = cols + step_columns + if max(np.abs(step).max(), np.abs(step_columns).max()) < 1e-12: + break fitted = grand + rows + cols - residual = np.where(mask, work - fitted, 0.0) - return np.where(mask, fitted, 0.0), residual + return np.where(mask, fitted, 0.0), np.where(mask, work - fitted, 0.0) def spectral_strip(matrix: np.ndarray, components: int) -> np.ndarray: |