Gaia DR3 EpochPhotometry Analyzer

Identifies variable astronomical objects in Gaia DR3 epoch photometry data using IRIS Embedded Python with a C shared library (gaiaV2.c) for maximum throughput. Built for InterSystems Employee Programming Challenge #1.

Challenge

Given 20 gzip-compressed CSV files from the Gaia DR3 epoch photometry archive (~380 MB compressed), find every astronomical source whose BP or RP flux changed by more than 100% over its observation period.

Each row is one source’s complete light curve stored as quoted array strings like "[1234.5,null,6789.0,...]". There are 48 columns per row; only 3 are needed.

For each qualifying source, output:

Invalid (null/NaN) flux values are ignored. Sources with non-positive minimum flux are excluded.

Solution

The entry point is an ObjectScript routine (RunScript.mac) that delegates all processing to a Python module (processV10.py) via IRIS Embedded Python. processV10.py is a thin ctypes wrapper that calls into a pre-compiled C shared library (gaiaV2.so) baked into the Docker image at /usr/local/lib/gaiaV2.so.

The C library (gaiaV2.c) handles the entire pipeline — gzip decompression, CSV parsing, flux computation, filtering, and output — using one pthread per file, all running in parallel.

Prerequisites

• git
• Docker Desktop

Build

Clone the repository and place the 20 Gaia EpochPhotometry .csv.gz files into data/in/:

git clone 
cd intersystems-challenge-GAIA-C

Start the IRIS container (this also compiles gaiaV2.so via the multi-stage Dockerfile):

docker-compose up --build -d

The Dockerfile compiles gaiaV2.c with gcc -O3 -march=native -funroll-loops in a Debian builder stage and copies the resulting .so to /usr/local/lib/gaiaV2.so in the IRIS image. The library is placed outside the volume-mount path so it is not shadowed when the repo directory is bind-mounted at runtime.

Run

Open an IRIS terminal:

docker compose exec iris iris session iris -U USER

Compile and run:

do $System.OBJ.Load("/home/irisowner/dev/https://github.com/isc-mdotti/intersystems-challenge-GAIA-C/blob/main/src/RunScript.mac","ck")
do ^RunScript

Output is written to data/out/results.csv.

How the evaluator checks your work

docker-compose exec iris iris session iris
USER>do ^RunScript

Architecture

Entry point — RunScript.mac

Run()
{
  Do ##class(%File).CreateDirectoryChain("/home/irisowner/dev/data/out")
  Set sys = ##class(%SYS.Python).Import("sys")
  Do sys.path."append"("/home/irisowner/dev/src")
Set startTime = $ZHOROLOG
Set proc = ##class(%SYS.Python).Import("processV10")
Do proc.run()
Set elapsed = $ZHOROLOG - startTime
Write "Elapsed time: ", elapsed, " seconds", !
}

Creates the output directory, appends src/ to the Python path, starts the $ZHOROLOG timer, then calls processV10.run() — the entire data pipeline — before stopping the timer.

Python wrapper — processV10.py

A 28-line ctypes shim. On first call it loads /usr/local/lib/gaiaV2.so and sets the run function signature:

_lib.run.restype  = ctypes.c_int
_lib.run.argtypes = [ctypes.c_char_p, ctypes.c_char_p]

Then calls lib.run(IN.encode(), OUT.encode()) where IN = /home/irisowner/dev/data/in and OUT = /home/irisowner/dev/data/out. The .so lives at /usr/local/lib rather than inside the repo so the Docker volume mount at /home/irisowner/dev cannot shadow it.

C library — gaiaV2.c

The core of the solution. Exported symbol: int run(const char *in_dir, const char *out_dir).

Pipeline per file (one pthread per .csv.gz):

1. Decompress — mmap(MAP_POPULATE) pre-faults the file into the page cache, then libdeflate decompresses the full gzip stream into a malloc’d buffer. MAP_POPULATE overlaps disk I/O with decompressor setup; libdeflate’s SIMD Huffman decoder is ~2× faster than zlib for single-stream inputs.
2. Skip metadata — the first ~365 lines beginning with # are ECSV metadata comments. The CSV header row follows. Both are skipped with a byte-scan loop.
3. Parse rows — for each data line, extract_fields pulls columns 1, 11, and 16 (source_id, bp_flux, rp_flux) in a single RFC-4180-compliant pass, treating quoted fields containing commas correctly. All other columns are discarded without storage.
4. Compute min/max — parse_float_array scans the bracketed array literal, skips null tokens, and accumulates min and max using fast_parse_double — a custom inline parser that avoids strtod’s locale overhead, giving ~10% per-file speedup.
5. Filter — rows are rejected unless both bp_min > 0 and rp_min > 0, and at least one band satisfies max >= min * 2 (equivalent to >100% change but avoids division for the ~99% of rows that don’t qualify).
6. Format — qualifying rows are formatted into a per-thread output buffer (outbuf) using snprintf. The fmt_double formatter strips trailing zeros to match Python’s repr-style output.

Main thread after join:

Assembles a single contiguous buffer: CSV header + each thread’s pre-formatted chunk. Writes the entire output in one write() loop — one syscall instead of one per row.

Largest-first dispatch: Files are sorted by st_size descending before thread launch. This puts the longest job on the critical path from the start, keeping all threads busy until the last file finishes and minimising tail latency.

Dockerfile

Two-stage build:

The .so and its libdeflate.so.0 runtime dependency are copied from builder into /usr/local/lib in the final image. No compiler toolchain is shipped in the IRIS layer.

Performance

Times measured with $ZHOROLOG on IRIS Community 2026.1.

Why the C version is faster

libdeflate vs polars gzip: libdeflate’s DEFLATE decoder uses SIMD Huffman decoding and processes the full stream in one shot — about 2× faster than zlib (which polars wraps) for single-stream inputs. With 20 files running in parallel this compounds.

mmap + MAP_POPULATE: Pre-faults all file pages into the page cache before decompression begins. For 20 concurrent threads this means 20 files are being read off disk in parallel at the OS level, overlapping I/O with decompressor initialization.

Per-thread output buffers: Each thread formats its own CSV rows into a private buffer. The main thread only concatenates. No mutex, no shared state, no fprintf-per-row serialisation.

fast_parse_double: Avoids strtod’s locale lookup and edge-case handling (hex floats, NaN/Inf literals) that never appear in this dataset. Accumulates digits into a uint64_t, then converts once to double.

Single write(): One write() syscall for the entire output file rather than one fprintf per row. Reduces kernel-crossing overhead and lets the OS schedule one large sequential write.

Profiling context (polars single-file benchmark)

read_csv:       0.386s  (5345 rows)   ← 94% of total
parse+minmax:   0.021s
filter:         0.002s
pct_change:     0.001s
total:          0.410s/file → ~1.0s with 20 threads

The bottleneck was always I/O + decompression, not computation. The C solution attacks the dominant cost directly.

Technologies

• InterSystems IRIS — runtime and ObjectScript entry point
• IRIS Embedded Python — bridge between ObjectScript and C via ctypes
• libdeflate — SIMD-accelerated gzip decompression
• pthreads — one thread per file for parallel I/O
• polars — installed but only used in earlier benchmark versions; not used in V10

Development

Start the IRIS container:

docker-compose up --build -d

Open an IRIS terminal:

docker compose exec iris iris session iris -U USER

Recompile RunScript.mac and run:

do $System.OBJ.Load("/home/irisowner/dev/https://github.com/isc-mdotti/intersystems-challenge-GAIA-C/blob/main/src/RunScript.mac","ck")
do ^RunScript

Rebuild from scratch (no Docker cache):

docker-compose down
docker system prune -f
docker-compose up --build -d

Note on .py vs .c changes: The repo directory is volume-mounted at /home/irisowner/dev. Python file edits (processV10.py) take effect immediately. Changes to gaiaV2.c require a full image rebuild (docker-compose up --build) because the .so is compiled at image-build time and placed outside the volume mount.

IRIS connection credentials (development only): user _SYSTEM, password SYS, port 52773.

Data Format

The input files are ECSV 1.0 (Enhanced CSV). Each begins with ~365 comment lines (#) describing the schema, then a column-name header row, then data rows starting at line 367.

Key columns:

For each source:

1. Extract all finite, non-null flux values from both arrays
2. Find min and max per band
3. Reject if bp_min <= 0 or rp_min <= 0
4. Pass if bp_max >= bp_min * 2 or rp_max >= rp_min * 2
5. Compute percentage_change = max((bp_max−bp_min)/bp_min, (rp_max−rp_min)/rp_min) × 100