Skip to content

feat(Sobol): add non-finite (NaN/Inf) result filtering across nboot, … - #261

Open
bmchun722 wants to merge 4 commits into
SciML:masterfrom
bmchun722:sciml-sobol-filter
Open

feat(Sobol): add non-finite (NaN/Inf) result filtering across nboot, …#261
bmchun722 wants to merge 4 commits into
SciML:masterfrom
bmchun722:sciml-sobol-filter

Conversation

@bmchun722

@bmchun722 bmchun722 commented Aug 11, 2026

Copy link
Copy Markdown

…S2, batch, and multioutput

Checklist

  • [o] Appropriate tests were added
  • [o] Any code changes were done in a way that does not break public API
  • [o] All documentation related to code changes were updated
  • [o] The new code follows the
    contributor guidelines, in particular the SciML Style Guide and
    COLPRAC.
  • [o] Any new documentation only uses public API

Additional context

Add any other context about the problem here.

PR Description

Feature Summary

When conducting Global Sensitivity Analysis (GSA) on complex physical or ODE models, specific parameter combinations may lead to non-convergence, stiffness failures, or numerical overflow/underflow, resulting in non-finite outputs (NaN / Inf). Without filtering, even a single NaN sample propagates through Monte Carlo variance decompositions and pollutes the entire Sobol output ($S_1, S_2, S_T$, $V(Y)$, and confidence intervals).

This PR introduces a robust, generalized sample filtering mechanism (_compact_all_y_nonfinite!) to GlobalSensitivity.Sobol (src/sobol_sensitivity.jl). It automatically identifies and filters out non-finite sample points while maintaining design matrix alignment across:

  • Bootstrap confidence interval calculations (nboot ≥ 1)
  • First-order ($S_1$) and Second-order ($S_2$) interaction indices (order = [0, 1, 2])
  • Batch evaluation mode (batch = true / false)
  • Multi-output and time-series models (vector / matrix evaluations)

In addition, SobolResult is updated to report:

  • VY: total output variance $V(Y)$ (pooled across valid samples)
  • n: count of valid, non-filtered design sample points ($nk = \text{count}(keep)$)

Memory Layout & Masking Mechanism

1. Unified Linear Memory Layout

In GlobalSensitivity.Sobol, model evaluations all_y are flattened into a contiguous 1D vector (for scalar outputs) or a 2D matrix of shape (num_outputs, total_evaluations) (for multioutput models).

Given:

  • $n$: number of design points per design matrix
  • $d$: number of input parameters
  • $nboot$: number of bootstrap replicates ($b = 0 \dots nboot - 1$)
  • $nblocks$: number of design sub-blocks per replicate
    • 1st-order (second_order = false): $nblocks = 2 + d$ ($A, B, AB_1, \dots, AB_d$)
    • 2nd-order (second_order = true): $nblocks = 2 + 2d$ ($A, B, A_{b1 \dots d}, B_{a1 \dots d}$)

The total number of evaluations across all bootstrap replicates and blocks is $N_{\text{total}} = nboot \times nblocks \times n$.
For a specific sample point $k \in {1 \dots n}$, bootstrap replicate $b$, and block $bl$, the linear column index in all_y is:

$$\text{idx}(b, bl, k) = b \times (nblocks \times n) + bl \times n + k$$

2. Unified Masking Logic with keep

To ensure mathematical consistency across all design blocks, bootstrap replicates, and output components, the masking logic is unified across all dimensions:

  1. A boolean mask keep = Base.fill(true, n) of size $n$ is initialized.
  2. Sample point $k$ is marked valid (keep[k] = true) if and only if the evaluation at $\text{idx}(b, bl, k)$ is finite.
# 1D Vector (Scalar output)
@inbounds for b in 0:(nboot - 1)
    b_offset = b * (nblocks * n)
    for bl in 0:(nblocks - 1)
        block_offset = b_offset + bl * n
        for k in 1:n
            if !isfinite(all_y[block_offset + k])
                keep[k] = false
            end
        end
    end
end

# 2D Matrix (Multi-output)
num_outputs = size(all_y, 1)
@inbounds for b in 0:(nboot - 1)
    b_offset = b * (nblocks * n)
    for bl in 0:(nblocks - 1)
        block_offset = b_offset + bl * n
        for k in 1:n
            col = block_offset + k
            for r in 1:num_outputs
                if !isfinite(all_y[r, col])
                    keep[k] = false
                    break
                end
            end
        end
    end
end

3. Block-wise Compact Copying

If $nk = \text{count}(keep) < n$, a compact buffer filtered of size $nboot \times nblocks \times nk$ (or (num_outputs, nboot * nblocks * nk) for 2D) is allocated. The valid $nk$ entries for each block and replicate are copied in-place while preserving the block sequence:

filtered = similar(all_y, eltype(all_y), nboot * nblocks * nk)
dest = 1
@inbounds for b in 0:(nboot - 1)
    b_offset = b * (nblocks * n)
    for bl in 0:(nblocks - 1)
        block_offset = b_offset + bl * n
        for k in 1:n
            if keep[k]
                filtered[dest] = all_y[block_offset + k]
                dest += 1
            end
        end
    end
end

The resulting filtered array retains the exact layout expected by gsa_sobol_all_y_analysis, with each block now having length $nk$. The valid sample count $nk$ is returned as SobolResult.n.


Verification & Testing

  1. Unit Tests Added (test/sobol_method.jl):
    • Added test functions (ishi_nan, ishi_nan_batch, ishi_linear_nan) that explicitly return NaN for subset bounds to verify filtering under 1st order, 2nd order, bootstrap runs (nboot = 5), batch mode, and multioutput setups.
  2. Full Package Test Suite (Pkg.test()):
    • Executed Pkg.test(): All 12 GSA test suites passed 100% (Sobol 42/42, Morris 23/23, Shapley 13/13, etc.).
  3. Runic Formatting:
    • Verified code formatting compliance via Runic.main(["--check", "src/sobol_sensitivity.jl"]).

This description and the code was generated with the assist of LLM

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant