Unit-Cell Electron-Density C++ Benchmark

Benchmark UnitCell.zDensity_G for the original NumPy implementation and the C++ kernel. The C++ kernel is intentionally the only accelerated path for this calculation; special electron-density callback models retain the Python/SciPy implementation.

Timings exclude one warmup call and measure the standard atomic-density expression only.

Setup

Start Jupyter from benchmarks/ with an installed orgui package. The benchmark uses a synthetic unit cell with two coherent domains, so it does not depend on external files.

import statistics
import time

import numpy as np
from matplotlib import pyplot as plt

from orgui import __version__, get_build_config
from orgui.datautils.xrayutils import CTRcalc, CTRuc


print("orGUI version", __version__)
print(get_build_config())
print(f"C++ density acceleration available: {CTRuc.HAS_CPP_ACCEL}")

orGUI version 1.5.1.dev72+g06d7563e4.d20260724
{'native_optimization': False, 'target_cpu': 'generic', 'allowed_instruction_set': [], 'native_arguments': '', 'cpp_compiler_id': 'msvc', 'cpp_compiler_version': '19.43.34810', 'host_system': 'windows', 'host_cpu_family': 'x86_64', 'host_cpu': 'x86_64', 'available': True}
C++ density acceleration available: True

Benchmark Helpers

The profile coordinate z is in Angstrom. h and k are reference-frame reciprocal coordinates in r.l.u.; the result is complex electron density in electrons per Angstrom cubed.

def make_cell():
    cell = CTRcalc.UnitCell([3.0, 4.0, 5.0], [90.0, 90.0, 90.0])
    cell.addAtom("C", [0.2, 0.3, 0.4], 0.12, 0.18, 0.8)
    cell.addAtom("O", [0.7, 0.6, 0.1], 0.21, 0.24, 0.5)
    cell.setEnergy(10000.0)
    shifted_domain = np.vstack((np.identity(3).T, [0.1, -0.2, 0.15])).T
    cell.coherentDomainMatrix = [cell.coherentDomainMatrix[0], shifted_domain]
    cell.coherentDomainOccupancy = [0.65, 0.35]
    return cell


def time_call(func, repeats=9, warmups=1):
    for _ in range(warmups):
        func()
    timings = []
    for _ in range(repeats):
        start = time.perf_counter()
        func()
        timings.append(time.perf_counter() - start)
    return statistics.median(timings)


def format_seconds(value):
    if value < 1e-3:
        return f"{value * 1e6:.1f} us"
    return f"{value * 1e3:.2f} ms"

Correctness Check

Before timing, verify that the C++ kernel matches the NumPy reference for a multi-domain cell.

if not CTRuc.HAS_CPP_ACCEL:
    raise RuntimeError("Build or install orGUI with the C++ extension to run this benchmark.")

cell = make_cell()
z = np.linspace(-2.0, 7.0, 2_000)
CTRuc.set_accel_backend("numpy")
reference = cell.zDensity_G(z, 0.37, -0.22)
CTRuc.set_accel_backend("cpp")
accelerated = cell.zDensity_G(z, 0.37, -0.22)
np.testing.assert_allclose(accelerated, reference, rtol=1e-12, atol=1e-12)
print("C++ zDensity_G matches NumPy for the multi-domain test cell.")

C++ zDensity_G matches NumPy for the multi-domain test cell.

Benchmark Electron-Density Profiles

Use the median runtime across repeated evaluations. The profile lengths span interactive plotting through dense optical-profile grids.

cell = make_cell()
results = []
for points in (200, 2_000, 20_000, 200_000):
    z = np.linspace(-2.0, 7.0, points)
    timings = {}
    for backend in ("numpy", "cpp"):
        CTRuc.set_accel_backend(backend)
        timings[backend] = time_call(lambda: cell.zDensity_G(z, 0.37, -0.22))
    results.append({
        "points": points,
        "numpy_s": timings["numpy"],
        "cpp_s": timings["cpp"],
        "speedup": timings["numpy"] / timings["cpp"],
    })

print(f"{'points':>10} {'numpy median':>16} {'cpp median':>16} {'cpp speedup':>13}")
for row in results:
    print(
        f"{row['points']:10d} {format_seconds(row['numpy_s']):>16} "
        f"{format_seconds(row['cpp_s']):>16} {row['speedup']:12.2f}x"
    )

    points     numpy median       cpp median   cpp speedup
       200         518.0 us         228.0 us         2.27x
      2000         866.9 us         754.1 us         1.15x
     20000          5.51 ms          4.42 ms         1.25x
    200000        121.83 ms         42.75 ms         2.85x

Benchmark Results

The logarithmic scale shows how the native kernel changes runtime as the number of sampled profile points grows.

labels = [f"{row['points']:,}" for row in results]
x = np.arange(len(results))
width = 0.36
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.bar(x - width / 2, [row['numpy_s'] for row in results], width, label="NumPy")
ax.bar(x + width / 2, [row['cpp_s'] for row in results], width, label="C++")
ax.set_yscale("log")
ax.set_xticks(x, labels)
ax.set_xlabel("z profile points")
ax.set_ylabel("Median runtime [s]")
ax.grid(axis="y", which="both", alpha=0.25)
ax.legend()
fig.tight_layout()