Skip to content

Commit 8825607

Browse files
fix(space): correct centre offset and unbias SphericalShell sampling (#38)
## Summary Two related bugs in `maws/space.py` that silently distort sampling for `Sphere` and `SphericalShell`: - **Centre offset ignored in `generator()`** ([maws/space.py:270](maws/space.py#L270), [maws/space.py:339](maws/space.py#L339)) — both classes returned coordinates around the origin instead of `self.centre`. Any caller passing a non-zero centre (e.g. the ligand COM in `maws2023.py`) was silently sampling in the wrong region of space. `Box.generator()` already does this correctly; matched the same pattern. - **Biased radial + polar-angle sampling in `SphericalShell.generator()`** — `r ~ Uniform(R_in, R_out)` over-weights the inner shell (no r² Jacobian) and `psi ~ Uniform(0, π)` over-weights the poles. `Sphere.generator()` was already corrected in a prior change; brought `SphericalShell` to the same volume-correct form: - `r = (u·(R_out³ − R_in³) + R_in³)^(1/3)` for u ∈ [0, 1] - direction via `cos(psi)` uniform in [-1, 1] These matter because MAWS energy/entropy averages over the sampled poses (`S(energies, beta=BETA)` in the first MAWS step) — biased sampling biases the result. ## Tests Existing `tests/test_space.py` only used `centre=[0,0,0]`, which masked the centre-offset bug entirely. Added tests that: - Pin samples within the radius/shell measured from a non-zero centre (deterministic, would fail before fix). - Statistically pin the radial mean (E[r] ≈ 8.04 vs. biased 7.5) and direction (E[(z/r)²] ≈ 1/3 vs. biased 0.50) at N=10_000 with a tolerance band that cleanly separates correct from biased. `pytest tests/test_space.py` — 17/17 pass after the fix. ## Test plan - [ ] CI passes on `main` - [ ] `pytest tests/test_space.py` locally — 17 passed - [ ] No callers regressed (Sphere/SphericalShell are not currently invoked in `maws2023.py`; `Box`/`Cube` paths unchanged) 🤖 Generated with [Claude Code](https://claude.com/claude-code)
1 parent 242ceaf commit 8825607

2 files changed

Lines changed: 88 additions & 8 deletions

File tree

maws/space.py

Lines changed: 14 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -269,9 +269,9 @@ def generator(self):
269269
rotation_angle = np.random.uniform(0, 2 * np.pi)
270270
result = np.array(
271271
[
272-
r * np.cos(phi) * np.sin(psi),
273-
r * np.sin(phi) * np.sin(psi),
274-
r * np.cos(psi),
272+
self.centre[0] + r * np.cos(phi) * np.sin(psi),
273+
self.centre[1] + r * np.sin(phi) * np.sin(psi),
274+
self.centre[2] + r * np.cos(psi),
275275
x,
276276
y,
277277
z,
@@ -332,15 +332,21 @@ def generator(self):
332332
"""
333333
axis = np.random.uniform(-1, 1, 3)
334334
x, y, z = axis / np.linalg.norm(axis)
335-
r = np.random.uniform(self.innerRadius, self.outerRadius)
335+
# Uniform-in-volume radial sampling: r ~ (u·(R_out³ − R_in³) + R_in³)^(1/3)
336+
u = np.random.uniform(0, 1)
337+
r = (u * (self.outerRadius**3 - self.innerRadius**3) + self.innerRadius**3) ** (
338+
1 / 3
339+
)
336340
phi = np.random.uniform(0, 2 * np.pi)
337-
psi = np.random.uniform(0, np.pi)
341+
# Uniform direction on the sphere: cos(psi) uniform in [-1, 1]
342+
cos_psi = np.random.uniform(-1, 1)
343+
psi = np.arccos(cos_psi)
338344
rotation_angle = np.random.uniform(0, 2 * np.pi)
339345
result = np.array(
340346
[
341-
r * np.cos(phi) * np.sin(psi),
342-
r * np.sin(phi) * np.sin(psi),
343-
r * np.cos(psi),
347+
self.centre[0] + r * np.cos(phi) * np.sin(psi),
348+
self.centre[1] + r * np.sin(phi) * np.sin(psi),
349+
self.centre[2] + r * np.cos(psi),
344350
x,
345351
y,
346352
z,

tests/test_space.py

Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -93,6 +93,26 @@ def test_sphere_sample_within_radius(self):
9393
distance = np.linalg.norm(pos)
9494
assert distance <= 10.0, f"Sample {pos} outside sphere radius"
9595

96+
def test_sphere_sample_offset_by_centre(self):
97+
"""Sphere samples are within radius of the (non-zero) centre, not the origin."""
98+
centre = np.array([50.0, -30.0, 12.0])
99+
sphere = Sphere(radius=10.0, centre=centre)
100+
for _ in range(50):
101+
sample = sphere.generator()
102+
pos = np.array(sample[:3])
103+
distance = np.linalg.norm(pos - centre)
104+
assert distance <= 10.0 + 1e-9, (
105+
f"Sample {pos} not within radius of centre {centre}"
106+
)
107+
108+
def test_sphere_is_in_membership(self):
109+
"""Sphere.is_in() correctly accepts inside points and rejects outside ones."""
110+
centre = np.array([10.0, 10.0, 10.0])
111+
sphere = Sphere(radius=5.0, centre=centre)
112+
assert bool(sphere.is_in(np.array([10.0, 10.0, 10.0])))
113+
assert bool(sphere.is_in(np.array([13.0, 10.0, 10.0])))
114+
assert not bool(sphere.is_in(np.array([100.0, 0.0, 0.0])))
115+
96116

97117
class TestSphericalShell:
98118
"""Tests for SphericalShell sampling space."""
@@ -112,6 +132,60 @@ def test_shell_sample_within_bounds(self):
112132
distance = np.linalg.norm(pos)
113133
assert 5.0 <= distance <= 10.0, f"Sample {pos} outside shell bounds"
114134

135+
def test_shell_sample_offset_by_centre(self):
136+
"""SphericalShell samples sit in the shell measured from centre, not origin."""
137+
centre = np.array([50.0, -30.0, 12.0])
138+
shell = SphericalShell(innerRadius=5.0, outerRadius=10.0, centre=centre)
139+
for _ in range(50):
140+
sample = shell.generator()
141+
pos = np.array(sample[:3])
142+
distance = np.linalg.norm(pos - centre)
143+
assert 5.0 - 1e-9 <= distance <= 10.0 + 1e-9, (
144+
f"Sample {pos} outside shell of centre {centre}"
145+
)
146+
147+
def test_shell_radial_distribution_is_uniform_in_volume(self):
148+
"""Shell sampling must be uniform per unit volume (∝ r²), not uniform in r.
149+
150+
For R_in=5, R_out=10 with the volume-correct CDF
151+
r = (u·(R_out³ − R_in³) + R_in³)^(1/3),
152+
E[r] = (3/4)·(R_out⁴ − R_in⁴)/(R_out³ − R_in³) ≈ 8.036.
153+
A naive `r ~ Uniform(R_in, R_out)` gives E[r] = 7.5 — well outside
154+
the statistical band of the correct sampler at N=10_000.
155+
"""
156+
rng = np.random.default_rng(0)
157+
np.random.seed(0) # generator() uses np.random
158+
R_in, R_out = 5.0, 10.0
159+
shell = SphericalShell(innerRadius=R_in, outerRadius=R_out, centre=[0, 0, 0])
160+
radii = np.array([np.linalg.norm(shell.generator()[:3]) for _ in range(10_000)])
161+
expected = (3 / 4) * (R_out**4 - R_in**4) / (R_out**3 - R_in**3)
162+
# Tight band: correct sampler lands within ~0.05 of 8.036; biased
163+
# sampler is at 7.5, well outside this band.
164+
assert abs(radii.mean() - expected) < 0.1, (
165+
f"E[r] = {radii.mean():.3f}, expected ~{expected:.3f} (uniform-in-r bias?)"
166+
)
167+
_ = rng # silence unused
168+
169+
def test_shell_direction_is_uniform_on_sphere(self):
170+
"""Shell directions must be uniform on the sphere (E[z²/r²] = 1/3).
171+
172+
With `psi ~ Uniform(0, π)` the polar angle is biased toward the
173+
poles, giving E[cos²(psi)] = 1/2 instead of 1/3.
174+
"""
175+
np.random.seed(1)
176+
R_in, R_out = 5.0, 10.0
177+
shell = SphericalShell(innerRadius=R_in, outerRadius=R_out, centre=[0, 0, 0])
178+
samples = np.array([shell.generator()[:3] for _ in range(10_000)])
179+
# Direction vectors (unit), then mean of z²
180+
radii = np.linalg.norm(samples, axis=1, keepdims=True)
181+
unit_z = (samples / radii)[:, 2]
182+
mean_z2 = float((unit_z**2).mean())
183+
# Uniform-on-sphere: 1/3 ≈ 0.333. Biased: 1/2 = 0.500. Tolerance
184+
# 0.04 cleanly separates the two at N=10_000.
185+
assert abs(mean_z2 - 1 / 3) < 0.04, (
186+
f"E[(z/r)²] = {mean_z2:.3f}, expected ~0.333 (polar-angle bias?)"
187+
)
188+
115189

116190
class TestNAngles:
117191
"""Tests for NAngles (torsion angle) sampling space."""

0 commit comments

Comments
 (0)