Skip to main content

How To Get Maximum Performance

This page collects the performance model, implementation notes, thread tuning, and optimized-workspace guidance for imcombiners. Exact timings vary with shape, dtype, rejection density, stack depth, memory layout, CPU, and thread settings, so treat printed numbers as local measurements rather than universal constants. See 2-D Image Benchmarks for measured image-combination results.

The examples below run on each render. Their environment is printed here so the measured tables can be distinguished from the saved benchmark snapshots.

Measurement environment Value
Measured (UTC) 2026-11-11T17:47:08+00:00
Python 3.12.16
Platform Linux-7.0.0-1019-azure-x86_64-with-glibc2.43
imcombiners / reducers / NumPy 0.1.2 / 0.4.0 / 2.4.6
Rayon threads / sigma grain 4 / 1024

1. Implementation Tricks

imcombiners is built around image-stack reductions where both algorithmic work and memory traffic matter. The following implementation choices are intended to change runtime and memory traffic, not mathematical results. Here, “vector” is a Rust term for a contiguous 1-D array.

  • Fuse output-only rejection and final reduction when it is exact. IRAF avoids duplicate work when the rejection path already has the retained samples needed for the requested final value. imc applies that work-elimination principle through supported mean/median output-only Rust paths.
  • Axis-first public layout. (N, *spatial) lets np.stack(images, axis=0) copy each image into a contiguous slab. See Why The Stack Axis Is First for the layout tradeoff.
  • Prepared float workspaces. Public calls validate and promote dtypes at the boundary; hot loops can prepare a C-contiguous float32/float64 workspace once and then use validate=False on already-checked data.
  • Fused output-only kernels. With diagnostics=None, one supported rejector and an unweighted mean/median combine can run inside one Rust call when no growth is requested. Clipping still uses multiple iterations. These paths avoid a stack-sized rejection-mask result and a separate NaN-filled combine workspace; diagnostics=None alone does not guarantee fusion.
  • Selective order statistics instead of full sort. Medians and percentiles select the required ranks without sorting all samples. Rejection paths apply selection to the surviving samples.
  • Fused variance/mean only when both are requested. return_mean=True computes variance and mean together, avoiding a second axis reduction.
  • Scratch reuse. Supported rejection paths reuse temporary buffers for pixel samples, masks, statistics, and rank pairs. Parallel tasks have their own scratch buffers; serial loops reuse one set. Diagnostic paths retain output masks and diagnostic arrays as well.
  • Raw contiguous slice indexing. After validation, the Rust kernels require C-contiguous NumPy arrays. They view the stack as one flat memory slice and compute offsets directly, avoiding repeated multidimensional ndarray indexing and shape/stride bookkeeping in inner loops.
  • Unrolled reductions with separate accumulators. Sum, mean, and variance scans process several stack samples per loop and keep separate running sums/counts before merging them at the end. This reduces loop overhead and short dependency chains in the CPU pipeline. For axis-0 image stacks, the Rust loop advances a flat slice index by H * W to read arr[:, y, x] without repeated multidimensional indexing.
  • Branch-light min/max scans. NaN-aware stack min/max reductions use the same numeric min/max semantics as the desired nanmin/nanmax behavior, so the inner loop does not need a separate branch just to ask whether the current output is still NaN. This is not “Rust is faster by default”; it is a specific way to avoid an unnecessary branch in a very small hot loop.
  • CCD clip without a full gain copy. ccdclip, ccdclip_mask, and ccdclip_combine pass gain into the Rust noise model and compare residuals in the supplied data units. Per-image scales and zeros describe any prior normalization; no gain-scaled copy of the full stack is needed. The fused combine uses the supplied sample values.
  • Mask-aware finite handling. Non-finite samples are folded into masks once, so later inner loops can avoid redundant finite checks where correctness permits.
  • Squared-spread sigma clipping. sigma is the dimensionless multiplier of the measured spread. For sigclip, the Rust loop compares (x - center)^2 > sigma^2 * spread^2, using the appropriate multiplier for each tail. Mask-only and output-only paths can avoid the square root needed to return diagnostic std. CCD clipping instead computes each image’s noise sigma and scans sorted tails using normalized residuals; its reported std is a separate reference diagnostic.
  • MAD sigma without a full sort. sigclip(stdfunc="mad") is primarily a robustness option, not a reason to choose MAD for speed. Its implementation still avoids a full sort: it computes 1.4826 * median(abs(x - clip_cen)) using the same selective-median scratch buffer strategy as other median-style kernels.
  • Rank-pair final combine. Output-only minmax and pclip paths gather retained values from sorted rank pairs into reused scratch for the final mean or median.
  • Reusable pclip residual scores. The pclip implementation follows IRAF-style rank-offset clipping. During that calculation it builds a temporary vector of normalized distance scores used to decide which sorted values to keep when nkeep restoration is needed. That vector is reused from scratch memory instead of allocated fresh for every output pixel.

Why The Stack Axis Is First

imcombiners uses (N, *spatial) stacks, so a normal image cube is shaped (N, H, W) and the combine/reject axis is axis 0. This is deliberate.

For C-order arrays, (H, W, N) makes each per-pixel stack contiguous, while (N, H, W) makes each input image contiguous. Constructing either layout copies the full cube:

stack0 = np.stack(images, axis=0)   # shape (N, H, W)
stackn = np.stack(images, axis=-1)  # shape (H, W, N)

With axis=0, each image is copied into one contiguous output slab. With axis=-1, each image occupies a strided [..., i] plane. The public axis-0 convention matches image-by-image input and keeps masks, weights, zero/scale values, growth, and rejection diagnostics aligned with the same stack axis. Assess layout changes by timing construction and reduction together; this page does not measure an axis-last implementation.

np.moveaxis(stack, 0, -1) returns a view; it does not rearrange memory. Obtaining contiguous axis-last storage generally requires a full copy with np.ascontiguousarray(np.moveaxis(stack, 0, -1)).

2. API Choice

The shortest rule is:

  • Use Standard Combiner or ndcombine() described in tutorials when you only need the final image.
  • Use Chained Combiner or diagnostics="simple"/"full" when you need last_reject, mask_rej, low, upp, nit, or output_flags, plus sample_flags.
  • Use direct kernel calls (import imcombiners.kernels as imck) only for trusted hot loops where the input layout is already guaranteed.

Direct imck calls use the same Rust kernels as the supported public routes. Prepared calls with validate=False skip repeated input checks and workspace conversion. They also avoid Combiner construction and its default workspace copy. Compare imck.sigclip_combine(..., validate=False) with Combiner(..., copy=False, validate=False).combine(..., diagnostics=None) to measure that API overhead on your workload. Use reducers directly for standalone 1-D reductions.

Workspace Requirements

The trusted Rust paths require:

  1. C-contiguous layout — arr.flags['C_CONTIGUOUS'] must be True.
  2. Float workspace dtype — float32 for uint8, uint16, int16, and float32 inputs; float64 for int32 and float64 inputs.

This function shows one example:

import numpy as np

def prepare_workspace(stack):
    stack = np.asarray(stack)
    if stack.dtype in (np.dtype("int32"), np.dtype("float64")):
        dtype = np.float64
    elif stack.dtype in tuple(map(np.dtype, ("uint8", "uint16", "int16", "float32"))):
        dtype = np.float32
    else:
        raise TypeError(f"Unsupported input dtype: {stack.dtype}")
    return np.ascontiguousarray(stack, dtype=dtype)
  • int32 and float64 → float64
  • uint8, uint16, int16, and float32 → float32

This prepares the dtype and layout only. Check shape and all other kernel arguments before using validate=False.

3. Public API Baseline

The baseline for output-only timing is the default Standard Combiner: Combiner(...).combine(..., diagnostics=None). It keeps the public object API, returns only the final image, and can use the fused reject-and-combine Rust path for supported mean/median cases.

This benchmark compares that baseline with three nearby public routes:

  • SC-None: Standard Combiner, Combiner(...).combine(..., diagnostics=None), output-only reference row.
  • nd-None: compact ndcombine(..., diagnostics=None), output only.
  • CC: Chained Combiner, Combiner(...).reject(...).combine(...), retained rejection state.
  • nd-simple: compact ndcombine(..., diagnostics="simple")[0], final image plus rejection diagnostics.

Diagnostic rows ask for different work, so they are grouped separately and compared with Chained Combiner. Output-only rows are compared with Standard Combiner.

import numpy as np
import imcombiners as imc
import time
def time_call(func, warmups=5, repeat=100):
    result = None
    for _ in range(warmups):
        result = func()

    times = []
    for _ in range(repeat):
        t0 = time.perf_counter()
        result = func()
        times.append(time.perf_counter() - t0)
    return result, float(np.median(times))


def comparison_text(seconds, reference_seconds):
    ratio = reference_seconds / seconds
    if np.isclose(ratio, 1.0, rtol=0.05):
        return "Near-tie"
    if ratio > 1.0:
        return f"{ratio:.2f}x faster"
    return f"{1.0 / ratio:.2f}x slower"


def print_timing_table(rows, reference):
    print(f"{'name':<10} {'work':<12} {'ms':>9}  comparison")
    ref_seconds = rows[reference]["seconds"]
    for name, row in rows.items():
        seconds = row["seconds"]
        if name == reference:
            comparison = ""
        else:
            comparison = f"{comparison_text(seconds, ref_seconds)} vs {reference}"
        print(f"{name:<10} {row['work']:<12} {seconds * 1e3:9.1f}  {comparison}")
rng = np.random.default_rng(20250311)
stack = np.clip(
    rng.normal(1000.0, 30.0, size=(31, 512, 512)),
    0,
    np.iinfo(np.uint16).max,
).astype("uint16")

kw_sc = dict(combine="median", reject="sigclip", sigma=(3.0, 3.0), maxiters=5, ddof=0)
sigc = imc.SigClip(sigma=(3.0, 3.0), maxiters=5, ddof=0)

timings = {}

timings["SC-None"] = {}
timings["SC-None"]["result"], timings["SC-None"]["seconds"] = time_call(
    lambda: imc.Combiner(stack).combine("median", rejectors=sigc, diagnostics=None)
)
timings["SC-None"]["work"] = "output-only"

timings["nd-None"] = {}
timings["nd-None"]["result"], timings["nd-None"]["seconds"] = time_call(
    lambda: imc.ndcombine(stack, **kw_sc, diagnostics=None)
)
timings["nd-None"]["work"] = "output-only"

timings["CC"] = {}
timings["CC"]["result"], timings["CC"]["seconds"] = time_call(
    lambda: imc.Combiner(stack).reject(sigc).combine("median")
)
timings["CC"]["work"] = "diagnostic"

timings["nd-simple"] = {}
timings["nd-simple"]["result"], timings["nd-simple"]["seconds"] = time_call(
    lambda: imc.ndcombine(stack, **kw_sc, diagnostics="simple")[0]
)
timings["nd-simple"]["work"] = "diagnostic"

np.testing.assert_allclose(timings["SC-None"]["result"], timings["nd-None"]["result"], rtol=1e-5)
np.testing.assert_allclose(timings["SC-None"]["result"], timings["CC"]["result"], rtol=1e-5)
np.testing.assert_allclose(timings["SC-None"]["result"], timings["nd-simple"]["result"], rtol=1e-5)

out_diagnostic_one_call = timings["nd-simple"]["result"]
t_diagnostic_one_call = timings["nd-simple"]["seconds"]
out_output_only_one_call = timings["nd-None"]["result"]
t_output_only_one_call = timings["nd-None"]["seconds"]
out_output_only_combiner = timings["SC-None"]["result"]
t_output_only_combiner = timings["SC-None"]["seconds"]
out_diagnostic_chained = timings["CC"]["result"]
t_diagnostic_chained = timings["CC"]["seconds"]

print_timing_table(
    {k: timings[k] for k in ("SC-None", "nd-None")},
    reference="SC-None",
)
print()
print_timing_table(
    {k: timings[k] for k in ("CC", "nd-simple")},
    reference="CC",
)
name       work                ms  comparison
SC-None    output-only       79.4  
nd-None    output-only       69.7  1.14x faster vs SC-None

name       work                ms  comparison
CC         diagnostic       180.7  
nd-simple  diagnostic       176.8  Near-tie vs CC

nd-simple and CC retain rejection diagnostics; the output-only rows request only the final image and can use fused rejection/combine kernels. These are different workloads, so their timing difference includes diagnostic work. Combiner() also uses copy=True by default, making an owned workspace copy after validation. Use copy=False only when sharing the caller’s workspace is acceptable: later in-place pipeline operations may modify it.

4. Prepared Data & validate=False

Prepare the workspace once, then reuse it inside the hot loop. The baseline uses the original uint16 stack; prepared calls use a C-contiguous float32 workspace built before timing. These ratios combine conversion and validation savings, with an additional copy saving for SC-prep:

  • CC-prep: Chained Combiner with validate=False; default copy=True remains.
  • SC-prep: output-only Standard Combiner with validate=False and copy=False.

Preparation time is excluded. Include it when evaluating a single combination; these rows measure repeated use of an existing workspace, not validation alone.

stack_ws = prepare_workspace(stack)

def output_only_combiner_no_validation(arr):
    return imc.Combiner(arr, copy=False, validate=False).combine(
        "median",
        rejectors=sigc,
        diagnostics=None,
    )


prepared = {}

prepared["CC"] = timings["CC"]
prepared["CC-prep"] = {}
prepared["CC-prep"]["result"], prepared["CC-prep"]["seconds"] = time_call(
    lambda: (
        imc.Combiner(stack_ws, validate=False)
        .reject(sigc)
        .combine("median")
    )
)
prepared["CC-prep"]["work"] = "diagnostic"

prepared["SC-None"] = timings["SC-None"]
prepared["SC-prep"] = {}
prepared["SC-prep"]["result"], prepared["SC-prep"]["seconds"] = time_call(
    lambda: output_only_combiner_no_validation(stack_ws)
)
prepared["SC-prep"]["work"] = "output-only"

np.testing.assert_allclose(prepared["CC"]["result"], prepared["CC-prep"]["result"], rtol=1e-5)
np.testing.assert_allclose(prepared["SC-None"]["result"], prepared["SC-prep"]["result"], rtol=1e-5)

out_chained_prepared = prepared["CC-prep"]["result"]
t_chained_prepared = prepared["CC-prep"]["seconds"]
out_output_only_prepared = prepared["SC-prep"]["result"]
t_output_only_prepared = prepared["SC-prep"]["seconds"]

print_timing_table(
    {k: prepared[k] for k in ("CC", "CC-prep")},
    reference="CC",
)
print()
print_timing_table(
    {k: prepared[k] for k in ("SC-None", "SC-prep")},
    reference="SC-None",
)
name       work                ms  comparison
CC         diagnostic       180.7  
CC-prep    diagnostic       178.4  Near-tie vs CC

name       work                ms  comparison
SC-None    output-only       79.4  
SC-prep    output-only       66.1  1.20x faster vs SC-None

SC-prep performs the same output-only combination as SC-None, while reusing the caller’s prepared workspace. CC-prep still copies that workspace before retaining rejection state. The result checks above verify agreement for this input; the printed times show the measured cost of each complete call.

5. Stack Number (0-th Axis)

For a stack shaped (N, H, W), runtime includes API checks, any workspace conversion/copy, and clipping/reduction work. Conversion and copying scale with the full stack size. Clipping also depends on iterations and retained samples, so its time need not grow linearly with N.

The following comparison repeats the prepared-workspace experiment with five images. It measures the combined effect of preparation, validation, and copy settings at each stack depth.

rng_small = np.random.default_rng(20250311)
stack_small = np.clip(
    rng_small.normal(1000.0, 30.0, size=(5, 512, 512)),
    0,
    np.iinfo(np.uint16).max,
).astype("uint16")
stack_small_ws = prepare_workspace(stack_small)

out_output_only_small, t_output_only_small = time_call(
    lambda: imc.Combiner(stack_small).combine(
        "median",
        rejectors=sigc,
        diagnostics=None,
    )
)

out_output_only_small_prepared, t_output_only_small_prepared = time_call(
    lambda: output_only_combiner_no_validation(stack_small_ws)
)

np.testing.assert_allclose(
    out_output_only_small,
    out_output_only_small_prepared,
    atol=1e-5,
)

depth_rows = {
    "N5-SC": {
        "work": "output-only",
        "seconds": t_output_only_small,
    },
    "N5-prep": {
        "work": "output-only",
        "seconds": t_output_only_small_prepared,
    },
    "N31-SC": {
        "work": "output-only",
        "seconds": t_output_only_combiner,
    },
    "N31-prep": {
        "work": "output-only",
        "seconds": t_output_only_prepared,
    },
}

print_timing_table(
    {k: depth_rows[k] for k in ("N5-SC", "N5-prep")},
    reference="N5-SC",
)
print()
print_timing_table(
    {k: depth_rows[k] for k in ("N31-SC", "N31-prep")},
    reference="N31-SC",
)
name       work                ms  comparison
N5-SC      output-only       12.1  
N5-prep    output-only       11.2  1.09x faster vs N5-SC

name       work                ms  comparison
N31-SC     output-only       79.4  
N31-prep   output-only       66.1  1.20x faster vs N31-SC

Compare the N=5 and N=31 ratios for these measured inputs. Repeat close results under similar machine load before treating a small difference as a useful improvement.

6. Thread Tuning

TipStart with defaults

Measure the relevant operation before changing thread counts or work grains.

Thread counts and work grains control different parts of scheduling:

  • RAYON_NUM_THREADS — how many threads Rayon uses inside Rust kernels. Set it before starting Python; Rayon reads it once at thread-pool initialization.
  • reducers parallel grains — work grains for axis reductions and generic sigma clipping. Sigma clipping uses the axis_order_median key returned by rd.get_parallel_grains(). These settings do not control CCD, linear, or IRAF rank rejection. Start with the built-in defaults and benchmark changes on the target machine.

Set the Rayon thread count from the shell before starting Python:

export RAYON_NUM_THREADS=8
python your_pipeline.py

Inspect reducers settings for generic reducer workloads:

import reducers as rd

print(rd.get_num_threads())
print(rd.get_parallel_grains())

# Optional machine/workload-specific tuning:
# python -m reducers.autotuner

There is no separate imcombiners parallel-control API. Measure the operation that matters to your pipeline before changing thread counts or grains.

Thread Benchmark Script

benchmark_threads.py measures a single-thread baseline and a range of Rayon thread counts in fresh subprocesses, so every candidate gets a clean thread pool. By default it runs sigclip_median on both (5, 512, 512) and (31, 512, 512) stacks and reports separate recommendations:

uv run --extra bench python benchmarks/benchmark_threads.py

The script prints an environment table before the timing results, including Python, package versions, OS/kernel details, machine type, and logical CPU count. Keep that block with any copied benchmark table so results can be compared later.

The reported values are a starting point for the measured workload, not a package-wide optimum. To tune for your own shape or operation:

uv run --extra bench python benchmarks/benchmark_threads.py \
  --op median --dtype float32 --n 31 --height 1024 --width 1024

Pass multiple --n values to cover several stack sizes in one run:

uv run --extra bench python benchmarks/benchmark_threads.py --n 5 31 63
Note

Default thread-count candidates. On machines with <= 32 logical CPUs the script tests every parallel thread count: 2, 3, …, os.cpu_count(). The baseline uses RAYON_NUM_THREADS=1; the worker still calls the same Rayon-backed implementation. Passing --threads 1 adds that count to the candidate list too.

On machines with > 32 logical CPUs, os.cpu_count() may not match your scheduler allocation. Pass --threads with explicit candidates, for example based on SLURM_CPUS_PER_TASK, or --force-test to get the sparse {2, 4, 8, ..., ncpu} list anyway.

Note

The benchmark uses the active reducers grain settings for reductions and generic sigma clipping. CCD, linear, and IRAF rank rejection use IMC’s internal policy. The script changes thread counts without overriding either policy.

More threads does not always mean faster. Repeated runs can recommend different near-best values because timing varies with operation, dtype, number of images, output shape, CPU allocation, memory bandwidth, and machine load. The tuner treats median times within --tie-tolerance of the fastest as near-ties and recommends the lower thread count. For a more stable estimate, increase --repeats and --warmups:

uv run --extra bench python benchmarks/benchmark_threads.py --repeats 15 --warmups 3

If your pipeline also calls threaded NumPy/BLAS code, consider pinning OPENBLAS_NUM_THREADS, OMP_NUM_THREADS, MKL_NUM_THREADS, and VECLIB_MAXIMUM_THREADS to 1 while benchmarking to avoid measuring nested oversubscription.