mlem

anri.index.mlem(d, cand, pred, ring_j, pos, scan, dims, f0, n_iter, vb, log=<built-in function print>, censor=0.0, accel=False)

Fit occupancies by MLEM: f <- f A^T(d / A f) / A^T 1, n_iter times, logging the Poisson deviance.

Empty bins are censored above censor (counts per bin, censored_ratio()); 0 for plain MLEM.

accel: SQUAREM (Varadhan & Roland 2008, Scand. J. Stat. 35, 335). Two MLEM steps f1, f2 from f; then f - 2 a r + a^2 v with r = f1 - f, v = f2 - 2 f1 + f and a = -norm(r) / norm(v) (at most -1: never shorter than plain MLEM), kept above 1e-3 f2 (MLEM multiplies: an occupancy at 0 could never come back), and one more MLEM step. If that fits worse than f1, f2’s step is kept instead, so the deviance never rises. Thin features (twins one or two voxels thick) converge in far fewer steps. n_iter counts MLEM steps either way (3 per SQUAREM step), so the cost is the same.

Parameters:

• d (Array) – [n_cells] histogram

• cand (Array) – See forward()

• pred (tuple) – See forward()

• ring_j (Array) – See forward()

• pos (Array) – See forward()

• scan (dict) – See forward()

• dims (tuple) – See forward()

• vb (int) – See forward()

• f0 (Array) – [Nv, K] starting occupancies

• n_iter (int) – Iterations

• log (Callable, default: <built-in function print>) – Progress messages

Returns:

jax.Array – [Nv, K] occupancies