Does Shinka Rediscover My SoA Tricks? A Python N-body Warmup

The last post ended with a promise: the circle-packing run proved ShinkaEvolve + Meridian + my Claude Code plan works for zero dollars, but I did not care whether a circle-packing search beats the literature. What I wanted was to point the pipeline at code I wrote and see whether it rediscovers the patterns I hand-tuned in the AoSoA post — for_each_multistream<8>, per-field AVX2 accumulators, the lambda-driven for_each that kills the iterator branch.

This post is not that experiment. This is the cheaper warmup before it.

The sanity check: given a small, naive Python HPC kernel whose optimal form I already know, does Shinka find it? If the answer is no on a kernel I can reason about in my head, there is no point pointing the tool at 500 lines of C++ template metaprogramming. If the answer is yes, the next run has a shot.

So I wrote the smallest meaningful HPC kernel I could think of, seeded it as a Python list of dicts, and ran 30 generations of Haiku 4.5 against it.

The Seed

An N-body simulation frame, stripped to the three operations that matter: integrate (pos += vel * dt, life -= dt), kinetic energy (sum(0.5 * m * v^2)), and cull (drop life <= 0).

The internal storage is a Python list of dicts, one per particle with keys x, y, z, vx, vy, vz, mass, life. The hot loop:

for p in particles:
    p['x'] += p['vx'] * dt
    p['y'] += p['vy'] * dt
    p['z'] += p['vz'] * dt
    p['life'] -= dt

Roughly the worst representation you can pick on a modern CPU. Every field access is a dict lookup, every op goes through the interpreter, nothing vectorizes. Naive seed: 48 ms for 20 frames at N=10000, or 2.4 ms per frame. Score 20.6 (score = 1 / min_seconds).

The benchmark: N=10000, 20 frames per rep, 5 reps, min wall time. Validation compares summary scalars plus one frame of kinetic energy against reference at atol=1e-6. The three ops and the init go inside an EVOLVE-BLOCK-START/END pair; the harness is locked outside.

Config for the run:

Parameter	Value
Generations	30
Islands	1
Archive size	20
Model	Claude Haiku 4.5 (claude-haiku-4-5-20251001)
Patch types	70% diff, 30% full rewrite
Routing	Meridian → Claude Code OAuth (zero API spend)

Why Closed-Form Init Matters

The single most transferable thing I learned on this run.

The first version of init_particles used np.random.RandomState(seed) with per-particle scalar draws in a for-loop. Looked fine, passed my unit tests. I launched Shinka, and the first SoA candidate came back as "wrong physics" — sums off by order 1.

The physics was not wrong. The candidate had replaced the per-particle scalar draws with one vectorized rng.random(N) per field. Mathematically correct — but rng.random(N) consumes the PRNG stream in a different order than 8 × N interleaved scalar draws. Same bits, different particles. The validator only sees sum_x off and reports "wrong physics." A correct candidate got rejected on representation ambiguity.

The fix: delete the RNG and use a closed-form deterministic recipe.

phi = 0.6180339887498949   # golden ratio fraction
sq2 = 0.4142135623730951   # sqrt(2) - 1
particles.append({
    'x':  math.fmod(fi * phi + s * 0.01, 1.0),
    'y':  math.fmod(fi * sq2 + s * 0.02, 1.0),
    'vx': 0.1 * math.sin(fi * 0.7 + s * 0.11),
    # ...
})

Pure algebra and trig. The particles look quasi-random because of the irrational multipliers, but every value is a deterministic function of index i. A candidate can compute them in a Python for-loop, a numpy vectorized expression, or a Numba JIT function and get bit-identical results per particle. I checked math.fmod vs np.fmod(float64, float64) and math.sin vs np.sin are bit-identical on this CPython+numpy install.

Without this fix the experiment would have been impossible. Every SoA candidate would have failed validation as "wrong physics" when all that changed was the arithmetic order of the init.

The transferable lesson: your reference values must be computable identically across representations, OR your validator must be representation-aware. The cleanest fix is the former. RNG-seeded init is a trap. Closed-form init is not.

Score Timeline

View larger ↗

Valid candidates in color, invalid as hollow markers. Gen 1 is the big one-shot jump. Gens 2-21 plateau at ~800 (42x). Gen 25 breaks it, gen 27 sets the final best at 1599 (80x), gen 29 lands on top of it from a different lineage.

View larger ↗

Three plateaus: 1x → 38x → 42x → 80x. The first gap took one generation. The second took 20. The third took one generation again, after a specific environmental change mid-run.

The 38x One-shot

Gen 1, the first diff Haiku proposed, replaced the list-of-dicts with a dict of numpy arrays and rewrote run_frame as four in-place vectorized ops: x += vx * dt, np.sum(0.5 * mass * (vx*vx + vy*vy + vz*vz)), boolean-mask cull. Score went from 20.6 to 763.6. 38.2x over the seed, from a single diff patch, on the very first attempt.

This is the standard SoA-vectorized pattern — the first thing you teach in a numpy course. Haiku produced it correctly the first time. The only "cheat" is that the closed-form init gave bit-identical particles under vectorization, so validation passed.

The 42x Plateau

After gen 1, Haiku spent thirteen generations refining the same SoA-numpy approach. Each step was small:

Gen	Speedup	Change
1	38.2x	dict-of-arrays, vectorized integrate + KE + cull
4	41.1x	np.sum(0.5*m*v2) → 0.5 * np.dot(m, vx*vx+vy*vy+vz*vz)
10	41.8x	hoist field references into Python locals before the op
21	42.3x	hoist everything incl. cull; explicit dict literal instead of comprehension

Gen 4 is the most interesting. Swapping np.sum(0.5 * mass * v_squared) for 0.5 * np.dot(mass, vx*vx + vy*vy + vz*vz) is a ~3% win because np.dot calls into BLAS and fuses the multiply-reduce with no intermediate allocation; np.sum(x*y) materializes x*y as a temporary and then reduces. BLAS beats generic numpy on this shape every time.

Gens 10 and 21 are pure bytecode micro-optim — LOAD_FAST on a local is cheaper than LOAD_ATTR + BINARY_SUBSCR on a dict field. Per-op negligible, over 20 × 3 × 10000 ops it adds up to a fraction of a percent.

This is the SoA-numpy local optimum. Thirteen generations to walk 38x → 42x is the diminishing-returns plateau when the representation is fixed and the hot loop is already vectorized. What broke the plateau was a different representation.

The Numba Subplot

Haiku tried to import numba six times in a row — gens 6, 8, 9, 11, 15, 19 — and crashed every time with No module named 'numba'. I had not installed Numba.

Six attempts is a strong signal. Once is a coincidence; twice is a preference; six times is the LLM telling me "the next step is JIT compilation and your environment is broken." The meta-summarizer even started naming candidates numba_jit_compilation_attempt by gen 15.

I installed Numba mid-run. Between gen 23 and gen 24, pip install numba into the venv (numba 0.65.0 + llvmlite 0.47.0).

And the first Numba candidate was slower than the plateau. Gen 24 scored 690 (34.5x) vs the 42.3x best. It used @numba.njit(parallel=True) with a prange loop over N particles. For N=10000, the cost of spawning OpenMP threads and coordinating the reduction exceeds the work. Classic "parallelized too small" mistake.

Gen 25 fixed it. It was a full rewrite with parent gen 21 — the last pre-Numba SoA-numpy best — and had not seen gen 24's failed attempt. It independently landed on @numba.njit(fastmath=True) with a serial scalar loop. Score: 1045, 52x. First real win.

Gen 27 (parent gen 25) fused the cull step into the JIT function: count alive in a first pass inside @njit, allocate output arrays at exactly that size, single copy pass. Score: 1599, 80x. New best, headline result.

The winning function looks roughly like this:

@numba.njit(fastmath=True)
def _frame(x, y, z, vx, vy, vz, mass, life, dt):
    n = x.shape[0]
    ke = 0.0
    alive = 0
    for i in range(n):
        x[i] += vx[i] * dt
        y[i] += vy[i] * dt
        z[i] += vz[i] * dt
        life[i] -= dt
        ke += 0.5 * mass[i] * (vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i])
        if life[i] > 0.0:
            alive += 1
    nx = np.empty(alive, dtype=x.dtype)  # ny, nz, nvx, ... same
    j = 0
    for i in range(n):
        if life[i] > 0.0:
            nx[j] = x[i]  # ... same for all fields
            j += 1
    return nx, ny, nz, nvx, nvy, nvz, nmass, nlife, ke

Three things do the work. Fusion: integrate and KE share the same index and reads; combining them halves cache misses. Serial, not parallel: no prange. Fastmath: lets LLVM reorder the accumulator and vectorize the inner FMA. At N=10000 × 8 × 8 bytes, the working set (~640 KB) is L2-resident.

To within a constant factor this is what I would have written in C with #pragma omp simd reduction(+:ke). Haiku got there from a Python list of dicts in 27 generations.

Two Independent Branches Converged on 80x

View larger ↗

Gen 27 (best, 80x) has parent gen 25 → gen 21. One clean chain: SoA-numpy → Numba serial + fastmath → + fused cull.

Gen 29 (79.8x) was a full rewrite with parent gen 18 — back on the SoA-numpy side of the archive, and the chain from gen 18 upward had never seen a Numba candidate. Gen 29 arrived at the same Numba structure from a completely different lineage.

The two candidates are within 4 in score and use the same primitives. Under this benchmark, "reach for @njit(fastmath=True) and fuse the cull" is a fixed point — if you arrive at the Numba family at all, you land on the same structure regardless of path. Strong signal that the search actually found the local optimum rather than hitting it by luck.

This is the only thing in the run that looks like independent evolutionary convergence rather than retrieval. The 38x one-shot is retrieval from training data. The 80x convergence from two separate branches is not.

Two Structured-array Regressions

Gens 12 and 20 both regressed to ~16x. Both replaced the dict-of-arrays with a single np.dtype structured array of 8 f8 fields.

This looks like SoA. It reads like SoA. It is not SoA. A structured numpy array with 8 f8 fields is internally a single contiguous buffer of 64-byte records — an AoS. particles['x'] returns a non-contiguous strided view with stride 64 bytes; element-wise ops degrade to gather/scatter and cache lines are 1/8 useful. ~3x slower than dict-of-arrays.

Haiku walked into this trap twice — gen 12 as a diff, gen 20 as an independent full rewrite. The archive rejected both, but the LLM did not learn. Retrieval-flavored mutations include patterns that look like the right move and are not, and the only way to reject them is to actually benchmark them.

What Went Wrong in the Invalid Half

View larger ↗

Out of 30 candidates, 22 were valid and 8 were not:

Failure mode	Count	Cause
No module named 'numba'	6	numba not in the venv (gens 6, 8, 9, 11, 15, 19)
Python scoping bug	2	UnboundLocalError / NameError in the evolve block
Wrong physics	0	none — closed-form init paid off

Six of the eight failures are my fault. Notable is the zero count on "wrong physics." On the circle-packing run the dominant failure was "right answer, wrong epsilon." Here, with closed-form init, it is "your environment is broken" and "the LLM wrote a typo" — both framework-level, not modeling. Closed-form init collapsed an entire category of spurious rejections.

Precision

View larger ↗

Every valid candidate matches the reference kinetic energy to at most 1 ULP — 1.11e-16, well inside atol=1e-6. Across 22 valid candidates spanning 80x in wall time, the Pareto frontier collapses to a horizontal line at "essentially exact." No speed-for-precision tradeoff.

The bit-identical candidates all use a reduction that walks indices 0..N-1 once, either as 0.5 * np.dot(mass, v_squared) or as a Numba scalar accumulator. The 1-ULP candidates all use np.sum(0.5 * mass * v_squared), which internally uses pairwise summation and lands 1 ULP from the seed's strict left-to-right accumulation.

I am calling this out because fastmath=True is the usual "does this break precision?" boogeyman, and on this kernel the answer is no. Numba's fastmath flag on the fused integrate+KE loop does not fire any reorder that visibly changes the 53-bit mantissa. 80x faster is free in precision terms. On a kernel with catastrophic cancellation the answer would be different; on a textbook N-body KE at N=10000 it is not.

Cost

Same as the previous post. Zero API spend.

Resource	Consumed
Wall-clock	~40 minutes
Per generation	~80 seconds
Anthropic API spend	$0.00
Claude Max messages	~30-50 (inside the 5h rolling quota)

Meridian routed every Haiku 4.5 call through my Claude Code OAuth token. The only change from the circle-packing run was the seed file.

What This Proves and What It Does Not

Proves:

1. On a small Python HPC kernel whose optimal form is a known pattern (SoA + vectorize + JIT), Haiku 4.5 + ShinkaEvolve finds that form in one shot for 38x and grinds out the rest of the SoA-numpy plateau in 13 more generations.
2. Given the right environment, Haiku finds Numba on its own (after six import failures), correctly identifies fastmath=True and single-threaded as the right flags, and discovers the fuse-cull-into-JIT trick that gets from 52x to 80x.
3. Two independent lineages converged on the same 80x answer. The search actually found the local optimum rather than hitting it by accident.
4. Closed-form deterministic init is the cleanest way to make a numerical-kernel validator representation-agnostic.

Does not prove:

1. ShinkaEvolve is going to rediscover my hand-tuned AoSoA tricks in C++. Different problem: the seed is larger, the optimal form depends on hardware and workload class, and patterns like for_each_multistream<N> and per-field accumulators are not in undergraduate numpy tutorials. Here, Haiku never had to invent anything — it retrieved the right pattern and applied it cleanly. On C++, retrieval might not be enough.
2. 80x on a Python N-body frame is impressive. It is not. Any competent HPC engineer gets 80x on this kernel in 30 minutes by hand. The point is not the 80x; the point is that the loop got there from a list of dicts with zero human intervention, correctly, for zero dollars, with a convergent lineage and in-spec precision.
3. Haiku is the right model for the C++ run. Here, retrieval-class mutations were enough. On C++, the harder step is going to be invention, and I may need Sonnet. The circle-packing near-miss at 99.36% of SOTA was also a Haiku+precision problem. Consistent signal about where Haiku's ceiling sits.

The most embarrassing learning: six generations wasted on missing Numba. Pre-install every JIT and optimizer the LLM might plausibly reach for — Numba, Cython, cffi, maybe Taichi. The cost of an unused package is zero. The cost of six generations crashing on import is six generations.

What's Next

This was the warmup. The actual experiment is still the one I promised: point ShinkaEvolve at the AoSoA container from the SoA-vs-AoS post, give it the multi-stream and per-field-accumulator benchmarks as the scoring function, and see whether it rediscovers for_each_multistream<8>. C++, seed an order of magnitude larger.

The questions I now know how to ask:

• How much is retrieval vs invention? Here it was mostly retrieval — SoA, np.dot, @njit(fastmath) are in every HPC Python tutorial. On C++, multi-stream unrolling for L3-bound reductions is a trick I picked up by reading perf counters, not from a class. If Shinka rediscovers it, that is invention.
• Does validator strictness matter? Here, closed-form init made the validator transparent. On C++ the "correct physics" check needs to be representation-agnostic across AoSoA stream counts. I know how to write that now.
• Does the lineage converge? Two branches converging on 80x was the most confidence-inspiring observation. If multiple C++ lineages independently discover multi-stream unrolling, I will believe the result.

I will report back.