From 1a2344f0292fec579a00fc950a20e783711bea08 Mon Sep 17 00:00:00 2001 From: Lehman Garrison Date: Wed, 21 Oct 2026 15:28:39 -0400 Subject: [PATCH 1/5] benchmarks: time potential evaluation up to 10,000 points At n = 1 and 10, each potential benchmark measures only the fixed cost of a call (unit handling and validation), and n = 1,000 is still mostly fixed cost for cheap potentials. Replace n = 10 with n = 10,000 so the per-point cost of each evaluation path shows. The classes with expensive potentials keep their explicit size caps. Generate the positions in a fixture, so that the timed region is the evaluation alone: at n = 10,000, generating them costs about as much as evaluating a cheap potential. With a fixed seed per test, every test also gets the same positions, regardless of which tests ran before it. Co-Authored-By: Claude Opus 5.5 --- tests/benchmarks/test_potentials_benchmark.py | 33 ++++++++++--------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/tests/benchmarks/test_potentials_benchmark.py b/tests/benchmarks/test_potentials_benchmark.py index c3181a156..6b4072417 100644 --- a/tests/benchmarks/test_potentials_benchmark.py +++ b/tests/benchmarks/test_potentials_benchmark.py @@ -28,27 +28,28 @@ class BenchmarkPotentialBase: tested (e.g., for expensive potentials that are too slow at large n_points). """ - @pytest.fixture(params=[1, 10, 1_000]) + @pytest.fixture(params=[1, 1_000, 10_000]) def n_points(self, request): return request.param - @pytest.fixture(scope="class") - def rng(self): - return np.random.default_rng(42) + @pytest.fixture + def q(self, potential, n_points): + # A fixture, so that generating the positions is not timed + return np.random.default_rng(42).normal(0, 10, size=(potential.ndim, n_points)) @pytest.mark.benchmark - def test_energy(self, potential, n_points, rng): - potential.energy(rng.normal(0, 10, size=(potential.ndim, n_points))) + def test_energy(self, potential, q): + potential.energy(q) @pytest.mark.benchmark - def test_gradient(self, potential, n_points, rng): - potential.gradient(rng.normal(0, 10, size=(potential.ndim, n_points))) + def test_gradient(self, potential, q): + potential.gradient(q) @pytest.mark.benchmark - def test_density(self, potential, n_points, rng): + def test_density(self, potential, q): if type(potential).__name__ in NO_DENSITY: pytest.skip("density not implemented for this potential") # type: ignore[call-arg] - potential.density(rng.normal(0, 10, size=(potential.ndim, n_points))) + potential.density(q) # ============================================================================ @@ -139,16 +140,16 @@ def potential(self): ) @pytest.mark.benchmark - def test_energy(self, potential, n_points, rng): - potential.energy(rng.normal(0, 10, size=(potential.ndim, n_points)), t=2.5) + def test_energy(self, potential, q): + potential.energy(q, t=2.5) @pytest.mark.benchmark - def test_gradient(self, potential, n_points, rng): - potential.gradient(rng.normal(0, 10, size=(potential.ndim, n_points)), t=2.5) + def test_gradient(self, potential, q): + potential.gradient(q, t=2.5) @pytest.mark.benchmark - def test_density(self, potential, n_points, rng): - potential.density(rng.normal(0, 10, size=(potential.ndim, n_points)), t=2.5) + def test_density(self, potential, q): + potential.density(q, t=2.5) @pytest.mark.skipif(not EXP_ENABLED, reason="requires EXP") From 0a350afe64531d5128d9fef2c39936dd3cbc8982 Mon Sep 17 00:00:00 2001 From: Lehman Garrison Date: Wed, 21 Oct 2026 15:28:49 -0400 Subject: [PATCH 2/5] benchmarks: add a hessian benchmark Hessians had no benchmark. One cheap potential is enough to track the cost of the evaluation machinery around the kernels. Co-Authored-By: Claude Opus 5.5 --- tests/benchmarks/test_potentials_benchmark.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/tests/benchmarks/test_potentials_benchmark.py b/tests/benchmarks/test_potentials_benchmark.py index 6b4072417..7dd1a2d53 100644 --- a/tests/benchmarks/test_potentials_benchmark.py +++ b/tests/benchmarks/test_potentials_benchmark.py @@ -69,6 +69,18 @@ class TestCanonicalPotentialsBenchmark(BenchmarkPotentialBase): pass +@pytest.fixture +def positions(n_points): + return np.random.default_rng(42).normal(0, 10, size=(3, n_points)) + + +@pytest.mark.benchmark +@pytest.mark.parametrize("n_points", [1, 1_000, 10_000]) +def test_hessian(positions): + # Just one representative potential + CANONICAL["HernquistPotential"].hessian(positions) + + # ============================================================================ # Special (file-based potentials with expensive setup) From 60852d26ff03e86ba2e8d93122646c443ef5571c Mon Sep 17 00:00:00 2001 From: Lehman Garrison Date: Wed, 21 Oct 2026 15:29:23 -0400 Subject: [PATCH 3/5] benchmarks: integrate n_orbits orbits in the integration benchmarks Both integration benchmarks took an n_orbits parameter but always integrated a single orbit, so each case ran four times under different names. Integrate n_orbits copies of the initial condition, and skip the cases with more than 10^6 orbit-steps, which are too slow under CodSpeed's simulation. Build the initial conditions in fixtures, so that only the integration is timed. Co-Authored-By: Claude Opus 5.5 --- tests/benchmarks/test_integrate_benchmark.py | 39 ++++++++++++++++---- 1 file changed, 32 insertions(+), 7 deletions(-) diff --git a/tests/benchmarks/test_integrate_benchmark.py b/tests/benchmarks/test_integrate_benchmark.py index a8aaed410..d0ff18d8f 100644 --- a/tests/benchmarks/test_integrate_benchmark.py +++ b/tests/benchmarks/test_integrate_benchmark.py @@ -1,4 +1,5 @@ import astropy.units as u +import numpy as np import pytest import gala.dynamics as gd @@ -10,21 +11,45 @@ def potential(): return gp.HernquistPotential(m=1e11, c=5, units="galactic") +def make_w0(v_y, n_orbits, n_steps): + if n_orbits * n_steps > 1_000_000: + pytest.skip("too slow under CodSpeed's simulation") + xyz = np.zeros((3, n_orbits)) + xyz[0] = 10.0 + vxyz = np.zeros((3, n_orbits)) + vxyz[1] = v_y + return gd.PhaseSpacePosition(xyz * u.kpc, vxyz * u.km / u.s) + + +@pytest.fixture +def circular_w0(n_orbits, n_steps): + return make_w0(138.25768647, n_orbits, n_steps) + + +@pytest.fixture +def eccentric_w0(n_orbits, n_steps): + return make_w0(0.25 * 138.25768647, n_orbits, n_steps) + + @pytest.mark.benchmark @pytest.mark.parametrize("Integrator", ["dop853", "leapfrog", "ruth4"]) @pytest.mark.parametrize("n_steps", [1, 10, 1_000, 100_000]) @pytest.mark.parametrize("n_orbits", [1, 10, 1_000, 100_000]) -def test_integrate_circular_orbit(potential, Integrator, n_steps, n_orbits): - w0 = gd.PhaseSpacePosition([10.0, 0, 0] * u.kpc, [0, 138.25768647, 0] * u.km / u.s) - potential.integrate_orbit(w0, dt=1 * u.Myr, n_steps=n_steps, Integrator=Integrator) +def test_integrate_circular_orbit( + potential, circular_w0, Integrator, n_steps, n_orbits +): + potential.integrate_orbit( + circular_w0, dt=1 * u.Myr, n_steps=n_steps, Integrator=Integrator + ) @pytest.mark.benchmark @pytest.mark.parametrize("Integrator", ["dop853", "leapfrog", "ruth4"]) @pytest.mark.parametrize("n_steps", [1, 10, 1_000, 100_000]) @pytest.mark.parametrize("n_orbits", [1, 10, 1_000, 100_000]) -def test_integrate_eccentric_orbit(potential, Integrator, n_steps, n_orbits): - w0 = gd.PhaseSpacePosition( - [10.0, 0, 0] * u.kpc, [0, 0.25 * 138.25768647, 0] * u.km / u.s +def test_integrate_eccentric_orbit( + potential, eccentric_w0, Integrator, n_steps, n_orbits +): + potential.integrate_orbit( + eccentric_w0, dt=1 * u.Myr, n_steps=n_steps, Integrator=Integrator ) - potential.integrate_orbit(w0, dt=1 * u.Myr, n_steps=n_steps, Integrator=Integrator) From d14a50dc5961e8200d8f46e9cf7de51505d6be28 Mon Sep 17 00:00:00 2001 From: Lehman Garrison Date: Wed, 21 Oct 2026 15:29:35 -0400 Subject: [PATCH 4/5] benchmarks: evaluate energy from Quantity and Orbit inputs The potential benchmarks pass plain arrays already in the potential's units, so they skip the conversions applied to other inputs. Benchmark energy with a Quantity, with an integrated Orbit, and through Orbit.energy(). The orbit has 10^6 points: at 10^4, the copies made converting it cost a few microseconds, hidden by the ~0.3 ms of fixed overhead per call, while at 10^6 they take most of the time: Orbit.energy() takes 69 ms, against 15 ms for the same positions as an array. Co-Authored-By: Claude Opus 5.5 --- tests/benchmarks/test_potentials_benchmark.py | 38 +++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/tests/benchmarks/test_potentials_benchmark.py b/tests/benchmarks/test_potentials_benchmark.py index 7dd1a2d53..baa0635e0 100644 --- a/tests/benchmarks/test_potentials_benchmark.py +++ b/tests/benchmarks/test_potentials_benchmark.py @@ -6,6 +6,7 @@ import pytest from gala._cconfig import EXP_ENABLED, GSL_ENABLED +import gala.dynamics as gd import gala.potential as gp from gala.units import SimulationUnitSystem, galactic @@ -200,3 +201,40 @@ def potential(self): coef_file=EXP_SINGLE_COEF_FILE, units=exp_units, ) + + +# ============================================================================ +# Input types other than plain arrays, which are converted before evaluation + + +@pytest.fixture +def positions_kpc(positions): + return positions * u.kpc + + +@pytest.mark.benchmark +@pytest.mark.parametrize("n_points", [1, 10_000]) +def test_energy_quantity(positions_kpc): + CANONICAL["HernquistPotential"].energy(positions_kpc) + + +@pytest.fixture(scope="module") +def orbit(): + # 1,000 orbits x 1,000 times: large enough that copies of the orbit cost more + # than the fixed overhead of a call + rng = np.random.default_rng(42) + w0 = gd.PhaseSpacePosition( + pos=rng.normal(0, 10, size=(3, 1000)) * u.kpc, + vel=rng.normal(0, 100, size=(3, 1000)) * u.km / u.s, + ) + return CANONICAL["HernquistPotential"].integrate_orbit(w0, dt=1.0, n_steps=999) + + +@pytest.mark.benchmark +def test_energy_orbit(orbit): + CANONICAL["HernquistPotential"].energy(orbit) + + +@pytest.mark.benchmark +def test_orbit_energy_method(orbit): + orbit.energy() From 365c12f10f973cff0adcb944a62cc7c48a9a90ca Mon Sep 17 00:00:00 2001 From: Lehman Garrison Date: Wed, 21 Oct 2026 14:49:28 -0400 Subject: [PATCH 5/5] benchmarks: time O2GF on one orbit under CodSpeed The O2GF benchmarks with 10 orbits took 49 and 33 minutes of the CodSpeed run, about 300 times their real time (11 s and 5.5 s), because most of O2GF's time is the Python-level fit of the toy potential. Each orbit is fit separately, so the cost is linear in the number of orbits, and timing 10 orbits measures the same code as timing one. Skip the multi-orbit cases under CodSpeed. They still run in the walltime scaling benchmarks, which plot the scaling with n_orbits in the docs and take about a minute. Co-Authored-By: Claude Opus 5.5 --- tests/benchmarks/test_actionangle_benchmark.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/tests/benchmarks/test_actionangle_benchmark.py b/tests/benchmarks/test_actionangle_benchmark.py index b7e6b8eab..342b1c4f3 100644 --- a/tests/benchmarks/test_actionangle_benchmark.py +++ b/tests/benchmarks/test_actionangle_benchmark.py @@ -1,3 +1,5 @@ +import os + import astropy.units as u import numpy as np import pytest @@ -35,7 +37,14 @@ def test_staeckel_fudge(potential, benchmark, n_points): @pytest.mark.parametrize("n_steps", [1_000, 10_000]) @pytest.mark.parametrize("n_orbits", [1, 10]) -def test_o2gf(potential, benchmark, n_orbits, n_steps): +def test_o2gf(request, potential, benchmark, n_orbits, n_steps): + # Codspeed takes an hour for 10 orbits; stick to 1. + codspeed = ( + request.config.getoption("--codspeed", False) or "CODSPEED_ENV" in os.environ + ) + if codspeed and n_orbits > 1: + pytest.skip("multiple orbits are only timed for the scaling benchmarks") + # integrate outside the timed region: integration has its own benchmarks orbits = potential.integrate_orbit( make_w0(potential, n_orbits), dt=1 * u.Myr, n_steps=n_steps