Two convention bugs in `Result.mueller_matrix` and `Experiment.set_vector` (v0.22.6)

[antorbus]

Two convention bugs in Result.mueller_matrix and Experiment.set_vector (v0.22.6)

Both verified on pyElli 0.22.6 / Python 3.14 / NumPy. Both have one-line fixes.
Neither affects pyElli's psi, delta, R, T, or jones_matrix_* outputs (those are derived from rho directly), but each silently corrupts downstream code that uses the affected functions as transformations.

---

Bug 1: Result.mueller_matrix returns the wrong Kronecker order in elli/result.py:228.

s_kron_s_star = np.einsum(
    "aij,akl->aikjl", np.conjugate(self.rho_matrix), self.rho_matrix
).reshape([self.rho_matrix.shape[0], 4, 4])

This is kron(conj(J), J) = J* ⊗ J (the column-major Kronecker), but the surrounding reshape is C-order (row-major). The conventions don't match. The defining property of the Mueller matrix, S_out = M @ S_in, fails.

The standard derivation (Goldstein 2003 Eq. 4.40; Tompkins–Irene 2005 §1.5; Wikipedia "Mueller calculus"; Cloude 1986; Aiello–Woerdman arXiv:1906.11198) uses the row-major vec–Kronecker identity vec_row(ABC) = (A ⊗ Cᵀ) vec_row(B), applied to Φ_out = J Φ Jᴴ with (Jᴴ)ᵀ = J*, giving vec_row(Φ_out) = (J ⊗ J*) vec_row(Φ). The correct Kronecker is kron(J, conj(J)), operands swapped relative to the current code.

Reproduce (pure NumPy, no pyElli sample needed):

import numpy as np

A = np.array([[1, 0, 0, 1], [1, 0, 0, -1], [0, 1, 1, 0], [0, 1j, -1j, 0]], dtype=complex)
J = np.array([[0.5, 0.1-0.2j], [0.05+0.3j, 0.3+0.2j]])

# input LCP per pyElli's S_3 = -2 Im(E_x E_y*)  ->  Jones = (1, +i)/sqrt(2)
E_in = np.array([1, 1j]) / np.sqrt(2)
S_in = np.array([1.0, 0.0, 0.0, 1.0])

def stokes(e):  # pyElli's forward formula, experiment.py:92-107
    return np.array([abs(e[0])**2 + abs(e[1])**2,
                     abs(e[0])**2 - abs(e[1])**2,
                     2 * np.real(e[0] * np.conj(e[1])),
                     -2 * np.imag(e[0] * np.conj(e[1]))])

S_direct  = stokes(J @ E_in)
M_pyelli  = (A @ np.kron(np.conj(J), J) @ np.linalg.inv(A)).real     # current code
M_correct = (A @ np.kron(J, np.conj(J)) @ np.linalg.inv(A)).real     # standard

print("S_direct      =", S_direct)
print("M_pyelli @ S  =", M_pyelli @ S_in)        # MISMATCH
print("M_correct @ S =", M_correct @ S_in)       # matches direct to ~1e-15

Output:

S_direct      = [0.44125  0.05875 -0.045    0.435  ]
M_pyelli @ S  = [0.08125  0.01875  0.075   -0.025  ]   # off by 0.36 on S_0
M_correct @ S = [0.44125  0.05875 -0.045    0.435  ]

100/100 random complex Jones matrices fail with the current code at tol 1e-10; 0/100 fail with the corrected order.

For an isotropic sample (diagonal J = diag(r_p, r_s)), the discrepancy collapses to a sign flip on M[2,3] and M[3,2] — equivalently, the returned Mueller matches Azzam–Bashara if and only if the user later flips the sign of Δ.

Fix (one line, elli/result.py:228):

# was:
s_kron_s_star = np.einsum("aij,akl->aikjl", np.conjugate(self.rho_matrix), self.rho_matrix).reshape(...)
# fix:
s_kron_s_star = np.einsum("aij,akl->aikjl", self.rho_matrix, np.conjugate(self.rho_matrix)).reshape(...)

After the fix, Result.mueller_matrix satisfies M @ S_in == Stokes(J @ Jones(S_in)) to machine precision and matches the textbook isotropic-block form with the standard Δ sign.

---

Bug 2: Experiment.set_vector(Stokes → Jones) inverts S_3 on round-trip in elli/experiment.py:129.

b = U / (2 * a) - 1j * V / (2 * a)

The minus sign on V is wrong under pyElli's own forward Stokes convention (experiment.py:92-107, S_3 = -2 Im(E_x E_y*)). Round-tripping any fully-polarised input with V ≠ 0 flips S_3.

Reproduce.

import numpy as np
import elli

exp = elli.Experiment(
    elli.Structure(front=elli.AIR, layers=[], back=elli.AIR),
    [550.0], 0.0,
    vector=np.array([1.0, 0.0, 0.0, 1.0]),    # LCP under pyElli's S_3 convention
)
print("input    .stokes_vector  =", exp.stokes_vector)     # [1. 0. 0. 1.]
print("from     .jones_vector   =", exp.jones_vector)      # (0.7071, -0.7071j)

j = exp.jones_vector
S_recovered = np.array([abs(j[0])**2 + abs(j[1])**2,
                        abs(j[0])**2 - abs(j[1])**2,
                        2 * np.real(j[0] * np.conj(j[1])),
                        -2 * np.imag(j[0] * np.conj(j[1]))])
print("re-stokes from jones    =", S_recovered)            # [1. 0. 0. -1.]   <-- S_3 flipped

Linear states (V = 0) round-trip cleanly because the buggy -1j V / (2a) term is zero.

Fix (one character, elli/experiment.py:129):

# was:
b = U / (2 * a) - 1j * V / (2 * a)
# fix:
b = U / (2 * a) + 1j * V / (2 * a)

After the fix, the round-trip identity Stokes(jones_from_set_vector(S)) == S holds exactly for every fully-polarised S.

---

Notes

• We discovered both while building a torch port of Solver4x4 and applying Result.mueller_matrix as a Stokes-domain transformation.

[MarJMue]

Hi @antorbus,

thank you very much for pointing out these mistakes and the very detailed report.
I will run a few checks myself and publish a fix as soon as possible.

Please let us know if there are any news about the Pytorch Solver4x4. That would be something we would like to try out as well.