handle edge case lo == hi when the nominal value is 0

[github-actions]

LineCableModels.jl/src/uq/UQ.jl

Lines 42 to 43 in 2502c78

	# TODO: handle edge case lo == hi when the nominal value is 0
	# Issue URL: https://github.com/Electa-Git/LineCableModels.jl/issues/31

module UQ

# Export public API
export sample, trial, mc, hist

# Module-specific dependencies
using ..Commons: BASE_FLOAT
import ..Commons: domain, PhaseDomain, ModalDomain,
	LineParamsDomain
using ..ParametricBuilder:
	MaterialSpec, PartSpec, CableBuilderSpec, SystemBuilderSpec, AbstractPositionSpec,
	PositionSpec, PositionGroupSpec, build, iterate, _spec, determinize
using ..Engine:
	EMTFormulation, compute!, LineParameters
using ..DataModel: get_outer_radius
using Measurements: Measurement, measurement, value, uncertainty
using Random, Statistics, DataFrames
using Distributions:
	Distributions, ContinuousUnivariateDistribution, Normal, Uniform, cdf, sampler
using StatsBase: fit, Histogram, normalize, quantile, ecdf
using LinearAlgebra


# Draw once from a "range-like" spec
#   spec :: Number                             → return as-is
#   spec :: AbstractVector                     → random element (uniform over indices)
#   spec :: (lo::Number, hi::Number, n::Int)   → given [lo, hi], interpret as ±1σ around μ = (lo+hi)/2, σ = (hi-lo)/2.
#   anything iterable                          → pick a random element
@inline function _rand_in(spec, distribution::Symbol)
	if spec isa Number
		return spec

	elseif spec isa AbstractVector
		@inbounds return spec[rand(1:length(spec))]

	elseif spec isa Tuple && length(spec) == 3 &&
		   spec[1] isa Number && spec[2] isa Number && spec[3] isa Integer
		lo, hi = spec[1], spec[2]
		# - Given [lo, hi], interpret as ±1σ around μ = (lo+hi)/2, σ = (hi-lo)/2.
		lo_f = float(lo)
		hi_f = float(hi)
		# TODO: handle edge case lo == hi when the nominal value is 0
		# Issue URL: https://github.com/Electa-Git/LineCableModels.jl/issues/31
		@assert hi_f > lo_f "hi must be greater than lo"
		μ = (lo_f + hi_f) / 2
		σ = (hi_f - lo_f) / 2

		if distribution === :normal
			# - :normal => Normal(μ, σ).
			return rand(Distributions.Normal(μ, σ))
		elseif distribution === :uniform
			# - :uniform  => Uniform(μ ± √3 σ) so std matches σ.
			d = √3 * σ
			return rand(Distributions.Uniform(μ - d, μ + d))
		else
			throw(
				ArgumentError(
					"unsupported distribution: $(distribution). Use :uniform or :normal",
				),
			)
		end

	else
		if Base.iterable(spec)
			vals = collect(spec)
			@inbounds return vals[rand(1:length(vals))]
		end
		return spec
	end
end

# Collapse a (spec, pct) pair → (value::Number, pct::Union{Nothing,Number})
"""
	_collapse_pair(sp::Tuple, distribution::Symbol; domain=nothing, max_tries::Int=10_000)

Collapse a (spec, pct) pair into `(value, pct_value)` by drawing once from the
"range-like" `spec` and `pct` using `_rand_in`.