Add geometry importance sampling and Russian roulette to egs_kerma - #1440
Open
mainegra wants to merge 22 commits into
Open
Add geometry importance sampling and Russian roulette to egs_kerma#1440mainegra wants to merge 22 commits into
mainegra wants to merge 22 commits into
Conversation
Implements geometry-based variance reduction for deep-penetration
calculations. Each region carries an importance I_r (default 1.0),
specified via 'region importances' inside the 'calculation geometry'
input block.
At every region crossing (detected in ausgab AfterTransport):
- Forward crossing (I_new > I_old): split into N copies, weight /= N,
where N is sampled stochastically so that <N> = I_new/I_old
- Backward crossing (I_new < I_old): Russian roulette with survival
probability p = I_new/I_old; survivor weight *= I_old/I_new
The FD estimator remains unbiased: FD scoring fires from the
pre-crossing position (in selectPhotonMFP) before any split or
roulette is applied.
Example input:
:start calculation geometry:
geometry name = iron_sphere
region importances = 1 1 2 2 4 4 8 8 ...
:stop calculation geometry:
Fix: AfterTransport was not enabled in setAusgabCall, so the splitting
and roulette code in ausgab() would never fire. Now enabled conditionally
when any region importance differs from 1.0.
Diagnostic output at the end of initScoring() reports whether importance
sampling is active and, if so, the importance range per geometry:
Geometry importance sampling: ON
iron_sphere 80 regions, importances [1.0, 64.0]
Adds a compact triplet-based input style for per-region importances: importance region ranges = ir_start ir_end I [ir_start ir_end I ...] Triplets are applied in order (later ones override earlier for overlapping ranges), making structured importance tables (e.g. doubling every N mfp) readable without enumerating all geometry regions. The legacy 'region importances' flat list is still accepted as a fallback for backward compatibility. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
When the importance ratio is large (e.g. exp(5) ~ 148 per step in the coarse section of the iron-sphere geometry), nsplit can exceed MXSTACK and abort the simulation. Instead, cap nsplit to the number of free slots before pushing copies. The estimator remains unbiased because the particle weight is divided by the actual nsplit used; only variance increases when the cap is active. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
$SELECT-PHOTON-MFP (selectPhotonMFP) is called once per photon MFP at
:PNEWENERGY:, not at each boundary crossing. This left prev_ir_imp
stale for the 2nd, 3rd, ... crossings within a single MFP.
Without this fix the importance ratio for each crossing was computed
relative to the region from the START of the MFP, not the immediately
preceding region. This caused two wrong behaviours mid-MFP:
- Spurious splitting within same-importance zones (ratio >1 when it
should be 1)
- Splitting instead of Russian roulette on inward crossings (e.g.
I=4 to I=2 gave ratio=2 split instead of ratio=0.5 roulette)
The inward-split case creates extra copies of primary photons (latch=0)
at lower-importance zones, skewing K_pri at shallow depths and changing
the apparent buildup factor.
Fix: update prev_ir_imp = ir at the end of every AfterTransport call
for photons, so each crossing is evaluated against its true predecessor.
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Geometry importance sampling must not act on primary photons (latch==0). The FD estimator scores K_pri from the current position to all scoring shells in a single call inside selectPhotonMFP(). When IS splits a primary photon, EGSnrc restarts transport for each copy, triggering another selectPhotonMFP() -> scoreInCV() call at the split point. This adds a spurious K_pri term (w * FD(p_split)) on top of the already-correct term scored at the photon origin (w * FD(p0)), causing K_pri to exceed the number of source photons — physically impossible. IS is now applied only to secondary photons (latch != 0), for which scoreInCV() does not score K_pri (!latch gate).
latch_ was a copy of the_stack->latch[np], which latch already holds from the AfterTransport block above.
Variance reduction settings are now reported in describeSimulation() after the geometry and calculation details, where all simulation description output belongs. initScoring() retains only the functional setAusgabCall(AfterTransport) logic.
When a secondary photon is split at a boundary, N-1 copies are pushed onto the stack. When each copy restarts transport it enters selectPhotonMFP() and calls scoreInCV(), scoring the FD contribution from the split point pb onward. But the original photon already called scoreInCV() from its MFP start ps, which fully accounts for the expected contribution from ps onward (including pb onward). The copies' calls are pure double-counting, overestimating K_total and inflating B. This bias was hidden in the original IS code because spurious within-zone splits (from stale prev_ir_imp) saturated the stack, capping nsplit to 0 actual copies. After the prev_ir_imp staleness fix the real per-zone splits create actual copies and the double-counting appears. Fix: mark each IS split copy with IS_COPY_FLAG (bit 30 of latch). selectPhotonMFP() clears the flag immediately and skips scoreInCV() for that one restart. All subsequent MFPs (after physical interactions) score normally — those represent new physics, not double-counting.
When a secondary photon is split at a boundary, N-1 copies are pushed onto the stack. When each copy restarts transport it enters selectPhotonMFP() and calls scoreInCV(), scoring the FD contribution from the split point pb onward. But the original photon already called scoreInCV() from its MFP start ps, which fully accounts for the expected contribution from ps onward (including pb onward). The copies' calls are pure double-counting, overestimating K_total and inflating B. This bias was hidden in the original IS code because spurious within-zone splits (from stale prev_ir_imp) saturated the stack, capping nsplit to 0 actual copies. After the prev_ir_imp staleness fix the real per-zone splits create actual copies and the double-counting appears. Fix: mark each IS split copy with IS_COPY_FLAG (bit 30 of latch). selectPhotonMFP() clears the flag immediately and skips scoreInCV() for that one restart. All subsequent MFPs (after physical interactions) score normally — those represent new physics, not double-counting.
FD is a VR technique and belongs in the variance reduction output block alongside geometry importance sampling. Prints ON/OFF with a per-geometry mapping of calculation geometry -> FD geometry. Co-Authored-By: Claude Sonnet 4.6
Most input files don't define any ausgab objects, which caused EGS_ObjectFactory::createObjects() to print a warning that confused users. Guard createAusgabObjects() with a non-destructive pre-check so the warning is only shown for callers where a missing section is genuinely unexpected (e.g. source definitions).
… count When the stack is too full to accommodate all nsplit copies, the previous code capped nsplit to the available slots and then divided the weight by that smaller number. A photon intended to carry weight w/148 ended up with weight w/10 if only 10 slots were free, producing 14.8× over-weighted outliers that dominated the variance in deep-penetration runs. Introduce n_actual = min(nsplit, avail) for the loop bound while keeping nsplit as the weight divisor. Copies that don't fit are simply dropped; this introduces a small downward bias only when the stack is genuinely full (when deep-penetration results are already unreliable), but eliminates the catastrophic variance spikes. Co-Authored-By: Claude Sonnet 4.6
…pies When avail < 2 (stack too full to push any copies), the weight reduction to w/nsplit was previously skipped, leaving the original particle with its full weight w — effectively un-split and massively overweighted relative to every other particle at that depth. Move the weight division outside the n_actual > 1 guard so it is always applied. IS_COPY_FLAG is kept inside the guard: with no copies there is no LIFO re-entry and the flag must not be set (it would suppress the next legitimate scoreInCV call instead). Co-Authored-By: Claude Sonnet 4.6
When the stack is too full to accommodate all nsplit copies, record the
event. Three counters flow through the full parallel state machinery
(storeState/readData/addState/resetCounter) so combined runs report
correctly: n_split_events and n_cap_events sum across workers;
max_stack_needed (= max(np + nsplit) over all capping events) takes
the per-worker maximum.
outputResults() prints:
IS splitting events: 1234567 total, 42 capped (0.003%)
Stack cap: largest deficit = 138 slots => recommend MXSTACK >= 2000138 (current: 2000000)
If no capping occurred the second line is omitted. If IS is not active
the entire block is suppressed.
Co-Authored-By: Claude Sonnet 4.6
Replace the scattered raw-format egsInformation calls in outputResults() with a single boxed section matching the ==...== style used in describeSimulation(). Items are grouped under Source particles, Scoring photons, and Importance sampling (IS) subheadings with aligned colons. Energy averages are computed before printing (using has_* guards), and the old "last case = %lld fluence = %g" debug line is replaced by proper labelled fields. Co-Authored-By: Claude Sonnet 4.6
The "Source fluence / particles" line now distinguishes between sources that return actual fluence per cm² (parallel beam, collimated) and those that return only a particle count (isotropic, beam, point, phase-space): if getFluence() == current_case the label reads "Source particles" with no units; otherwise "Source fluence" with cm^-2.
describeSimulation(): remove the per-geometry FD announcement (---> Scoring using forced detection) from the calculation geometry loop. The information is already present in the Variance reduction section (FD ON/OFF with geom -> fd_geom pairs); having it in both places is redundant. Co-Authored-By: Claude Sonnet 4.6
When the stack is totally saturated (avail = 0, n_actual = 0), the previous code still divided the particle weight by nsplit without creating any copies. For aggressive IS schemes (exponential, factor ~148 per 5 mfp), this rapidly drove particle weights to zero, leaving 2M near-zero-weight particles clogging the stack and doing expensive FD ray-traces for negligible statistical contribution. Replace the unconditional weight division with Russian roulette when n_actual == 0: kill with probability (nsplit-1)/nsplit, survive with probability 1/nsplit at original weight. This is unbiased — the expected post-RR weight equals w/nsplit — and immediately clears the stack pressure instead of accumulating zero-weight particles. Survivors retain their weight and will be split properly once room opens up. The n_actual > 0 path (partial capping) is unchanged. Co-Authored-By: Claude Sonnet 4.6
In scoreInCV(), the fallback branch for zero-attenuation regions
(sigma == 0, e.g. vacuum) returned exp_Att = 1.0, dropping the
geometric path-length factor t_sc[i] (or t_sc_tot for the total).
The correct limit of exp_Lambda*(1-exp(-mu*t))/mu as mu->0 is
exp_Lambda*t, so:
per-region: exp_Lambda*(1-exp_CV)/mu_cv → exp_Lambda*t_sc[i]
total: exp_Lambda_to_CV*(1-exp_CV)/mu_cv → exp_Lambda_to_CV*t_sc_tot
The same exp_Att drives weightedExpAtt, so the fluence FD estimator
is also corrected. Non-vacuum regions (sigma > 0) are unaffected.
Co-Authored-By: Claude Sonnet 4.6
mainegra
requested review from
ftessier and
rtownson
and removed request for
a team
July 29, 2026 18:26
Without this fix, geometry importance sampling (IS) silently did nothing when `score primaries = no` was the default: the latch-flagging block was gated on `if (score_primaries)`, so all photons kept latch=0, and the IS guard `if (latch && ...)` in AfterTransport was never true. Promote `imp_active` from a local variable to a class member (initialised in the constructor, set in initScoring). Condition latch-flagging on `score_primaries || imp_active` so that IS and primary-kerma scoring are fully independent of each other. Co-Authored-By: Claude Sonnet 4.6
…ight When the importance-sampling block splits a photon at a forward boundary crossing, it pushes n_actual-1 copies onto the stack and sets IS_COPY_FLAG on all of them and on the original. The flag suppresses one scoreInCV() call, so that a copy re-entering :PNEWENERGY: does not book a forced-detection score at a boundary, which is not a collision site. That is correct for every copy but one. PHOTON does not re-read np after $AUSCALL($TRANAUSA) (egsnrc.mortran:6526); it falls straight through to the rest of the :PTRANS: loop. The last copy pushed is therefore the particle that continues the original's in-flight free path. It never passes through :PNEWENERGY:, never calls selectPhotonMFP(), and so never consumes its flag at the split boundary. The flag survives to that particle's first genuine post-interaction free path and suppresses a legitimate score there. The result is one lost FD score per split event: a buildup-factor deficit linear in the splitting rate. Two importance maps differing only in step size (e^5 vs e^2.5) disagreed by 10% at 15 mfp, and the deficits scaled as 6.86 against a splitting-rate ratio of 6.83. Splitting and Russian roulette are unbiased by construction, so a map-dependent answer is necessarily a bug. Clearing the flag on the top of stack after the push loop removes the map dependence: the two maps move from 12.8 sigma apart to 1.7 sigma, and both agree with two independent calculations that use no importance sampling at all (30.671 and 30.217 against 30.304 and 30.328). The original at index np keeps its flag. It is a leftover that will be popped fresh at :PNEWENERGY: and would otherwise book the spurious boundary score the flag exists to prevent. Co-Authored-By: Claude Sonnet 5 <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.
Summary
egs_kerma: forward boundarycrossings stochastically split photons (⟨N⟩ = I_new/I_old); backward
crossings apply Russian roulette (p = I_new/I_old).
primary photon exclusion (latch==0 → skip IS),
prev_ir_impupdatedunconditionally at AfterTransport, IS_COPY_FLAG (bit 30 of latch) on both
original and copies to suppress FD re-scoring at split boundaries, and
nsplit/n_actual separation (weight always divided by nsplit; when n_actual=0,
RR with survival probability 1/nsplit).
from the optical-depth integral, producing ~1.5% bias at shallow depths.
Test plan
BUF_Fe_1MeV_xoshiro256pp_IS-exp.egsinp; verifyconverged BUF values are similar to published ANSI/ANS iron-sphere data at shallow depths.
as primary-kerma branch at shallow depths.
BUF_Fe_1MeV_xoshiro256pp_IS-x2per5mfp.egsinp; confirm it fails toconverge beyond ~25 mfp (known limitation of constant-factor schemes).