fix eme r2l lossy projection

src/meow/eme/common.py

@@ -110,8 +110,10 @@ def compute_interface_s_matrix(
 
     # enforce S@S.H is diagonal: HACK!
     if enforce_lossy_unitarity:
-        U, s, V = np.linalg.svd(S)
-        S = np.diag(s) @ U @ V
+        # Project to a contractive matrix by clipping singular values to <= 1.
+        U, s, Vh = np.linalg.svd(S, full_matrices=False)
+        s_clipped = np.minimum(s, 1.0)
+        S = U @ np.diag(s_clipped) @ Vh
 
     # ensure reciprocity: HACK?
     if enforce_reciprocity:

src/meow/eme/propagate.py

@@ -81,7 +81,8 @@ def r2l_matrices(
     """Right to left S-matrices."""
     Ss = [pairs[-1]]
 
-    for p in pairs[-1::-1]:
+    # Already seeded with the rightmost pair; only prepend remaining pairs.
+    for p in pairs[-2::-1]:
         Ss.append(_connect_two(p, Ss[-1], sax_backend))
 
     return Ss[::-1]

src/tests/test_eme_fixes.py

@@ -0,0 +1,71 @@
+from typing import Any, cast
+
+import numpy as np
+import pytest
+import sax
+
+import meow.eme.common as eme_common
+import meow.eme.propagate as eme_propagate
+from meow.mode import Mode
+
+
+def test_r2l_matrices_does_not_duplicate_last_pair(
+    monkeypatch: pytest.MonkeyPatch,
+) -> None:
+    calls: list[tuple[str, str]] = []
+
+    def fake_connect_two(l: Any, r: Any, sax_backend: sax.Backend) -> str:  # noqa: ARG001
+        calls.append((l, r))
+        return f"({l}>{r})"
+
+    monkeypatch.setattr(eme_propagate, "_connect_two", fake_connect_two)
+
+    pairs = cast(list[sax.STypeMM], ["p0", "p1", "p2"])
+    matrices = eme_propagate.r2l_matrices(pairs, sax_backend="klu")
+
+    assert calls == [("p1", "p2"), ("p0", "(p1>p2)")]
+    assert matrices == ["(p0>(p1>p2))", "(p1>p2)", "p2"]
+
+
+def test_enforce_lossy_unitarity_projects_to_contractive_matrix(
+    monkeypatch: pytest.MonkeyPatch,
+) -> None:
+    left = cast(Mode, object())
+    right = cast(Mode, object())
+
+    def fake_inner_product_conj(a: Any, b: Any) -> float:
+        if a is left and b is left:
+            return 1.0
+        if a is right and b is right:
+            return 1.0
+        if a is left and b is right:
+            return 0.01
+        if a is right and b is left:
+            return 10.0
+        msg = "unexpected mode pair"
+        raise AssertionError(msg)
+
+    monkeypatch.setattr(eme_common, "inner_product_conj", fake_inner_product_conj)
+
+    S_no, _ = eme_common.compute_interface_s_matrix(
+        [left],
+        [right],
+        enforce_lossy_unitarity=False,
+        ignore_warnings=False,
+    )
+    S_yes, _ = eme_common.compute_interface_s_matrix(
+        [left],
+        [right],
+        enforce_lossy_unitarity=True,
+        ignore_warnings=False,
+    )
+
+    s_no = np.linalg.svd(S_no, compute_uv=False)
+    assert float(s_no.max()) > 1.0
+
+    U, s, Vh = np.linalg.svd(S_no, full_matrices=False)
+    expected = U @ np.diag(np.minimum(s, 1.0)) @ Vh
+    assert np.allclose(S_yes, expected)
+
+    s_yes = np.linalg.svd(S_yes, compute_uv=False)
+    assert float(s_yes.max()) <= 1.0 + 1e-12