Add Python interface, refactor structure, enhance CI, optimize performance, and add testing framework - #2
Open
wahln wants to merge 38 commits into
Open
Add Python interface, refactor structure, enhance CI, optimize performance, and add testing framework#2wahln wants to merge 38 commits into
wahln wants to merge 38 commits into
Conversation
Mechanical cleanups along code that runs once or more per particle step. Verified bit-identical against a serial reference run of the omc_dosxyz smoke test (WATER phantom, 2000 histories, nsplit=20). - howfar()/hownear() in both user codes decoded the region number with three integer divisions. Reusing the quotient of the first division brings that down to two, each of which the compiler pairs with its own remainder into a single instruction. - initHistory() in omc_matrad located the entry voxel with three linear scans over the boundary arrays, costing up to isize+jsize+ksize iterations per primary history. Replaced with a binary search, which additionally clamps instead of running past the end of bounds[] if the position is ever out of range. - Replaced pow() with integer exponents by explicit multiplication in the hot physics (rayleigh, pair, compton, photo, msdist, brems, moller, bhabha, annih, electron). Compilers normally fold the squares already, but the cubes and fourth powers are only folded under unsafe-math flags, so brems in particular was making real libm calls inside a rejection loop. Grouping is preserved with parentheses so the result is unchanged. Init-time call sites are left alone deliberately. - The RNG scale factor 2^-24 was a struct member, so setRandom() loaded it from thread-local storage on every call. It is now a macro. No measurable wall-clock change on the small water-phantom smoke test; the value here is removing guaranteed-avoidable work from the hot path and not depending on optimizer heuristics for it. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Two independent speedups, both leaving results bit-identical on a serial
reference run of the omc_dosxyz smoke test.
rayleigh() looked up its starting bin with
ibin = (int)rnno0*dwi;
where the cast binds to rnno0 alone. Since rnno0 < 1 there, ibin was
always 0, so the i_array acceleration table built by initRayleighData()
was never actually used and the search below it degenerated into a linear
scan of the whole 100 entry CDF from element 0 on every Rayleigh event.
Casting the product restores the intended O(1) lookup. Results do not
change here because the linear scan was still converging on the same bin,
only slowly.
On the build side:
- -fno-math-errno where the compiler supports it. GCC and Clang otherwise
have to assume errno is inspected after every sqrt() and so cannot emit
a bare sqrt instruction, which matters because uphi21/uphi32, mscat,
spinRejection and msdist call it on every step. This is not -ffast-math
and does not change any value.
- Link time optimization on optimized configurations, behind the new
OMPMC_LTO option (default ON, silently skipped if the toolchain cannot
do it). ausgab(), howfar() and hownear() live in the user code but are
called from the transport loop in ompmc.c, so without LTO each call is
an opaque cross translation unit call on the hottest path in the
program. Where LTO is on and the toolchain ships gcc-ar/gcc-ranlib,
CMAKE_AR and CMAKE_RANLIB are pointed at them, otherwise plain binutils
ar warns "plugin needed to handle lto object" on every link.
Measured on the smoke test (2000 histories, WATER phantom, serial,
GCC 13.1, -O3): 13.3 s before, 11.7 s after the Rayleigh fix, 10.4 s with
the build changes on top, so roughly 22 percent overall. Both omc_dosxyz
and the omc_matrad MEX file build clean with LTO.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
All three are of the same kind: an index or an assignment that does not
do what the surrounding code plainly intends.
mscat() computed the width of the tabulated u bin as
u = ums_array[base + k];
du = ums_array[base + k] - u;
reading the same element twice, so du was identically zero. That disabled
the sub-bin interpolation immediately below it and pinned the sampled
multiple scattering angle to the lower edge of whichever bin was drawn,
biasing the angular distribution. Corrected to read element k+1, which is
the bin width EGSnrc uses. k is at most MXU_MS-1 here -- ims_array is only
populated up to that index -- so k+1 stays inside the slab.
initRegions() in both user codes, when the global ecut is below the PEGS
ae for a medium, printed a warning about falling back to the PEGS value
and then never assigned it, leaving region.ecut[i] as whatever malloc
returned and using uninitialised memory as an electron transport cutoff.
The neighbouring pcut branch assigns correctly; this now matches it. The
warning text said "pcut" in the ecut branch, fixed too.
initRegions() also wrote region.rhof[0], pcut[0] and ecut[0] inside the
loop over regions when a voxel is vacuum, so every vacuum voxel clobbered
region 0 and kept its own entries uninitialised. Changed to index i.
Effect on results: the ecut and vacuum bugs only fire on configurations
that do not occur in the bundled test cases. The mscat fix does change
sampling. Checked with 400k histories on the WATER phantom: the central
axis depth dose agrees with the previous behaviour to within 1.6 percent
everywhere above 5 percent of dmax, i.e. within the statistical noise at
that history count, and dmax moves by one 0.5 cm voxel. A proper
validation against EGSnrc reference data is still worth doing.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
There was no test infrastructure at all, so this adds a CTest based one:
cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build
ctest --test-dir build
Two tests are registered. "unit" runs tests/test_ompmc.c, a dependency
free executable with a small assertion harness that links against
ompmc_core and covers the pieces that do not need the PEGS and cross
section data loaded:
- pwlfEval/pwlfInterval, including that the interval index stays
consistent with what the evaluator returns, since every cross section
and stopping power lookup depends on that
- setSpline/spline: reproduces the tabulated values at the nodes and
tracks a smooth function between them
- heap_sort: ordering plus the index permutation it returns
- the RANMAR generator: range, reproducibility for a seed, that
different seeds diverge, that every value lands exactly on the 2^-24
grid the integer state implies, and a five sigma check on the mean.
All of these cross the NRANDOM refill boundary.
- selectAzimuthalAngle: the unit circle invariant and both first moments
- the voxel geometry helpers
- kn_sigma0: positive and monotonically falling with energy
"smoke" runs the existing omc_dosxyz smoke test end to end from the
repository root, which is where the data folders it needs live.
To make the geometry helpers reachable from the tests, the region number
decode and the entry voxel search moved out of the two user codes into
static inline functions in omc_utilities.h: omcDecodeRegion() and
omcFindVoxelIndex(). Both were duplicated across omc_dosxyz and
omc_matrad, and omcDecodeRegion() had four copies. They stay inline
because howfar() and hownear() call the decode on every step.
The decode has a differential test against the three division form the
user codes used before this branch, over a whole grid, so the rewrite is
pinned rather than just asserted.
Tests are on by default and can be turned off with -DOMPMC_BUILD_TESTS=OFF.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
rayleigh() takes imed and applied it to pmax, but indexed i_array, fcum,
c_array, xgrid and b_array with no medium offset at all, even though each
is one flat allocation holding media.nmed slabs back to back. Every
Rayleigh event therefore sampled the momentum transfer distribution of
whichever medium happened to be first, scaled by the correct pmax for the
medium the photon was actually in.
The tables are not interchangeable. For the four media of TG119 loaded
from 700icru, bone's normalised CDF departs from air's by up to 0.109 in
absolute probability and 98 of its 100 i_array entries differ.
The effect on sampling is not a blur but a sign error. Rayleigh scattering
off a higher Z medium is broader, because the atomic form factor extends
to larger momentum transfer, so mean cos(theta) must come out lower for
bone than for air. Measured over 40000 samples per point:
30 keV 50 keV 100 keV
air bone air bone air bone
before 0.717 0.828 0.854 0.897 0.952 0.956
after 0.717 0.710 0.854 0.807 0.952 0.905
i.e. before the fix bone came out more forward peaked than air at every
energy, the opposite of the physics. Air, lung and tissue were nearly
indistinguishable from one another for the same reason, their only
variation coming from pmax.
Fixing the index also unmasks a latent out of bounds read. pmax can come
out marginally above 1 when lgle lands on the last energy interval, whose
interpolation coefficients are copied from the previous one and so
extrapolate; ibin would then reach RAYCDFSIZE-1 and ibin+1 read one past
the medium's slab, or past the whole allocation for the last medium. That
was masked while ibin was pinned at 0, so it is bounded here.
Also corrects initRayleighData() writing i_array[i*MXRAYFF + 0] where
every other access uses the RAYCDFSIZE stride. Harmless today because both
are 100, wrong if either changes.
Single medium results are unaffected: the WATER smoke test is bit
identical, since the offset is zero for medium 0.
Adds tests/test_media_data.c, registered with CTest as "media" and run
from the repository root because it needs the PEGS and cross section data.
It checks the per medium layout invariants of the Rayleigh tables and
asserts the bone versus air ordering above at all three energies. Verified
to fail on the unfixed code at every energy, not just to pass on the fixed
one.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The dose from one beamlet occupies a small fraction of the grid, but every accumulation, output and reset step walked all of it, and a Dij run does that nbeamlets times over. For 200 beamlets on a 2.25M voxel dose grid that is on the order of 90 GB of memory traffic spent almost entirely on zeros. ausgab() now records each voxel that receives energy in a list, and accumEndep(), accumulateResults(), the threshold scan, the sparse population and the per beamlet reset all iterate that list. This is equivalent rather than approximate: for an untouched voxel every one of those steps reduces to writing back the zero that is already there. Measured on omc_matrad with the bundled BOXPHANTOM fixture refined to 240^3 = 13.8M voxels, 25 beamlets, 10 batches, 32 threads, best of three: 15.23 s before, 10.60 s after, so 30 percent. The saving scales with gridsize times nbeamlets, and is invisible on the unrefined 48^3 fixture, which is about a thousand times too small to show it. The scoring arrays, ausgab() and accumEndep() move out of the two user codes into src/omc_score.c, since both had identical copies. Along the way: - resetBeamScore() clears accum_endep2 as well. The memset it replaces cleared only accum_endep, so one beamlet's squared sums were carried into the next and the variance output was wrong for every beamlet after the first. - score.ensrc += ein ran inside the parallel history loop unsynchronised in both user codes. It goes through scoreSource() now. omc_dosxyz actually reports that number. - The list is sorted before use, because a CSC column needs ascending row indices. test_omc_matrad_mex.m now checks that property, which nothing did before. Note on what is *not* here: combining deposits per voxel in thread local state before the atomic, to cut the contention on ausgab(), was implemented and then removed. Measured against the same baseline on omc_dosxyz it cost about ten percent, consistently, across an interleaved three way comparison of no change / list only / list plus combining. An electron does not stay in one voxel for long enough runs of steps to pay for touching thread local state on every call. The atomic contention is real but this is not the way to fix it. On omc_dosxyz, which has a single beamlet and so gains nothing from the list, the change measures somewhere between neutral and 5 percent slower; the machine's run to run spread is about the same size, so it cannot be resolved better than that here. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
region.pcut and region.ecut were arrays the size of the whole geometry,
but initRegions() filled them with max(global cut, the medium's PEGS
threshold) -- a value that depends only on the medium. Every voxel of a
given medium held an identical copy. They are now two MXMED+1 tables
inside struct Region, indexed by medium + 1 so VACUUM (-1) lands on slot
0, reached through the new regionPcut()/regionEcut() accessors.
Region data for a dose grid drops accordingly:
61x61x150 (WATER phantom) 15.6 MB -> 6.7 MB
144^3 (refined fixture) 83.6 MB -> 35.8 MB
240^3 387.1 MB -> 165.9 MB
Results are bit identical: a serial omc_dosxyz run on the WATER phantom
matches the previous commit exactly, which it must, since the values
looked up are the same ones.
On speed, honestly: I expected this to pay off in the transport loop and
could not measure that it does. Both omc_dosxyz at 32 threads and
omc_matrad on the 13.8M voxel grid land inside the run to run spread of
this machine. Two reasons it is smaller than it looks: voxel traversal is
spatially coherent rather than random, so the prefetcher was already
handling the extra arrays well, and the accessor turns one independent
load into a dependent pair (med[irl], then the table). It is kept for the
footprint, which is real, and for the two fixes below.
- The "global cut is below PEGS's value" warnings were emitted inside
the loop over regions, i.e. once per voxel. On a multi-million voxel
grid a single misconfigured cut meant millions of identical lines
through printf or mexPrintf. They are now printed once per medium.
- initRegions() validates the material index. It indexes the cut-off
tables now, so a bad index in a phantom file or in matRad's cubeMatIx
would read past them; previously it only produced strange physics.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
struct Stack was eleven parallel arrays -- iq, ir, e, x, y, z, u, v, w,
dnear, wt -- each MXSTACK long. Nothing ever sweeps a field across
particles; the transport code works on one particle at a time, indexed by
stack.np. So the eleven arrays become one array of struct Particle, 80
bytes each, and initStack() does one allocation instead of eleven.
What that changes about locality, per thread:
live particles cache lines, SoA cache lines, AoS
1 11 2
4 11 5
8 11 10
20 31 25
40 51 50
and 4 KB pages touched drops from 11, one per allocation, to 1. The case
that repeats most is a single particle read and written over and over
while it is transported, which is the top row.
Results are bit identical: a serial omc_dosxyz run on the WATER phantom
matches the previous commit exactly.
On speed: no measurable change, on either omc_dosxyz at 32 threads or
omc_matrad on the 13.8M voxel grid. Both land inside this machine's run to
run spread. The likely reason is that the stack is small and hot enough to
sit in L1 either way, and eleven pages is not enough to trouble the TLB.
Committed for the layout and because eleven allocations becoming one is
worth having on its own; not as a speedup.
That is now three predictions in a row -- ausgab write combining, the
region cut-offs, and this -- where the transport loop did not behave the
way the cache analysis said it would. Worth profiling before picking the
next one.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
- pair(): the angle rejection loop exited on rtest <= rejfactor OR theta < pi instead of AND (De Morgan slip porting EGSnrc's UNTIL), so the Motz-Olsen-Koch rejection function was effectively never applied and unphysical angles above pi could be accepted. Also restore EGSnrc's 1.02 safety factor on the estimated maximum. - initSpinData(): the last energy bin of q1ce_ms0 was filled from q1ce_ms1, corrupting the spin correction to the first MS moment in the topmost energy interval. - omc_matrad: the batch variance swapped the nbatch and nbatch-1 divisors, leaving a spurious mean^2/(nbatch*(nbatch-1)) term that floors every voxel's relative uncertainty near 10% at 10 batches. - omc_matrad initHistory(): the lower z clamp read ybounds, and the 2*DBL_MIN nudges were absorbed entirely next to any normal boundary coordinate; use nextafter() towards the opposite face instead. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
photon() re-evaluated the mean free path interval, its coefficients and the Rayleigh correction on every voxel crossing although they change only with the medium; cache them per medium, refreshed at the interaction point (which also stops a rare mix of the new medium's table with the old medium's interval when an interaction lands within SGMFP of a medium boundary). howfar()/hownear() spent two integer divisions per call decoding the region number; a thread-local irl -> (ix,iy,iz) memo in the shared utilities makes the decode free in steady state, with howfar() staging the neighbour's indices - known without division - whenever it truncates a step to a voxel face. howfar() also caches the reciprocal direction cosines so voxel marching multiplies instead of divides. The memo half is verified bit-identical to the previous code; the reciprocals change last-ulp rounding only (total dose difference shrank from 1.1% at 4k histories to 0.20% at 15k, i.e. pure statistics). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Photon flights are sampled against a majorant inverse mean free path built from per-medium maximum density ratios (region.rhof_max, filled by initRegions), jumping in one step regardless of voxel boundaries and accepting collisions with probability sigma(local)/sigma_max. Per-medium tables are evaluated once per photon line; the voxel-by-voxel howfar() marching is gone from photon transport, which now locates points with the new user-code callback regionIndex(x,y,z) - one multiplication per axis on the uniform grids both user codes use, binary search otherwise. The photon splitting VRT keeps its survivor Russian roulette but each sub-photon now delta-tracks an independent flight from the common start with a stratified first majorant hop; the previous scheme transported one photon to nsplit pre-sampled optical depths, which requires exactly the ray integration Woodcock tracking avoids. The stratified eta values partition (0,1), so the ensemble still samples the free-path distribution without bias. Validated against the ray-tracing code on a uniform water null test and on the heterogeneous PROSTATE phantom (4 seeds each, 200k histories): total dose agrees to +0.05% at 1.2 sigma and per-depth-slice deviations have rms d/sigma = 1.06 with no depth trend. Photon-only transport time drops about 1.4x; total runtime at clinical ecut is unchanged because electron transport dominates it. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
An electron of total energy below the new esave threshold whose residual CSDA range is shorter than the perpendicular distance to the closest region boundary cannot leave its voxel; deposit its remaining energy on the spot and stop transporting it. Rejected positrons still emit their annihilation photons. The technique's approximation is that bremsstrahlung the electron would have radiated below esave is absorbed locally, so esave should stay modest (~2 MeV). esave (total energy, MeV) comes from the input file, or from the optional mcOpt.esave field in the matRad user code; absent or zero disables the technique, so existing setups are unchanged. The dosxyz smoke test enables it to exercise the branch. Validated on the PROSTATE phantom (4 seeds, 200k histories each, esave 2.0): total dose -0.03% at -0.6 sigma, per-depth-slice deviations consistent with noise. Buys roughly 10-20% of total runtime on a 3 mm grid at ecut 0.521; the gain is limited by the boundary-crossing algorithm, whose electrons sit close to voxel faces where the residual range cannot undercut the perpendicular distance, and grows with voxel size and ecut. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
- The interaction routines (compton, pair, moller, bhabha, annih, rannih, brems) wrote into stack.p[np+1] without any bound check; they now verify the slot exists first, like the photon splitting loop always did. - Input keys are stored trimmed and matched exactly with strcmp; the substring match used before let a short key answer for a longer one depending on storage order. Lines without a 'key = value' form are skipped instead of handing strcpy a NULL. - The matRad user code fails with a clear message when a required MC options field is missing instead of passing NULL into the MATLAB API and crashing the session. The default spectrum file path no longer overwrites the buffer pointer it just allocated with a string literal. - initSpinData() checks every fread and validates the endianness marker. Doing so exposed that the binary spin file was opened in text mode, which on Windows translates CRLF byte pairs inside the records and stops at the first 0x1A byte: Windows builds have always read silently corrupted spin tables. Open with "rb". The measured dose impact of the corruption on a 6 MV prostate case is within statistical noise (total -0.01%), and timing is unchanged; Linux and macOS builds never translated and are unaffected. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
RANMAR was 24-bit and seeded per thread, so with dynamic OpenMP scheduling no two runs were comparable. Philox4x32-10 is a pure function of the 'rng seeds' key and a 128-bit counter whose high half holds the global history index, set by the new setRandomHistory() at the start of every history. Each history owns its own stream regardless of which thread runs it: repeated runs are now byte-identical, even across different thread counts (verified on the smoke case with 1x, 3x and default threads). Values are (word + 0.5)*2^-32, so the stream is finer than the old 2^-24 grid, mean-exact, and never exactly 0 or 1. The round function is checked against the published Random123 test vectors in the unit tests. Dose agreement with RANMAR verified on the water smoke case with eight independent seeds per generator: total dose -0.02% at -0.7 sigma, per-depth-slice d/sigma rms 1.11 over 150 slices. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Every per-energy table (esig/ededx/tmxs/... over MXEKE, gmfp/gbr/cohe/ pmax over MXGE) stored its slope and intercept in two separate arrays, so each pwlfEval() touched two cache lines. The pairs now live interleaved in one array -- entry i at [2*i] and [2*i+1] -- and pwlfEval() takes a single coefficient pointer. The per-medium mapping pairs (eke0/eke1, ge0/ge1) are indexed by medium only and stay separate. Purely a data-layout change: the smoke-case dose output is byte- identical to the previous commit (checkable at full precision now that runs are reproducible). Wall-time gain is below this machine's noise floor. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
A new electron of total energy below e_rr (MeV) survives with probability 1/f_rr and carries f_rr times its weight, otherwise it is removed without depositing, keeping the expectation value of the dose exact. The roulette is played once per electron, at its creation point: after the split-photon interaction compaction in photon() and after moller()/bhabha() in electron(). Positrons are left alone so their annihilation photons stay smooth. Input keys e_rr/f_rr (matRad: mcOpt fields of the same name); disabled by default. Verified unbiased on the water smoke case, 8 seeds per setting (e_rr=1.0 f=5: +0.011% at 0.3 sigma; e_rr=1.0 f=20: +0.021% at 0.5 sigma; e_rr=0.7 f=10: -0.019% at -0.6 sigma). The efficiency study argues for leaving it off in this regime: at 6 MV, 3 mm voxels, ecut 0.521 and nsplit 20 the technique saves only 4-8% wall time while raising the per-slice dose variance by 17-131%, an efficiency (1/variance/time) loss of 13-53% -- sub-MeV electrons are individually cheap here but carry a significant share of the local dose. The knob exists for regimes with a different balance (coarser grids, higher ecut). Also dedents the allocation lines the pwlf interleaving script left misaligned. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
MSVC only implements OpenMP 2.0, which for C compilation units does not accept declaring the loop variable inside the for statement of a parallel for (error C3015); the accumulateResults() loop now declares it beforehand like the history loops already do. Also makes the double-to-int truncation in mscat() explicit, which was the one C4244 warning MSVC raised in changed code. Verified with a local MSVC Release build: all four targets compile and the smoke case runs -- with dose output byte-identical to the MinGW build, and about twice as fast (interleaved: ~3.9 s vs ~8.8 s with gcc 13), so MSVC is the better choice for Windows binaries. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
MinGW GCC does not use the fast UCRT math routines: its bundled software log/exp are about 5x slower per call, and the transport samples -log(rng) for every photon flight segment. OMPMC_WITH_OPENLIBM fetches openlibm v0.8.7 (MIT licensed) at configure time, builds it as a static archive and links it ahead of the toolchain's math library on the core target, so every user code -- including the MEX file -- resolves libm calls from it without shipping an extra DLL. Measured on the smoke case with MinGW GCC 13: ~15% faster overall (libm is about a fifth of the runtime; the remaining gap to MSVC is emulated TLS and libgomp overhead, not math). Dose is physically unchanged: 8 seeds openlibm vs mingwex give a total difference of +0.0000% at 0.0 sigma and per-slice d/sigma rms 0.01 -- last-ulp rounding only, with run-to-run outputs still byte-identical. The option stays OFF by default (it needs CMake 3.25 and network access at configure time) and is switched on for the windows-x64-mingw CI job. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
MinGW GCC emulates thread-local storage on Windows: every access to a __thread or threadprivate variable is a __emutls_get_address() call, and the RNG state, the particle stack and the geometry memo are all thread local. llvm-mingw (clang targeting the same MinGW runtime) uses native TLS and, combined with OMPMC_WITH_OPENLIBM, measured about 2x faster than MinGW GCC on the smoke case -- within ~20% of MSVC, from a fully non-proprietary toolchain. Output is byte-identical to the gcc+openlibm build, so the existing physics validation carries over. The new windows-x64-llvm-mingw job downloads a pinned llvm-mingw release (which brings its own mingw32-make), points CMake at its clang explicitly so a gcc found on the image cannot shadow it, and ships libomp.dll -- the binaries' only non-system dependency -- with the artifacts. Both omc_dosxyz and the MEX file build and pass the local test suite with this toolchain; the MinGW GCC job stays for coverage of that toolchain. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Clang emits the same __kmpc_* calls that MATLAB's Intel runtime
(bin/win64/libiomp5md.dll) implements -- LLVM's libomp is a fork of
it -- and MATLAB ships the import library right next to the DLL. When
the MEX file is compiled with clang on Windows, link it against that
instead of libomp: one OpenMP runtime in the MATLAB process instead of
two KMP-family copies, which trip Intel's duplicate-runtime check
("OMP: Error #15"), and no libomp.dll to ship with the MEX file. The
one entry point newer clang emits that MATLAB's runtime predates,
__kmpc_dispatch_deinit, becomes a documented no-op shim -- runtimes
before LLVM 19 reclaim the dispatch buffers at thread teardown instead.
GCC keeps its libgomp: it emits GOMP_* calls and Intel's Windows
runtime has no GOMP compatibility layer; the two runtimes share no
symbols and coexist. Verified with MATLAB R2025b: the full MEX test
passes, twice in one session, with libiomp5md.dll the only OpenMP
import.
This also fixes a bug the openlibm change introduced for every MinGW
MEX build: openlibm's objects carry explicit dllexports, which switches
GNU-style linkers out of export-everything mode, leaving mexFunction
unexported and MATLAB reporting a missing gateway function. A .def file
now states the one export a MEX file needs, for GCC and clang alike;
both variants pass the MATLAB test again.
Also retires two stale claims in BUILDING.md: linking MATLAB's
libiomp5 does work from clang (the old warning predated an ABI-matched
compiler), and results no longer depend on the thread count since the
counter-based RNG.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
A process that aborts rather than failing a test took its own diagnosis with it: pytest redirects the file descriptors, so the abort message, any OpenMP runtime error and pytest's own faulthandler traceback all went nowhere. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
A test process killed by a signal leaves nothing behind: what it wrote to stdout is still in a buffer when it goes. macOS keeps a report with the termination reason and a backtrace of every thread. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The folders read out of the input table went into 128 byte arrays, but a value is up to BUFFER_SIZE-1 = 255 characters, and removeSpaces() copies until the source ends. A Python wheel under a temporary directory -- .../T/cibw-run-0wegr2r6/cp39-macosx_arm64/venv-test-arm64/lib/python3.9/ site-packages/ompmc/data/ is 148 characters -- ran twenty bytes past the end of the array. The overflow has always been there; only the macOS wheels showed it, because Apple's clang turns on -fstack-protector-strong by default and its __stack_chk_fail aborts through the system log rather than stderr, which is why the test process died without printing anything. The manylinux and Windows paths both stay under 128 characters. Size those buffers by what the input table can hold, and assemble the file paths with snprintf() so the length no longer has to be reasoned about. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
This pull request introduces comprehensive CI/CD automation for the project by adding new GitHub Actions workflows and Dependabot configuration. These changes ensure that builds, tests, and releases are automated across multiple platforms, and that dependencies for GitHub Actions are kept up to date.
Continuous Integration and Build Automation:
.github/workflows/build.ymlto automate building and testing the project on Windows, Linux (x64 and ARM), and macOS (Intel and ARM), including matrix builds for various toolchains and MATLAB/Octave integration. This workflow also collects and uploads build artifacts and generates code coverage reports.Python Wheels and Distribution:
.github/workflows/wheels.ymlto build, test, and publish binary Python wheels and source distributions for theomc_pythonuser code across all major platforms and Python ABIs, including automated upload to PyPI using trusted publishing.Dependency Management:
.github/dependabot.ymlto enable Dependabot for GitHub Actions, ensuring that workflow action versions are automatically updated with scheduled pull requests and grouped by update type.