tboudreaux/MeanField
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Preserve completed weak-response study while retaining Q3 baseline

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-18 14:27:15 -04:00
parent 211be17190
commit 6ceda61bdd
3 changed files with 525 additions and 18 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
experiments/plot_shape_surface_comparison.py
View File

@@ -28,6 +28,8 @@ def main():
raise ValueError("This comparison requires n=1.5, J=0, and Q3 followed by Q4 displacement.")
with (directory / 'weak_shape_physical_surface_nodes.csv').open() as stream:
rows = list(csv.DictReader(stream))
with (directory / 'newton.csv').open() as stream:
accepted = sum(int(row['accepted']) for row in csv.DictReader(stream))
points = np.array([[float(row['unit_radial_'+a]) for a in 'xyz'] for row in rows])
radii = np.array([float(row['physical_radius']) for row in rows])
if not np.isfinite(points).all() or not np.isfinite(radii).all() or np.any(radii <= 0):
@@ -35,20 +37,21 @@ def main():
latitude = np.arcsin(np.clip(points[:, 2], -1, 1))
longitude = np.arctan2(points[:, 1], points[:, 0])
deviation = 100 * (radii / np.mean(radii) - 1)
data.append((longitude, latitude, deviation, len(rows)))
data.append((longitude, latitude, deviation, len(rows), accepted))
limit = max(float(np.max(np.abs(row[2]))) for row in data)
if limit == 0:
limit = 1e-12
fig, axes = plt.subplots(1, 2, figsize=(12, 5), subplot_kw={'projection': 'mollweide'}, layout='constrained')
for axis, (longitude, latitude, deviation, count), title in zip(axes, data, ('Q3 displacement / surface', 'Q4 displacement / surface')):
fig.get_layout_engine().set(rect=(0, .07, 1, .93))
for axis, (longitude, latitude, deviation, count, accepted), title in zip(axes, data, ('Q3 displacement / surface', 'Q4 displacement / surface')):
picture = axis.scatter(longitude, latitude, c=deviation, s=14, cmap='RdBu_r', vmin=-limit, vmax=limit, linewidths=0)
axis.set_title(f'{title}\n{count} radial nodes')
axis.set_title(f'{title}\n{accepted} accepted updates; {count} radial nodes')
axis.grid(True, alpha=.25)
axis.tick_params(labelsize=8)
fig.colorbar(picture, ax=list(axes), orientation='horizontal', shrink=.65, pad=.08,
label='Radius relative to its unweighted nodal mean (%)')
fig.suptitle('Nonrotating n = 1.5: surface samples after the bounded order trial', fontsize=14)
fig.text(.5, .01, 'Common color scale. Nodal observations only; neither endpoint is certified by this plot.', ha='center', fontsize=9)
fig.text(.5, .01, 'Common color scale; different nodal sampling. Neither endpoint is certified by this plot.', ha='center', fontsize=9)
args.output.parent.mkdir(parents=True, exist_ok=True)
fig.savefig(args.output, dpi=200)
plt.close(fig)