# HEALPix Plotter - AI Coding Agent Instructions
## Project Overview
This is a Rust CLI application for visualizing HEALPix sky maps in Mollweide, Hammer, and Gnomonic projections. It reads astronomical data from FITS files and generates publication-quality plots with customizable colormaps, scaling, and output formats (PDF/PNG).
## Architecture & Key Components
### Core Modules (`src/`)
- **`plot.rs`**: Main plotting logic, handles Mollweide projection rendering and layout
- **`healpix.rs`**: HEALPix coordinate system utilities and data reading
- **`colormap.rs`**: Colormap management with 80+ built-in colormaps (matplotlib + custom)
- **`scale.rs`**: Data scaling transformations (linear, log, symlog, asinh, histogram)
- **`fits.rs`**: FITS file parsing and column extraction
- **`colorbar.rs`**: Colorbar rendering and tick formatting
- **`layout.rs`**: Figure layout calculations for borders and positioning
- **`projection.rs`**: Coordinate projection mathematics
- **`mollweide.rs`**: Mollweide-specific projection implementation
- **`render/`**: Output format rendering (PDF via Cairo, PNG via image crate)
### Data Flow
1. Parse CLI args → 2. Read FITS file metadata → 3. Load HEALPix data → 4. Apply scaling → 5. Project to Mollweide → 6. Render pixels → 7. Add colorbar/borders → 8. Save output
## Development Workflow
### Building & Running
```bash
cargo build # Compile with warnings
cargo run -- [args] # Run with CLI args
cargo run -- --help # See all options
```
### Testing
```bash
cargo test # Run unit tests (currently broken - see below)
```
**⚠️ Known Issues:**
- Tests fail due to `scale_value()` API changes (missing `Option<&HistogramScale>` parameter)
- Several unused imports and variables (run `cargo fix` to auto-clean)
- Dead code in `colorbar.rs` (`apply_lightness`, `sample_distortion`)
### Adding Colormaps
Colormaps are precomputed 256-entry RGB LUTs generated by `tools/generate_colormaps.py`:
```bash
cd tools && python3 generate_colormaps.py
```
This creates `.rs` files in `colormap/` that get `include!()`'d in `colormap.rs`.
## Key Patterns & Conventions
### Colormap Access
```rust
use map2fig::get_colormap;
let cmap = get_colormap("viridis"); // Case-insensitive, panics on invalid name
let available = available_colormaps(); // Vec of names
```
### Scaling Values
```rust
use map2fig::scale::{Scale, scale_value};
let normalized = scale_value(value, min, max, Scale::Log, NegMode::Unseen, None);
```
### Pixel Rendering
All output goes through the `PixelSink` trait:
- `PngSink`: Direct image buffer writing
- `CairoRasterSink`: Cairo context rendering
- `CairoImageSink`: Cairo image surface
### Error Handling
- FITS parsing: Uses `fitsrs` crate, panics on invalid files
- Invalid colormaps: Panics with list of available names
- Missing data: Filters out UNSEEN values, panics if no valid pixels
## CLI Usage Examples
### Basic plot
```bash
cargo run -- -f cosmoglobe.fits -o map.pdf
```
### Custom scaling and colors
```bash
cargo run -- -f data.fits --log --min 1e-6 --max 1e-3 -c plasma --gamma 0.8
```
### Histogram equalization
```bash
cargo run -- -f data.fits --hist --min 0.1 --max 0.9
```
## File Organization Notes
- `colormap/`: Auto-generated LUT files (do not edit manually)
- `tools/`: Python scripts for colormap generation
- `assets/`: Fonts and static assets
- Sample FITS files included for testing
## Common Tasks
- **Add new colormap**: Edit `tools/generate_colormaps.py`, run it, update `colormap.rs` includes
- **Fix scaling**: Modify `scale.rs` functions, update all `scale_value()` calls
- **Change projection**: Update `mollweide.rs` or `projection.rs` math
- **Add output format**: Implement new `PixelSink` in `render/`
## Dependencies
- `cdshealpix`: HEALPix coordinate math
- `fitsrs`: FITS file reading
- `cairo-rs`: PDF rendering
- `image`: PNG/image processing
- `clap`: CLI argument parsing
## ✅ SUCCESSFUL OPTIMIZATIONS (Completed)
**Tier 1: Direct Binary Reading for Float32 Columns (Feb 16, 2025):**
- **Achievement:** Eliminated FITS DataValue enum conversion by reading float32 binary data directly
- **Result: 3.4× speedup (71% improvement)** on large files (6.41s → 1.88s)
- **Method:** Parse FITS headers, find column offset, read binary float32 directly, convert to f64 in tight loop
- **File:** `src/fits.rs` new functions: `try_read_float32_column_fast()`, `parse_tform()`, `find_binary_table_data_offset()`
- **Impact:** Reclaimed 70% of execution time from FITS type conversion overhead; now GPU rendering is the new bottleneck
- **Backward compatible:** Falls back to slow path for non-float32 columns automatically
**Tier 1.1 & 1.2: Memory I/O and Percentile Optimization (Feb 16, 2025):**
- **Tier 1.1:** Eliminated `Vec<DataValue>` intermediate buffer in sparse FITS column extraction
- Result: 30-35% speedup by reducing random memory access patterns
- File: `src/fits.rs` lines 95-155
- **Tier 1.2:** Memory optimization for huge maps - streaming percentile computation
- **Problem:** nside=8192 allocated 45 GB (14.5× file size) due to double vector allocation
- **Root cause:** compute_mollweide_scale() + _plot_mollweide_pdf_impl() both allocated 806M pixel vectors
- **Solution:** Streaming percentile (sample 10M pixels max = 80 MB instead of 6.4 GB)
- **Result: 79% memory reduction (~45 GB → 9.4 GB) on nside=8192**
- **Bonus: 49% faster (39.2s → 20.08s)** due to eliminating double sort
- File: `src/plot/mollweide.rs` (new `compute_percentile_from_map()` + hybrid logic)
- Status: ✅ Tested on 25 MB, 193 MB, 577 MB, 3.1 GB files - all work perfectly with linear memory scaling
- **MmapFitsReader enabled for memory-mapped I/O**
- Result: 20-21% additional speedup by eliminating kernel memcpy overhead
- File: `src/fits.rs` lines 63-65 (1-line change)
- **Combined Results Summary (Tier 1 + 1.1 + 1.2):**
- **nside=8192 (3.1 GB file):**
- Memory: 45 GB → 9.4 GB (79% reduction, 5× improvement) ✅
- Speed: 39.2s → 20.08s (49% faster) ✅
- **Linear memory scaling now achieved:** 2-3× file size (vs 14.5× before)
- **Ready for production:** All file sizes tested and working perfectly
See `TIER1_OPTIMIZATION_SUCCESS.md` for initial Tier 1 analysis, and `TIER1_MEMORY_FIX.md` for Tier 1.2 memory optimization details.
**Tier 5: Prefetch Hints for Downsampling Inner Loop (Feb 17, 2026 - ✅ SUCCESSFUL):**
- **Bottleneck Identified:** Downsampling accounts for 75.93% of runtime (5.7s of 7.5s total)
- Root cause: 3.2 billion random memory accesses across 806M-pixel array, 82% CPU stall rate
- Memory system can't prefetch random access patterns; CPU pipeline starves waiting for cache misses
- **Solution Implemented:** x86_64 explicit prefetch hints (`_mm_prefetch`) in inner loop
- Prefetch 2 iterations ahead in downsampling loop to hide latency
- File: `src/healpix.rs` lines 1297-1330 (downgrade_healpix_map_xyf_parallel)
- x86_64 specific with fallback for other architectures
- **Result:** **+3.2% wall-clock improvement** (7.502s → 7.263s on nside=8192 map)
- Prefetch visible cost: 7.68% in perf call-graph (uses previously-idle CPU stall window)
- More stable: ±0.192s std dev vs baseline ±0.205s
- Validated with 5 real benchmark runs using Hyperfine
- **Key Insight:** Amdahl's Law applies: can't fix latency on 82% stall window without introducing _some_ cost
- **Documentation:** See `PREFETCH_OPTIMIZATION_RESULTS.md` and `DOWNSAMPLING_OPTIMIZATION_SESSION_FEB2026.md`
**Tier 5.1: Spatial Tiling for Cache Locality (Feb 18, 2026 - ❌ FAILED):**
- **Hypothesis:** Process pixels in 256×256 spatial tiles per HEALPix face to improve cache locality
- **Implementation:** Attempted tile-based iteration with ~3000 Rayon tasks
- **Result:** **-12.3% regression** (7.263s → 8.156s baseline 7.502s)
- **Root Causes:**
1. Task overhead from 3000 sub-tasks exceeds spatial grouping benefit
2. HEALPix NESTED indexing (peano/morton Z-order curve) defeats spatial locality assumption
3. Prefetch already solved the latency problem; gains were imaginary
4. Amdahl's Law in reverse: once one bottleneck fixed, others become proportionally larger
- **Revert:** `git checkout src/healpix.rs` successfully restored prefetch-only version
- **Lesson Learned:** Don't attempt iteration reorganization on already well-optimized loops; measure first
- **Documentation:** See `TILING_OPTIMIZATION_FAILURE_ANALYSIS.md`
### Remaining Optimization Opportunities
**Hard Limits:**
- Theoretical minimum bandwidth-limited time: 3.1 GB ÷ 9.1 GB/s = 0.34 seconds
- Current: 7.263 seconds (prefetch optimized)
- Remaining potential: ~95% of wall time (but heavily constrained by algorithm)
**Priority 1 (5-10× speedup potential):**
- **GPU Acceleration:** Downsampling is embarrassingly parallel (3.2B ops, no dependencies)
- Approach: CUDA or HIP implementation of downgrade_healpix_map functions
- Difficulty: HIGH (new toolchain, CUDA SDK dependency)
- Benefit: Likely 5-10× speedup based on GPU capability
**Secondary Options (Difficult, Low ROI After Prefetch):**
- **Vectorize Mollweide math with SIMD:** 15-25% theoretical, but trigonometric ops already LLVM-optimized; measured gain would be <5%
- **Ring-ordered processing:** Would break NESTED semantics, requires flag/documentation
- **Multi-socket systems:** Rayon already parallelizes well; current bottleneck is bandwidth-limited (per-thread, not total)
**⛔ Failed Approaches - Do Not Retry:**
- Tier 3 (original): Downgrade-during-parsing (25% slower due to per-pixel coordinate conversion overhead)
- Tier 5.1: Spatial tiling (12% regression due to task overhead, prefetch already solved latency)
---
## ⛔ KNOWN FAILED OPTIMIZATIONS (Do Not Retry)
**Tier 3: Downgrade-During-Parsing (Feb 2025):** Attempted to fuse downsampling into FITS load phase to avoid 50M→12M vector allocation. **RESULT: 25% SLOWER (6.41s → 8.04s)** due to exponentially expensive per-pixel coordinate conversions (pix2ang_nest + ang2pix). Added ~50 CPU cycles per pixel × 50M pixels = 2.5B cycle overhead. Memory allocation savings (~6%) were trivial vs transcendental math costs. **Amdahl's Law lesson:** Cannot optimize 6% of total time by adding work to 39% of total time.
See `TIER3_OPTIMIZATION_FAILURE_ANALYSIS.md` for detailed analysis. **DO NOT retry downgrade-during-parsing.** Focus instead on Tier 1 (direct column reading, 30-40% gain).
**F32 Precision Reduction (Feb 15, 2026):** Attempted to speed up math by casting f64→f32→f64 or using native f32 arithmetic. **RESULT: Both approaches were SLOWER by 2-3.7%** due to conversion overhead exceeding any math speedup. Math is only 11.8% of CPU time and is already well-optimized by LLVM. The real bottleneck is Mollweide projection algorithm (77.5%) and Cairo rasterization (3.57× slower than PNG).
See `docs/F32_OPTIMIZATION_RESULTS.md` for full analysis. **DO NOT attempt precision reduction again.**
**Tier 5.1: Spatial Tiling for Cache Locality (Feb 18, 2026):** Attempted to improve cache performance by processing pixels in 256×256 spatial tiles per HEALPix face. **RESULT: -12.3% REGRESSION (7.263s → 8.156s)** despite hypothesis of better cache locality. Root causes:
1. Task overhead: 3000 sub-tasks added scheduler overhead that exceeded any data locality benefit
2. HEALPix Morton curve defeats spatial grouping: NESTED indexing follows Z-order curve, not rectangular spatial order
3. Wrong problem to solve: Prefetch already eliminated the latency bottleneck; tiling tried to solve a non-existent problem
4. Amdahl's Law in reverse: Fixing 75% bottleneck makes remaining 25% proportionally larger; attempting to optimize wrong place causes regression
**Key Lesson:** After fixing primary bottleneck with prefetch (7.68% visible cost, 3.2% net gain), further iteration reorganization backfires. Memory bandwidth, not cache locality, is the actual constraint for random-access workloads.
See `TILING_OPTIMIZATION_FAILURE_ANALYSIS.md` for detailed root cause analysis and `DOWNSAMPLING_OPTIMIZATION_SESSION_FEB2026.md` for full session context. **DO NOT retry spatial reorganization on this workload.**</content>
<parameter name="filePath">/home/dwatts/projects/map2fig/.github/copilot-instructions.md