⚡ Bolt: Vectorize Newton-Raphson M-step in MMLE estimator - #383
⚡ Bolt: Vectorize Newton-Raphson M-step in MMLE estimator#383seonghobae wants to merge 2 commits into
Conversation
|
👋 Jules, reporting for duty! I'm here to lend a hand with this pull request. When you start a review, I'll add a 👀 emoji to each comment to let you know I've read it. I'll focus on feedback directed at me and will do my best to stay out of conversations between you and other bots or reviewers to keep the noise down. I'll push a commit with your requested changes shortly after. Please note there might be a delay between these steps, but rest assured I'm on the job! For more direct control, you can switch me to Reactive Mode. When this mode is on, I will only act on comments where you specifically mention me with New to Jules? Learn more at jules.google/docs. For security, I will only act on instructions from the user who triggered this task. |
📝 WalkthroughWalkthroughMMLE M-step의 항목별 Newton-Raphson 갱신을 활성 항목 배열 기반의 벡터화 반복으로 변경했습니다. Hessian이 유효하지 않거나 수렴한 항목은 후속 반복에서 제외합니다. ChangesMMLE Newton 반복 벡터화
Estimated code review effort: 2 (Simple) | ~15 minutes Possibly related PRs
🚥 Pre-merge checks | ✅ 5✅ Passed checks (5 passed)
✨ Finishing Touches📝 Generate docstrings
🧪 Generate unit tests (beta)
Comment |
There was a problem hiding this comment.
Actionable comments posted: 2
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Inline comments:
In `@python/fast_mlsirm/estimators/mmle.py`:
- Around line 135-142: Update the computation around eta, p, w, and resid to
process active item indices in fixed-size batches rather than allocating arrays
for all active items at once. Preserve the existing vectorized operations within
each batch and write each batch’s done results back to the corresponding
original item indices; alternatively, enforce an explicit upper bound on
user-derived work dimensions before allocation.
- Around line 126-167: Update the Newton vectorization guidance in
.jules/bolt.md at lines 44-46 to conditionally allow vectorizing the independent
item updates in the ULS2PLM 2PL M-step, as implemented in mmle.py lines 126-167.
State that the model equations, observed mask, parameter bounds, and termination
criteria must remain unchanged; no direct code change is required at the mmle.py
site.
🪄 Autofix (Beta)
Fix all unresolved CodeRabbit comments on this PR:
- Push a commit to this branch (recommended)
- Create a new PR with the fixes
ℹ️ Review info
⚙️ Run configuration
Configuration used: Organization UI
Review profile: CHILL
Plan: Pro Plus
Run ID: d60d7a2f-258c-47d3-9c10-763f8d8d8634
⛔ Files ignored due to path filters (1)
uv.lockis excluded by!**/*.lock
📒 Files selected for processing (2)
.jules/bolt.mdpython/fast_mlsirm/estimators/mmle.py
| # Optimized M-step: Vectorize Newton-Raphson across all active items simultaneously | ||
| # to avoid slow Python loops and per-item array overhead. | ||
| active = np.ones(n_items, dtype=bool) | ||
| nodes_sq = nodes * nodes | ||
|
|
||
| for _ in range(25): | ||
| if not active.any(): | ||
| break | ||
|
|
||
| eta = a_new[active, None] * nodes[None, :] + b_new[active, None] | ||
| p = _sigmoid(eta) | ||
|
|
||
| n_act = n_iq[active] | ||
| r_act = r_iq[active] | ||
|
|
||
| w = n_act * p * (1.0 - p) | ||
| resid = r_act - n_act * p | ||
|
|
||
| g_a = resid @ nodes - ridge_a * a_new[active] | ||
| g_b = resid.sum(axis=1) - ridge_b * b_new[active] | ||
|
|
||
| h_aa = -(w @ nodes_sq) - ridge_a | ||
| h_bb = -w.sum(axis=1) - ridge_b | ||
| h_ab = -(w @ nodes) | ||
|
|
||
| det = h_aa * h_bb - h_ab * h_ab | ||
| det_ok = np.abs(det) >= 1e-12 | ||
|
|
||
| da = np.zeros_like(g_a) | ||
| db = np.zeros_like(g_b) | ||
|
|
||
| da[det_ok] = (h_bb[det_ok] * g_a[det_ok] - h_ab[det_ok] * g_b[det_ok]) / det[det_ok] | ||
| db[det_ok] = (h_aa[det_ok] * g_b[det_ok] - h_ab[det_ok] * g_a[det_ok]) / det[det_ok] | ||
|
|
||
| a_new[active] -= da | ||
| b_new[active] -= db | ||
| a_new[active] = np.clip(a_new[active], 1e-3, 10.0) | ||
|
|
||
| done = (~det_ok) | ((np.abs(da) + np.abs(db)) < 1e-8) | ||
| if done.any(): | ||
| active_idx = np.where(active)[0] | ||
| active[active_idx[done]] = False |
There was a problem hiding this comment.
📐 Maintainability & Code Quality | 🟡 Minor | ⚡ Quick win
🧩 Analysis chain
🏁 Script executed:
command -v codegraph || true
codegraph --help 2>&1 | head -80 || trueRepository: ContextualWisdomLab/fast-mlsirm
Length of output: 212
🏁 Script executed:
#!/bin/bash
set -e
printf '%s\n' '--- repository guidance ---'
for f in AGENTS.md; do
if [ -f "$f" ]; then cat -n "$f"; fi
done
printf '%s\n' '--- changed-file summary ---'
git diff --stat -- python/fast_mlsirm/estimators/mmle.py .jules/bolt.md
printf '%s\n' '--- estimator diff ---'
git diff --unified=80 -- python/fast_mlsirm/estimators/mmle.py
printf '%s\n' '--- bolt guidance ---'
cat -n .jules/bolt.md | sed -n '35,55p'
printf '%s\n' '--- related files ---'
git ls-files | grep -E '(^|/)(test.*mmle|mmle|MLSIRM|MLS2PLM|AGENTS|.*fit.*stat)' || true
printf '%s\n' '--- relevant symbols and references ---'
rg -n -C 3 'Newton|observed|n_iq|r_iq|ridge_a|ridge_b|a_new|b_new|MLS2PLM|MLSIRM|fit statistic' \
python tests docs .jules 2>/dev/null | head -400 || trueRepository: ContextualWisdomLab/fast-mlsirm
Length of output: 44683
🏁 Script executed:
#!/bin/bash
set -e
printf '%s\n' '--- Python MMLE implementation ---'
cat -n python/fast_mlsirm/estimators/mmle.py | sed -n '1,240p'
printf '%s\n' '--- MMLE estimator tests ---'
cat -n tests/test_estimator_mmle.py | sed -n '1,320p'
printf '%s\n' '--- Rust MMLE implementation ---'
cat -n crates/mlsirm-core/src/mmle.rs | sed -n '1,300p'
printf '%s\n' '--- MMLE design document ---'
cat -n docs/mmle_marginal_lsirm_design.md | sed -n '1,135p'
printf '%s\n' '--- formula compilation ---'
cat -n docs/papers/mmle-lsirm-formula-compilation.md | sed -n '1,180p'
printf '%s\n' '--- recent commits ---'
git log --oneline -8 -- python/fast_mlsirm/estimators/mmle.py .jules/bolt.mdRepository: ContextualWisdomLab/fast-mlsirm
Length of output: 43361
🏁 Script executed:
#!/bin/bash
set -e
printf '%s\n' '--- vectorization commit ---'
git show --stat --oneline c103514
git show --format=fuller --no-ext-diff --unified=80 c103514 -- \
python/fast_mlsirm/estimators/mmle.py .jules/bolt.md tests/test_estimator_mmle.py
printf '%s\n' '--- parent implementation ---'
git show c103514^:python/fast_mlsirm/estimators/mmle.py | \
cat -n | sed -n '105,190p'
printf '%s\n' '--- commit files ---'
git diff-tree --no-commit-id --name-status -r c103514Repository: ContextualWisdomLab/fast-mlsirm
Length of output: 18528
🏁 Script executed:
#!/bin/bash
python3 - <<'PY'
import numpy as np
from numpy.polynomial.hermite_e import hermegauss
def sigmoid(x):
return 1.0 / (1.0 + np.exp(-np.clip(x, -35.0, 35.0)))
def fit(y, observed, vectorized, n_nodes=11, max_iter=40, seed=1):
y = np.asarray(y, dtype=np.float64)
observed = np.asarray(observed, dtype=bool)
y_filled = np.where(observed, y, 0.0)
obs_f = observed.astype(np.float64)
nodes, raw = hermegauss(n_nodes)
weights = raw / raw.sum()
rng = np.random.default_rng(seed)
p_item = (y_filled * obs_f).sum(0) / np.clip(obs_f.sum(0), 1.0, None)
p_item = np.clip(p_item, 0.02, 0.98)
a = np.ones(y.shape[1]) + 0.01 * rng.standard_normal(y.shape[1])
b = np.log(p_item / (1.0 - p_item))
trace = []
for _ in range(max_iter):
logit = nodes[:, None] * a[None, :] + b[None, :]
log_p1 = -np.logaddexp(0.0, -logit)
log_p0 = -np.logaddexp(0.0, logit)
pos = (y_filled * obs_f) @ log_p1.T
neg = ((1.0 - y_filled) * obs_f) @ log_p0.T
joint = pos + neg + np.log(weights)[None, :]
max_j = joint.max(axis=1, keepdims=True)
stab = np.exp(joint - max_j)
denom = stab.sum(axis=1, keepdims=True)
posterior = stab / denom
trace.append(float((max_j[:, 0] + np.log(denom[:, 0])).sum()))
n_iq = obs_f.T @ posterior
r_iq = (obs_f * y_filled).T @ posterior
a_new, b_new = a.copy(), b.copy()
if vectorized:
active = np.ones(len(a), dtype=bool)
nodes_sq = nodes * nodes
for _ in range(25):
if not active.any():
break
eta = a_new[active, None] * nodes[None, :] + b_new[active, None]
p = sigmoid(eta)
n_act, r_act = n_iq[active], r_iq[active]
w = n_act * p * (1.0 - p)
resid = r_act - n_act * p
ga = resid @ nodes - 1e-3 * a_new[active]
gb = resid.sum(axis=1) - 1e-3 * b_new[active]
haa = -(w @ nodes_sq) - 1e-3
hbb = -w.sum(axis=1) - 1e-3
hab = -(w @ nodes)
det = haa * hbb - hab * hab
ok = np.abs(det) >= 1e-12
da = np.zeros_like(ga)
db = np.zeros_like(gb)
da[ok] = (hbb[ok] * ga[ok] - hab[ok] * gb[ok]) / det[ok]
db[ok] = (haa[ok] * gb[ok] - hab[ok] * ga[ok]) / det[ok]
a_new[active] -= da
b_new[active] -= db
a_new[active] = np.clip(a_new[active], 1e-3, 10.0)
done = (~ok) | ((np.abs(da) + np.abs(db)) < 1e-8)
if done.any():
idx = np.where(active)[0]
active[idx[done]] = False
else:
for i in range(len(a)):
ai, bi = a[i], b[i]
for _ in range(25):
eta = ai * nodes + bi
p = sigmoid(eta)
w = n_iq[i] * p * (1.0 - p)
resid = r_iq[i] - n_iq[i] * p
ga = float((resid * nodes).sum()) - 1e-3 * ai
gb = float(resid.sum()) - 1e-3 * bi
haa = -float((w * nodes * nodes).sum()) - 1e-3
hbb = -float(w.sum()) - 1e-3
hab = -float((w * nodes).sum())
det = haa * hbb - hab * hab
if abs(det) < 1e-12:
break
da = (hbb * ga - hab * gb) / det
db = (haa * gb - hab * ga) / det
ai -= da
bi -= db
ai = float(np.clip(ai, 1e-3, 10.0))
if abs(da) + abs(db) < 1e-8:
break
a_new[i], b_new[i] = ai, bi
a, b = a_new, b_new
return a, b, np.asarray(trace)
rng = np.random.default_rng(202608)
worst = np.zeros(3)
cases = 0
for n_persons, n_items, n_nodes in [(1, 1, 3), (20, 7, 7), (80, 25, 11), (120, 40, 41)]:
for case in range(20):
y = rng.integers(0, 2, size=(n_persons, n_items)).astype(float)
observed = rng.random(y.shape) > (0.95 if case == 0 else 0.35)
if not observed.any():
observed[0, 0] = True
y[~observed] = np.nan
old = fit(y, observed, False, n_nodes=n_nodes)
new = fit(y, observed, True, n_nodes=n_nodes)
diffs = np.array([
np.max(np.abs(old[0] - new[0])),
np.max(np.abs(old[1] - new[1])),
np.max(np.abs(old[2] - new[2])),
])
worst = np.maximum(worst, diffs)
cases += 1
print(f"cases={cases}")
print(f"max_abs_diff_a={worst[0]:.3e}")
print(f"max_abs_diff_b={worst[1]:.3e}")
print(f"max_abs_diff_loglik={worst[2]:.3e}")
print("equivalent_within_1e-12=", bool(np.all(worst < 1e-12)))
PYRepository: ContextualWisdomLab/fast-mlsirm
Length of output: 287
🏁 Script executed:
#!/bin/bash
set -e
printf '%s\n' '--- paper/design history ---'
git log --oneline --all -- docs/mmle_marginal_lsirm_design.md docs/papers/mmle-lsirm-formula-compilation.md | head -20
printf '%s\n' '--- design document at vectorization parent ---'
git cat-file -e c103514^:docs/mmle_marginal_lsirm_design.md && \
git show c103514^:docs/mmle_marginal_lsirm_design.md | sed -n '1,20p'
printf '%s\n' '--- pure-Python Newton equivalence probe ---'
python3 - <<'PY'
import math
import random
def sigmoid(x):
x = max(-35.0, min(35.0, x))
return 1.0 / (1.0 + math.exp(-x))
def scalar(n_iq, r_iq, nodes, a, b):
aa, bb = list(a), list(b)
for i in range(len(aa)):
ai, bi = aa[i], bb[i]
for _ in range(25):
ga = gb = haa = hbb = hab = 0.0
for q, node in enumerate(nodes):
p = sigmoid(ai * node + bi)
n = n_iq[i][q]
w = n * p * (1.0 - p)
resid = r_iq[i][q] - n * p
ga += resid * node
gb += resid
haa -= w * node * node
hbb -= w
hab -= w * node
ga -= 1e-3 * ai
gb -= 1e-3 * bi
haa -= 1e-3
hbb -= 1e-3
det = haa * hbb - hab * hab
if abs(det) < 1e-12:
break
da = (hbb * ga - hab * gb) / det
db = (haa * gb - hab * ga) / det
ai -= da
bi -= db
ai = max(1e-3, min(10.0, ai))
if abs(da) + abs(db) < 1e-8:
break
aa[i], bb[i] = ai, bi
return aa, bb
def batched(n_iq, r_iq, nodes, a, b):
aa, bb = list(a), list(b)
active = [True] * len(aa)
nodes_sq = [x * x for x in nodes]
for _ in range(25):
if not any(active):
break
active_indices = [i for i, flag in enumerate(active) if flag]
updates = []
for i in active_indices:
eta = [aa[i] * node + bb[i] for node in nodes]
p = [sigmoid(x) for x in eta]
w = [n * pp * (1.0 - pp) for n, pp in zip(n_iq[i], p)]
resid = [r - n * pp for r, n, pp in zip(r_iq[i], n_iq[i], p)]
ga = sum(x * node for x, node in zip(resid, nodes)) - 1e-3 * aa[i]
gb = sum(resid) - 1e-3 * bb[i]
haa = -sum(x * node_sq for x, node_sq in zip(w, nodes_sq)) - 1e-3
hbb = -sum(w) - 1e-3
hab = -sum(x * node for x, node in zip(w, nodes))
det = haa * hbb - hab * hab
ok = abs(det) >= 1e-12
if ok:
da = (hbb * ga - hab * gb) / det
db = (haa * gb - hab * ga) / det
else:
da = db = 0.0
updates.append((i, da, db, ok))
for i, da, db, ok in updates:
aa[i] -= da
bb[i] -= db
aa[i] = max(1e-3, min(10.0, aa[i]))
if (not ok) or abs(da) + abs(db) < 1e-8:
active[i] = False
return aa, bb
rng = random.Random(202608)
worst_a = worst_b = 0.0
cases = 0
for n_items, n_nodes in [(1, 3), (7, 7), (25, 11), (40, 41)]:
for case in range(100):
nodes = [rng.uniform(-12.0, 12.0) for _ in range(n_nodes)]
n_iq, r_iq = [], []
for _ in range(n_items):
counts = [0.0 if case == 0 else rng.uniform(0.0, 100.0)
for _ in range(n_nodes)]
responses = [n * rng.random() for n in counts]
n_iq.append(counts)
r_iq.append(responses)
a = [rng.uniform(1e-3, 10.0) for _ in range(n_items)]
b = [rng.uniform(-20.0, 20.0) for _ in range(n_items)]
sa, sb = scalar(n_iq, r_iq, nodes, a, b)
va, vb = batched(n_iq, r_iq, nodes, a, b)
worst_a = max(worst_a, *(abs(x-y) for x, y in zip(sa, va)))
worst_b = max(worst_b, *(abs(x-y) for x, y in zip(sb, vb)))
cases += 1
print(f"cases={cases}")
print(f"max_abs_diff_a={worst_a:.3e}")
print(f"max_abs_diff_b={worst_b:.3e}")
print("same_within_1e-12=", worst_a < 1e-12 and worst_b < 1e-12)
PYRepository: ContextualWisdomLab/fast-mlsirm
Length of output: 1568
🏁 Script executed:
#!/bin/bash
set -e
printf '%s\n' '--- Rust parity test ---'
if [ -f tests/test_rust_parity.py ]; then
cat -n tests/test_rust_parity.py | sed -n '1,320p'
else
printf '%s\n' 'tests/test_rust_parity.py is absent'
fi
printf '%s\n' '--- MMLE parity references ---'
rg -n -C 4 'mmle|MMLE|fit_mmle_2pl|rust.*numpy|numpy.*rust|parity' \
tests python crates | head -500 || true
printf '%s\n' '--- vectorization commit tests ---'
git diff c103514^ c103514 --name-only -- tests
printf '%s\n' '--- available Python runtimes/dependencies ---'
python3 - <<'PY'
try:
import numpy
print("numpy_available", numpy.__version__)
except Exception as exc:
print("numpy_unavailable", type(exc).__name__)
PYRepository: ContextualWisdomLab/fast-mlsirm
Length of output: 46038
.jules/bolt.md의 Newton 벡터화 지침을 조건부로 수정하십시오.
python/fast_mlsirm/estimators/mmle.py 변경은 기존 ULS2PLM 2PL M-step의 독립적인 item update만 벡터화합니다. 모델 공식, observed 마스크, 범위 제한, 종료 조건은 유지됩니다. 문헌 근거는 docs/mmle_marginal_lsirm_design.md와 docs/papers/mmle-lsirm-formula-compilation.md에 있습니다. 관련 테스트는 tests/test_estimator_mmle.py와 tests/test_rust_parity.py에 있습니다.
권장 문서 변경
- **Action:** Replace outer Python loops over independent elements with a fully vectorized inner loop.
+ **Action:** Before vectorizing an estimator update, document the paper-supported scope and verify that the formula, observed-mask semantics, bounds, stopping behavior, and Rust–NumPy parity are unchanged. If these conditions hold, replace outer Python loops over independent elements with a fully vectorized inner loop.📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
| # Optimized M-step: Vectorize Newton-Raphson across all active items simultaneously | |
| # to avoid slow Python loops and per-item array overhead. | |
| active = np.ones(n_items, dtype=bool) | |
| nodes_sq = nodes * nodes | |
| for _ in range(25): | |
| if not active.any(): | |
| break | |
| eta = a_new[active, None] * nodes[None, :] + b_new[active, None] | |
| p = _sigmoid(eta) | |
| n_act = n_iq[active] | |
| r_act = r_iq[active] | |
| w = n_act * p * (1.0 - p) | |
| resid = r_act - n_act * p | |
| g_a = resid @ nodes - ridge_a * a_new[active] | |
| g_b = resid.sum(axis=1) - ridge_b * b_new[active] | |
| h_aa = -(w @ nodes_sq) - ridge_a | |
| h_bb = -w.sum(axis=1) - ridge_b | |
| h_ab = -(w @ nodes) | |
| det = h_aa * h_bb - h_ab * h_ab | |
| det_ok = np.abs(det) >= 1e-12 | |
| da = np.zeros_like(g_a) | |
| db = np.zeros_like(g_b) | |
| da[det_ok] = (h_bb[det_ok] * g_a[det_ok] - h_ab[det_ok] * g_b[det_ok]) / det[det_ok] | |
| db[det_ok] = (h_aa[det_ok] * g_b[det_ok] - h_ab[det_ok] * g_a[det_ok]) / det[det_ok] | |
| a_new[active] -= da | |
| b_new[active] -= db | |
| a_new[active] = np.clip(a_new[active], 1e-3, 10.0) | |
| done = (~det_ok) | ((np.abs(da) + np.abs(db)) < 1e-8) | |
| if done.any(): | |
| active_idx = np.where(active)[0] | |
| active[active_idx[done]] = False | |
| ## 2026-08-02 - Vectorize iterative Newton steps across independent dimensions | |
| **Learning:** During the M-step of MMLE estimation, performing a Python loop over `n_items` to run an iterative inner loop (like Newton-Raphson) creates severe overhead. Independent scalar iterative calculations (like updating item difficulties) are very slow in Python. | |
| **Action:** Before vectorizing an estimator update, document the paper-supported scope and verify that the formula, observed-mask semantics, bounds, stopping behavior, and Rust–NumPy parity are unchanged. If these conditions hold, replace outer Python loops over independent elements with a fully vectorized inner loop. Use a boolean array to track active components (e.g. `active = np.ones(n_items, dtype=bool)`) and process only the active subset in each iteration, stopping early for components that converge via boolean masking. This transforms slow Python scalar iterations into fast C-level numpy operations. |
📍 Affects 2 files
python/fast_mlsirm/estimators/mmle.py#L126-L167(this comment).jules/bolt.md#L44-L46
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In `@python/fast_mlsirm/estimators/mmle.py` around lines 126 - 167, Update the
Newton vectorization guidance in .jules/bolt.md at lines 44-46 to conditionally
allow vectorizing the independent item updates in the ULS2PLM 2PL M-step, as
implemented in mmle.py lines 126-167. State that the model equations, observed
mask, parameter bounds, and termination criteria must remain unchanged; no
direct code change is required at the mmle.py site.
Source: Coding guidelines
| eta = a_new[active, None] * nodes[None, :] + b_new[active, None] | ||
| p = _sigmoid(eta) | ||
|
|
||
| n_act = n_iq[active] | ||
| r_act = r_iq[active] | ||
|
|
||
| w = n_act * p * (1.0 - p) | ||
| resid = r_act - n_act * p |
There was a problem hiding this comment.
🩺 Stability & Availability | 🟠 Major | 🏗️ Heavy lift
활성 항목을 배치 처리하거나 작업 배열 크기를 제한하십시오.
라인 135-142는 eta, p, w, resid를 모두 (n_active, n_nodes) 크기로 생성합니다. n_items와 n_nodes에는 상한이 없습니다. 많은 항목을 가진 입력에서는 이 배열들이 동시에 큰 피크 메모리를 사용하고 프로세스를 OOM으로 종료할 수 있습니다.
활성 항목 인덱스를 고정 크기 배치로 나누어 처리하십시오. 각 배치에서는 현재 벡터화 연산을 유지하고, done 결과를 원래 항목 인덱스에 다시 반영하십시오. 또는 할당 전에 명시적인 작업 요소 상한을 검증하십시오.
As per coding guidelines, “bound user-derived array dimensions” 규칙을 적용했습니다.
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In `@python/fast_mlsirm/estimators/mmle.py` around lines 135 - 142, Update the
computation around eta, p, w, and resid to process active item indices in
fixed-size batches rather than allocating arrays for all active items at once.
Preserve the existing vectorized operations within each batch and write each
batch’s done results back to the corresponding original item indices;
alternatively, enforce an explicit upper bound on user-derived work dimensions
before allocation.
Source: Coding guidelines
Closed during maintainer review. This change moves additional MMLE Newton-Raphson computation into NumPy/Python, which conflicts with the repository requirement that numerical computation execute through the Rust backend. A replacement must implement the M-step in Rust, provide CPU multithreading and GPU execution paths, preserve the cited model equations and termination criteria, include parameter-recovery validation, and satisfy 100% docstring/coverage gates. The unresolved batching review also remains valid.