Learning to Speak C & Cython: My GSoC Summer with Astropy

The summer is officially over. I am staring at a remarkably clean Git branch, my laptop didn't literally take off into orbit (though the CPU fans certainly tried a few times during local CI builds), and I somehow know what git rebase -i does without having to Google it in a cold sweat.

If you'd asked me back in May what I was going to be doing, I would have confidently told you I was going to "write tests for Astropy's C extensions." It sounded so neat. So contained. But open source doesn't really work like that. I came in thinking I was just going to write tests, and somewhere along the way, I ended up learning how the actual machinery underneath the Python abstraction works, how maintainers think about architecture, and how to safely catch C-level memory panics without taking down the entire interpreter.

So, here is the real story of what I did for the last few months, what broke, how we fixed it, and where the project stands now.

So, what was I actually supposed to do?

Astropy is a beast of a library. The Python-facing API is incredibly robust and beautifully documented. But underneath all those pretty Python classes is a complex, mixed-language architecture. The library relies heavily on compiled C, and Cython extensions to handle the performance-critical hot-paths.

The problem? That compiled layer was a massive testing blind spot. Before this summer, these performance-critical extensions were almost entirely tested indirectly, meaning they were only validated by calling the high-level Python wrappers. That is a risky abstraction. If a regression happens deep inside the C code, the Python layer sitting above it can accidentally mask it. You wouldn't know something was fundamentally broken until a downstream package started acting weird.

My project goal was to build a dedicated, de novo test suite that bypassed the public API completely and exercised each compiled extension module directly. This wasn't just for code coverage. It was an absolute prerequisite for three massive Astropy milestones:

• The APE Split: Proving the extensions are stable enough to live in a completely separate package.
• Meson Migration: Having a safety net to catch silent regressions when the build system transitions.
• Python 3.15 Limited API: Safely testing internal C-refactors for the upcoming free-threaded builds.

The mission was to harden Astropy's core stability. The reality was peeling away layers of Python abstractions until I was staring at raw memory buffers.

Turns out, "testing C extensions" is not quite as simple as it sounds

I started my summer in astropy/table, which my mentors flagged as a solid starting point. My first target was test_np_utils.py, writing a low-level suite for the join_inner Cython extension.

The thing is, Cython engines don't want your high-level Table objects. They want raw data. I had to build a helper function, _make_join_inputs, which acted as a pre-processor. It took two arrays, concatenated them, sorted them with argsort(kind="stable"), and used a boolean diffs array to find the exact boundaries of unique keys. It output the exact np.intp bindings that Cython demanded.

Once I could feed the raw engine, I built strictly typed dataclasses (ArrayMaskPair and ExpectedResults) to hold the outputs without making the pytest matrices unreadable with massive tuples. I tested standard overlaps, but the real fun was the Cartesian edge cases. I threw O(N²) expansions at the engine, where duplicate keys joined with duplicate keys, just to prove the C-engine could handle the explosive memory allocation without crashing. I also verified that it properly handles np.nan as a unique entity (since nan != nan in Python).

And all this stands as proof that I can actually write code that works, phew!

The C-Slots

Things got significantly more complicated when I moved to test_column_mixins.py. The objective here was to test the Cython __getitem__ routing in pure isolation.

I couldn't just pass a standard array. I had to create a "Shim Strategy" (so fancy), so I wrote MinimalColumn and MinimalMaskedColumn classes that inherited directly from the Cython mixins (_ColumnGetitemShim, _MaskedColumnGetitemShim) and mapped them directly onto raw NumPy arrays using .view() casting. By overriding the .data property to return a pure ndarray view, I stopped the Cython engine from returning raw memoryview buffers, which would have hard-crashed the tests. We were hitting the tp_as_mapping->mp_subscript C-slot directly.

But this is where I hit my first major hurdle.

While I was writing a test for structured arrays, indexing it with a single string, the structured dtype suddenly just dropped, and it returned a 1D array of the underlying field's type. I spent hours going down this rabbit hole, utterly convinced my test was wrong. Eventually, I found the machinery underneath: a literal trapdoor in base_getitem.

if dtype_kind == 'V' and isinstance(item, STRING_TYPES):
    return self.data[item]

Because it accessed the raw .data array, it inherited NumPy's default behavior and lost the structure. Rather than cementing this unintended behavior by writing a passing test for it, I dropped the test case and opened Issue #19827.

This was my first realization that tests aren't just a safety net; when written at this level, they act like a microscope.

And then things started breaking in interesting ways

Open source isn't just writing code, getting a green checkmark, and merging. The messy middle is where the actual work happens. I mean everyone knows that right? Yeah no sometimes people get a reality check; people is me, I am people.

• The Masked Array isinstance Trap: I proudly swapped isinstance() checks for strict type(result) is ... assertions, which instantly caused my masked array tests to fail. It turns out standard ndarray slicing strips subclasses, but np.ma.MaskedArray.__getitem__ inherently preserves them. The C-slot was actively returning the MinimalMaskedColumn shim.
• The CI Gatekeeper: While working on the mixins, I used an @override decorator. Python 3.12 CI matrices immediately failed because typing_extensions wasn't in the base test environment. Instead of bloating pyproject.toml with a new dependency just for one test file, I made the architectural choice to drop the decorator completely. CI went green.
• Git History looking like a crime scene: At one point, a .uv.lock file triggered a check-added-large-files pre-commit hook failure. My branch had some incredibly messy commits. So I took a couple a deep breaths, paced around, and after a midly dramatic breakdown, I ran git restore --staged, and fired up a local interactive rebase (git rebase -i HEAD~6) via Nano to squash the chaos into a single, clean production commit (see told ya it was simple :) )

Somewhere along the way, I stopped just writing tests

As the summer progressed, the project evolved. The APE split required not just test coverage, but rigorous structural boundaries. I started writing .pyi type stubs to define the Python-facing boundaries of compiled extensions so static type checkers could understand them without executing the C-engine, not to forget the Proof of Concept, all in a days work (I'm lying, I almost cried trying the machete strategy).

The XML C-Extension Stubs (iterparser_iterparser.pyi)

When I looked at the C struct for the _iterparser extension, it looked like the parser accepted fd, buffersize, file, and buffer. But looking at the C struct is like looking at a map and realizing the map is wrong.

By analyzing the C-API execution (PyArg_ParseTupleAndKeywords), I discovered that file and buffer were strictly internal C-state variables generated after crossing the Python boundary. The Python __init__ only actually takes fd and buffersize. Because the C-engine didn't natively support keyword argument routing here, I had to enforce positional-only arguments using the PEP 570 / marker in Python.

Even better, we had a friendly discussion about the fd.read callable. My mentor pointed out read(self->file, self->buffer...) in the C code. But tracing the Python C-API showed a Py_BuildValue("(n)", buffersize) call directly before PyObject_CallObject. The C-engine was packing a C Py_ssize_t into a 1-element Python tuple. We proved the Python callable only ever receives one integer, finalizing the strict protocol signature as Callable[[int], bytes].

Dimensionality Locking in _convolve

For _convolve.pyx, we needed to ensure that arrays passed to the C engine matched perfectly before runtime. I used bounded TypeVar aliases (_D1, _D2, _D3) to enforce a strict mathematical rule at the static analysis level. The result, array_to_convolve, and kernel arrays can be 1D, 2D, or 3D, but the linter now mathematically locks them so all three must have the exact same dimensionality.

The unit_list_proxy.c Dependency Injection

This was easily one of the coolest things I worked on. astropy/wcs/src/unit_list_proxy.c had a massive circular dependency problem. It was relying on astropy.units.UnitBase via PyImport_ImportModule at runtime.

We discovered that cunit wasn't a normal attribute; it was a C-level mutable proxy array. Standard Python @property wrappers broke the memory linkage, meaning the C-array never received updates.

Following my mentor Clément's strategy, we scrapped the Python wrappers and built a Dependency Injection architecture. We implemented a _setup_unit_class(PyObject* unit_class) function, dynamically exposed to _wcs. It runs exactly once during Python initialization to cache the Unit class pointer in a static C variable, completely severing the heavy PyImport dependency. We updated the getitem C-slot to natively yield Unit objects (like CompositeUnit) back to Python, and made the setitem slot use duck-typing (parse_strict="warn") so it gracefully emits a UnitsWarning for garbage FITS strings instead of hard-crashing. I was super elated when I finally succeeded and to see the CI checks green? Goodness was that rewarding, so I treated myself to an ice-cream sundae!

Numerical sanity and the "Grid Trap"

I also wrote isolated suites for the astropy.timeseries periodogram C-extensions. Bypassing the high-level BoxLeastSquares API, I fed raw np.float64 memory buffers directly into the Cython boundary (_impl.pyx).

I built a SyntheticTransit fixture to inject a perfect transit dip into a simulated light curve to prove the C-engine's best_objective algorithm peaked accurately.

But math is cruel :(

My tests were failing to recover a mathematically perfect period of 2.0. Why? Because I generated my frequency grid using np.linspace(1.5, 2.5, 100). If you do the floating-point math, the step size is ~0.0101. The number 2.0 literally did not exist in the array (the closest it got was 2.005). I was stuck fighting my own numerical grid. I changed it to 101 steps, and voila it hit the target perfectly!

We also dealt with a fascinating accuracy tradeoff. When using a coarse grid (1,000 points, oversample=10), the discrete phase-binning of the BLS algorithm caused the edges of the transit to smear, changing the recovered depth from a perfect 0.5 to 0.478. My mentor challenged the 5% tolerance. By fixing the linspace bug and adjusting the points, we tightened the mathematical accuracy to 1e-4 while balancing the matrix to keep the isolated test executing in under 300ms so we didn't anger the CI limits.

I ran into a similar trapdoor testing the Lomb-Scargle Cython implementation (cython_impl.pyx). I tried to test the Cython-level t.ndim != 1 error. But when I passed a 2D array for t and 1D arrays for y, np.broadcast_arrays panicked and raised a Python ValueError before the C-engine ever saw it. The fix? I deliberately passed identical 2D arrays for all parameters so it survived the Python broadcasting check and successfully crashed directly on the C-boundary. Hurrayy!

The parts that didn't go according to plan

Of course, not everything matched the original proposal. My initial timeline had me delivering tests for astropy/erfa and astropy.coordinates towards the end of the summer (Deliverables 5 and 6).

The reality of open-source struck: astropy/erfa had already been spun out into the standalone pyerfa package, and astropy.coordinates was delegating its math to that external library. You can't write isolated C-extension tests for an architecture that doesn't live in your repository anymore.

Rather than viewing this as a failure, I contacted my mentors, pivoted, and spent that time significantly deepening the coverage and architectural stability of the .pyi stubs and the XML parsers.

So, where did we end up? (Current State)

The goal was to make the compiled layer testable, and I am leaving behind an architecture to ensure it stays that way.

Merged work & contributions

• Table C-extensions: merged isolated tests for join_inner (PR #19458), documented jointype integer mappings (PR #19468), and added direct C-slot tests for _column_mixins (PR #19806).
• XML iterparser: built structural 4-tuple boundary tests for the _iterparser C-extension and caught memory buffer behaviors (PR #19922).
• Type stubs: defined the static type boundaries for _iterparser.pyi (PR #20006, positional-only) and _convolve (PR #20093, dimensionality locking).
• Time parse extension: overhauled the _parse_times.c boundary test using strict type(parser) is np.ufunc assertions and byte arrays ("S24") (PR #19875).
• Quantity unit extraction: fixed unit extraction in structured_to_unstructured for StructuredQuantity (PR #19106).
• APE split expansion: contributed to the proof of concept for the APE split, expanding the implementation section in Clement's fork (PR #3) to define what the low-level tests look like.

Still open, waiting on review

• WCS dependency injection: refactored unit_list_proxy.c to cache the Unit class pointer, severing the PyImport overhead (PR #20072).
• Timeseries math engines: built exact synthetic recovery matrices for Box Least Squares (bls.c / _impl.pyx, PR #20151) and Lomb-Scargle (cython_impl.pyx, PR #20222).

The Developer Guide

Perhaps the most lasting piece of work I did was writing and merging testing_extensions.rst (PR #20239) into Astropy's core documentation. This document establishes the architectural standards for making low-level compiled modules testable. It outlines the rules for bypassing Python wrappers, strictly avoiding isinstance() at the compiled boundary, structural validation, safe exception handling using pytest.raises substring matching, and .pyi architecture. If anyone has any questions regarding the project I think this would serve as a good starting point (before the actual coding part of course).

What's left after I leave?

The work on astropy/erfa was pivoted to the external pyerfa repository, so testing that boundary remains a future task for whoever manages that standalone package.

Beyond that, the primary future work is maintenance. When a new Cython extension is updated or a C-boundary changes, these isolated tests will immediately flag it. The testing_extensions.rst guide exists specifically so future contributors don't have to stumble through the same trapdoors I did. Future contributors can use this suite as a blueprint to expand coverage into newer subpackages.

Sidenote: I plan on staying (long time). I'm not going to disappear just because GSoC is over but rather continue to contribute to this epic organization. I'll keep working, raising PRs, beefing with GitHub and most importantly, I'll help guide all those who stumble across Astropy just as Clément and Nathaniel were always there for me!!

What I actually took away from this

I came into GSoC thinking I was going to write tests. Somewhere along the way, I ended up learning how to actually reason about a mature codebase.

I learned that you have to read the C code before writing the Python test. I learned to question unexpected behavior rather than blindly cementing it into a test. I learned that CI is a judgmental gatekeeper, Git history tells a story, and numerical tolerances matter a whole lot more than just getting the test to pass. I learned that maybe breaking a huge problem into separate sections wasn't such a bad idea after all and it's this lesson that I'll be applying to my life as well.

Most importantly, I learned how maintainer-level review actually works. The feedback I got wasn't just "fix this bug," it was "let's structure this so it's statically typed, mathematically locked, and future-proof." And it really shaped my journey. I mean the Reem from January decided to open her very first PR and make huge changes in 4 different files, which ultimately was pretty overwhelming and she closed the PR. My point is that I have grown as a contributor and that makes me so very proud.

A HUGE thank you

All this is just because I had the BEST mentors EVER, Clément Robert (@neutrinoceros), and Nathaniel Starkman (@nstarman).

Clément, thank you for never once making me feel dumb for asking the same question twice, and for having the patience to let me actually sit with a hard bug instead of just handing me the fix. For taking time out from your busy schedule every week, just so that I could ask you all the questions I saved, and to go through my code, suggest changes and help me understand all the whys. You didn't just correct my code, you really changed how I think about coding. I came into this summer knowing how to write Python and I'm ending it knowing how to reason about the machine underneath.

Nathaniel, thank you for being just as generous with your time, and just as willing to slow down and explain something properly even when the answer must have felt obvious from where you were standing.

Having two mentors who both cared this much about getting it right, and not just getting it merged, made all the difference. And I couldn't be more grateful!

To the OpenAstronomy and Astropy community: thank you for letting me poke around the raw engine of your library.

Until next time,
Reem <3