From 9c4961a1c079a01bff9771caabfdc57c47ea3467 Mon Sep 17 00:00:00 2001 From: Ryan Senne <50930199+rsenne@users.noreply.github.com> Date: Mon, 18 May 2026 08:10:32 -0400 Subject: [PATCH 1/4] Adds type promotion in all .*intensity methods --- src/univariate/hawkes/hawkes_process.jl | 7 +- .../inhomogeneous/intensity_functions.jl | 10 +- .../inhomogeneous/intensity_methods.jl | 91 ++++++------ src/univariate/poisson/poisson_process.jl | 7 +- test/hawkes.jl | 10 ++ test/inhomogeneous_poisson_process.jl | 131 ++++++++++++++++++ test/poisson_process.jl | 13 ++ 7 files changed, 214 insertions(+), 55 deletions(-) diff --git a/src/univariate/hawkes/hawkes_process.jl b/src/univariate/hawkes/hawkes_process.jl index a0e9032..ebfaacd 100644 --- a/src/univariate/hawkes/hawkes_process.jl +++ b/src/univariate/hawkes/hawkes_process.jl @@ -186,9 +186,12 @@ function ground_intensity(hp::HawkesProcess, h::History, t) return hp.μ + (hp.α * activation / exp(hp.ω * t)) end -function integrated_ground_intensity(hp::HawkesProcess, h::History, tmin, tmax) +function integrated_ground_intensity( + hp::HawkesProcess{T}, h::History, tmin, tmax +) where {T} + U = promote_type(T, typeof(tmin), typeof(tmax)) times = event_times(h, h.tmin, tmax) - integral = 0 + integral = zero(U) for ti in times # Integral of activation function. 'max(tmin - ti, 0)' corrects for events that occurred # inside or outside the interval [tmin, tmax]. diff --git a/src/univariate/poisson/inhomogeneous/intensity_functions.jl b/src/univariate/poisson/inhomogeneous/intensity_functions.jl index 069e302..dc0170e 100644 --- a/src/univariate/poisson/inhomogeneous/intensity_functions.jl +++ b/src/univariate/poisson/inhomogeneous/intensity_functions.jl @@ -51,8 +51,9 @@ struct PolynomialIntensity{R<:Real,L} <: ParametricIntensity end function (f::PolynomialIntensity)(t) - η = f.coefficients[1] - t_power = one(t) + U = promote_type(eltype(f.coefficients), typeof(t)) + η = convert(U, f.coefficients[1]) + t_power = one(U) for i in 2:length(f.coefficients) t_power *= t η += f.coefficients[i] * t_power @@ -217,11 +218,12 @@ struct PiecewiseConstantIntensity{R<:Real} end function (f::PiecewiseConstantIntensity)(t) + U = promote_type(eltype(f.rates), typeof(t)) idx = searchsortedlast(f.breakpoints, t) if idx == 0 || idx == length(f.breakpoints) - return zero(eltype(f.rates)) + return zero(U) end - return f.rates[idx] + return convert(U, f.rates[idx]) end function Base.show(io::IO, f::PiecewiseConstantIntensity) diff --git a/src/univariate/poisson/inhomogeneous/intensity_methods.jl b/src/univariate/poisson/inhomogeneous/intensity_methods.jl index db28af5..4f62b47 100644 --- a/src/univariate/poisson/inhomogeneous/intensity_methods.jl +++ b/src/univariate/poisson/inhomogeneous/intensity_methods.jl @@ -11,20 +11,19 @@ Analytical integral for polynomial intensity functions with identity link. For log link, we use numerical integration. =# function integrated_intensity( - f::PolynomialIntensity, lb::T, ub::T, config::IntegrationConfig -) where {T} + f::PolynomialIntensity{R}, lb::T, ub::T, config::IntegrationConfig +) where {R,T} + U = promote_type(R, T) if f.link === :identity - # Analytical integral for polynomial - result = zero(promote_type(eltype(f.coefficients), typeof(lb), typeof(ub))) + result = zero(U) for (i, coef) in enumerate(f.coefficients) power = i result += coef * (ub^power - lb^power) / power end return result else - # Numerical integration for log link using config integrand(t, p) = f(t) - prob = IntegralProblem(integrand, (lb, ub)) # 1D domain (lb, ub) + prob = IntegralProblem(integrand, (lb, ub)) integral = solve( prob, config.solver; @@ -32,7 +31,7 @@ function integrated_intensity( reltol=config.reltol, maxiters=config.maxiters, ) - return integral.u + return convert(U, integral.u) end end @@ -44,12 +43,13 @@ For polynomials, we sample densely and add a margin. function intensity_bound( f::PolynomialIntensity{R}, t::T; - lookahead::T=one(T), # default if you don't have history - n_samples::Int=10 * length(f.coefficients), # scale with degree + lookahead::T=one(T), + n_samples::Int=10 * length(f.coefficients), ) where {R,T} + U = promote_type(R, T) ts = range(t, t + lookahead; length=n_samples) max_val = maximum(f(ti) for ti in ts) - return (max_val * T(1.1), lookahead) + return (U(max_val * 1.1), U(lookahead)) end function intensity_bound( @@ -66,15 +66,15 @@ end ## ExponentialIntensity optimizations -# Analytical integral for exponential intensity: ∫ a*exp(b*t) dt = (a/b)*(exp(b*b) - exp(b*a)) +# Analytical integral for exponential intensity: ∫ a*exp(b*t) dt = (a/b)*(exp(b*ub) - exp(b*lb)) function integrated_intensity( - f::ExponentialIntensity, lb::T, ub::T, config::IntegrationConfig -) where {T} + f::ExponentialIntensity{R}, lb::T, ub::T, config::IntegrationConfig +) where {R,T} + U = promote_type(R, T) if abs(f.b) < 1e-10 - # b ≈ 0, treat as constant - return f.a * (ub - lb) + return U(f.a * (ub - lb)) end - return (f.a / f.b) * (exp(f.b * ub) - exp(f.b * lb)) + return U((f.a / f.b) * (exp(f.b * ub) - exp(f.b * lb))) end #= @@ -85,17 +85,19 @@ If b < 0 (decreasing), max is at left endpoint. If b ≈ 0 (constant), use constant bound. =# function intensity_bound(f::ExponentialIntensity{R}, t::T; lookahead::T=one(T)) where {R,T} + U = promote_type(R, T) + margin = U(1.05) if f.b > 1e-10 # Increasing: max at t + lookahead - B = f(t + lookahead) * T(1.05) # small margin + B = f(t + lookahead) * margin elseif f.b < -1e-10 # Decreasing: max at t - B = f(t) * T(1.05) + B = f(t) * margin else # Approximately constant - B = f.a * T(1.05) + B = U(f.a) * margin end - return (B, lookahead) + return (U(B), U(lookahead)) end function intensity_bound( @@ -112,16 +114,16 @@ end Analytical integral for sinusoidal intensity: ∫ (a + b*sin(ω*t + φ)) dt =# function integrated_intensity( - f::SinusoidalIntensity, lb::T, ub::T, config::IntegrationConfig -) where {T} + f::SinusoidalIntensity{R}, lb::T, ub::T, config::IntegrationConfig +) where {R,T} + U = promote_type(R, T) linear_part = f.a * (ub - lb) if abs(f.ω) < 1e-10 - # ω ≈ 0, sin term is approximately constant sin_part = f.b * sin(f.φ) * (ub - lb) else sin_part = -(f.b / f.ω) * (cos(f.ω * ub + f.φ) - cos(f.ω * lb + f.φ)) end - return linear_part + sin_part + return U(linear_part + sin_part) end #= @@ -130,48 +132,43 @@ Upper bound for sinusoidal intensity. Maximum is a + |b| (when sin = 1 if b > 0, or sin = -1 if b < 0). =# function intensity_bound(f::SinusoidalIntensity{R}, t::T, h::History) where {R,T} - B = f.a + abs(f.b) - L = typemax(T) # Bound holds for all time + U = promote_type(R, T) + B = U(f.a + abs(f.b)) + L = typemax(U) # Bound holds for all time return (B, L) end ## PiecewiseConstantIntensity optimizations function integrated_intensity( - f::PiecewiseConstantIntensity, lb::T, ub::T, config::IntegrationConfig -) where {T} - # Find the intervals that overlap with [lb, ub] + f::PiecewiseConstantIntensity{R}, lb::T, ub::T, config::IntegrationConfig +) where {R,T} + U = promote_type(R, T) start_idx = searchsortedlast(f.breakpoints, lb) end_idx = searchsortedlast(f.breakpoints, ub) - # Clamp to valid range start_idx = max(1, min(start_idx, length(f.rates))) end_idx = max(1, min(end_idx, length(f.rates))) if start_idx == end_idx - # Entire interval [lb, ub] is within a single constant region if start_idx > 0 && start_idx <= length(f.rates) - return f.rates[start_idx] * (ub - lb) + return U(f.rates[start_idx] * (ub - lb)) else - return zero(eltype(f.rates)) + return zero(U) end end - # Sum over multiple intervals - integral = zero(eltype(f.rates)) + integral = zero(U) - # First partial interval if start_idx > 0 && start_idx <= length(f.rates) integral += f.rates[start_idx] * (f.breakpoints[start_idx + 1] - lb) end - # Complete intervals in between for i in (start_idx + 1):(end_idx - 1) if i > 0 && i <= length(f.rates) integral += f.rates[i] * (f.breakpoints[i + 1] - f.breakpoints[i]) end end - # Last partial interval if end_idx > 0 && end_idx <= length(f.rates) integral += f.rates[end_idx] * (ub - f.breakpoints[end_idx]) end @@ -180,17 +177,17 @@ function integrated_intensity( end function intensity_bound(f::PiecewiseConstantIntensity{R}, t::T, h::History) where {R,T} + U = promote_type(R, T) idx = searchsortedlast(f.breakpoints, t) if idx == 0 || idx >= length(f.breakpoints) - # Outside the domain - return (zero(R), typemax(T)) + return (zero(U), typemax(U)) end - # Return current rate (exact) and time to next breakpoint - # This avoids rejections in Ogata's algorithm within constant pieces - current_rate = f.rates[idx] - time_to_next_breakpoint = f.breakpoints[idx + 1] - t + # Return current rate (exact) and time to next breakpoint. + # This avoids rejections in Ogata's algorithm within constant pieces. + current_rate = U(f.rates[idx]) + time_to_next_breakpoint = U(f.breakpoints[idx + 1] - t) return (current_rate, time_to_next_breakpoint) end @@ -215,9 +212,10 @@ end function intensity_bound( f::LinearCovariateIntensity{R}, t::T; lookahead::T=one(T), n_samples::Int=100 ) where {R,T} + U = promote_type(R, T) ts = range(t, t + lookahead; length=n_samples) max_intensity = maximum(f(ti) for ti in ts) - return (max_intensity * T(1.1), lookahead) + return (U(max_intensity * 1.1), U(lookahead)) end function intensity_bound( @@ -252,7 +250,8 @@ end function intensity_bound(f::F, t::T; lookahead::T=one(T), n_samples::Int=100) where {F,T} ts = range(t, t + lookahead; length=n_samples) max_intensity = maximum(f(ti) for ti in ts) - return (max_intensity * T(1.1), lookahead) + U = promote_type(typeof(max_intensity), T) + return (U(max_intensity * 1.1), U(lookahead)) end function intensity_bound( diff --git a/src/univariate/poisson/poisson_process.jl b/src/univariate/poisson/poisson_process.jl index 4b13ad9..f6896cf 100644 --- a/src/univariate/poisson/poisson_process.jl +++ b/src/univariate/poisson/poisson_process.jl @@ -66,9 +66,10 @@ mark_distribution(pp::PoissonProcess, t) = mark_distribution(pp) # For simulate_ intensity(pp::PoissonProcess, m, t, h) = intensity(pp, m) log_intensity(pp::PoissonProcess, m, t, h) = log_intensity(pp, m) -function ground_intensity_bound(pp::PoissonProcess, t::T, h) where {T<:Real} - B = ground_intensity(pp) - L = typemax(T) +function ground_intensity_bound(pp::PoissonProcess{R}, t::T, h) where {R,T<:Real} + U = promote_type(R, T) + B = U(ground_intensity(pp)) + L = typemax(U) return (B, L) end diff --git a/test/hawkes.jl b/test/hawkes.jl index abc6fda..009d058 100644 --- a/test/hawkes.jl +++ b/test/hawkes.jl @@ -55,3 +55,13 @@ params_est = (model_est.μ, model_est.α, model_est.ω) @test logdensityof(hp, h) ≈ sum(log.(hp.μ .+ (hp.α .* [0, exp(-hp.ω), exp(-hp.ω * 2) + exp(-hp.ω * 3)]))) - integral + +# Type promotion in integrated_ground_intensity (accumulator was Int before) +hp32 = HawkesProcess(0.5f0, 0.1f0, 1.0f0) +h32 = History(Float32[1, 2], 0.0f0, 5.0f0) +@test typeof(integrated_ground_intensity(hp32, h32, 0.0f0, 5.0f0)) === Float32 +# Empty history exercises the no-loop path where Int(0) used to leak +h32_empty = History(Float32[], 0.0f0, 5.0f0) +@test typeof(integrated_ground_intensity(hp32, h32_empty, 0.0f0, 5.0f0)) === Float32 +# Mixed: Float32 params + Float64 endpoints → Float64 +@test typeof(integrated_ground_intensity(hp32, h32, 0.0, 5.0)) === Float64 diff --git a/test/inhomogeneous_poisson_process.jl b/test/inhomogeneous_poisson_process.jl index 9d18257..50c92e1 100644 --- a/test/inhomogeneous_poisson_process.jl +++ b/test/inhomogeneous_poisson_process.jl @@ -1148,3 +1148,134 @@ end @test pp_est.mark_dist.value === nothing end end + +@testset "Type promotion in intensity methods" begin + config = IntegrationConfig() + h32 = History(Float32[0.5, 1.5, 2.5], 0.0f0, 3.0f0) + + @testset "Callable: PolynomialIntensity" begin + f64 = PolynomialIntensity([1.0, 0.5]) + f32 = PolynomialIntensity(Float32[1.0, 0.5]) + f_int_const = PolynomialIntensity([2]) # length-1, Int coeffs + f_f64_const = PolynomialIntensity([2.0]) # length-1, Float64 coeffs + + @test typeof(f64(2.0)) === Float64 + @test typeof(f32(2.0f0)) === Float32 # Float32 preserved + @test typeof(f64(2.0f0)) === Float64 # mixed → promote + @test typeof(f32(2.0)) === Float64 # mixed → promote + @test typeof(f_int_const(1.5)) === Float64 # length-1 used to drop typeof(t) + @test typeof(f_f64_const(1.0f0)) === Float64 + end + + @testset "Callable: PiecewiseConstantIntensity" begin + f32 = PiecewiseConstantIntensity(Float32[0, 1, 2], Float32[1, 2]) + + @test typeof(f32(0.5f0)) === Float32 # Float32 preserved + @test typeof(f32(0.5)) === Float64 # mixed → promote (was Float32 bug) + @test typeof(f32(-1.0)) === Float64 # out-of-domain promoted + @test typeof(f32(-1.0f0)) === Float32 + @test f32(-1.0f0) === 0.0f0 # value preserved + end + + @testset "integrated_intensity: uniform return type" begin + f_p32 = PolynomialIntensity(Float32[1.0, 0.5]) + f_e32 = ExponentialIntensity(1.0f0, 0.1f0) + f_s32 = SinusoidalIntensity(3.0f0, 1.0f0, Float32(2π), 0.0f0) + f_pw32 = PiecewiseConstantIntensity(Float32[0, 1, 2], Float32[1, 2]) + + @test typeof(PointProcesses.integrated_intensity(f_p32, 0.0f0, 2.0f0, config)) === Float32 + @test typeof(PointProcesses.integrated_intensity(f_e32, 0.0f0, 2.0f0, config)) === Float32 + @test typeof(PointProcesses.integrated_intensity(f_s32, 0.0f0, 2.0f0, config)) === Float32 + @test typeof(PointProcesses.integrated_intensity(f_pw32, 0.0f0, 2.0f0, config)) === Float32 + + # Mixed types: Float32 params + Float64 endpoints → Float64 + @test typeof(PointProcesses.integrated_intensity(f_p32, 0.0, 2.0, config)) === Float64 + @test typeof(PointProcesses.integrated_intensity(f_e32, 0.0, 2.0, config)) === Float64 + end + + @testset "intensity_bound: tuple components share type" begin + f_p32 = PolynomialIntensity(Float32[1.0, 0.5]) + f_e32 = ExponentialIntensity(1.0f0, 0.1f0) + f_s32 = SinusoidalIntensity(3.0f0, 1.0f0, Float32(2π), 0.0f0) + f_pw32 = PiecewiseConstantIntensity(Float32[0, 1, 2], Float32[1, 2]) + + # Same-type Float32 round-trip across all four parametric intensity types + for tup in ( + PointProcesses.intensity_bound(f_p32, 0.0f0; lookahead=1.0f0), + PointProcesses.intensity_bound(f_e32, 0.0f0; lookahead=1.0f0), + PointProcesses.intensity_bound(f_s32, 0.0f0, h32), + PointProcesses.intensity_bound(f_pw32, 0.5f0, h32), + ) + @test typeof(tup[1]) === typeof(tup[2]) === Float32 + end + + # Mixed-type: Float64 params, Float32 t → both Float64 + f_p64 = PolynomialIntensity([1.0, 0.5]) + tup = PointProcesses.intensity_bound(f_p64, 0.0f0; lookahead=1.0f0) + @test typeof(tup[1]) === typeof(tup[2]) === Float64 + + # Generic fallback for arbitrary callables + custom = t -> 1.0f0 + sin(t) + tup = PointProcesses.intensity_bound(custom, 0.0f0; lookahead=1.0f0) + @test typeof(tup[1]) === typeof(tup[2]) + end + + @testset "Process-level Float32 round-trip" begin + f32 = PolynomialIntensity(Float32[1.0, 0.5]) + pp32 = InhomogeneousPoissonProcess(f32) + + @test typeof(ground_intensity(pp32, 0.0f0, h32)) === Float32 + @test typeof(integrated_ground_intensity(pp32, h32, 0.0f0, 2.0f0)) === Float32 + tup = ground_intensity_bound(pp32, 0.0f0, h32) + @test typeof(tup[1]) === typeof(tup[2]) === Float32 + end + + @testset "ForwardDiff: intensity callables" begin + # λ(t) = p₀ + p₁·t → ∂λ/∂p₀=1, ∂λ/∂p₁=t + g = ForwardDiff.gradient(p -> PolynomialIntensity(p)(2.0), [1.0, 0.5]) + @test g ≈ [1.0, 2.0] + + # λ(t) = a·exp(b·t) → ∂λ/∂a=exp(b·t), ∂λ/∂b=a·t·exp(b·t) + g = ForwardDiff.gradient(p -> ExponentialIntensity(p[1], p[2])(2.0), [1.0, 0.1]) + @test g[1] ≈ exp(0.2) + @test g[2] ≈ 1.0 * 2.0 * exp(0.2) + + # λ(t) = a + b·sin(ω·t + φ) → ∂λ/∂a = 1 + g = ForwardDiff.gradient( + p -> SinusoidalIntensity(p[1], p[2], p[3], p[4])(0.25), [3.0, 1.0, 2π, 0.0] + ) + @test all(isfinite, g) + @test g[1] ≈ 1.0 + end + + @testset "ForwardDiff: end-to-end through process intensity methods" begin + # Gradient flows through both ground_intensity (per-event evaluation) + # and integrated_ground_intensity (analytical and numerical paths). + # Built as Λ - Σ log λ(tᵢ) — same shape as the unmarked log-likelihood + # used by the fitter, but without the mark-distribution call so the + # test isolates the intensity plumbing. + h = History([0.5, 1.2, 2.7], 0.0, 4.0) + nll(pp, h) = + integrated_ground_intensity(pp, h, h.tmin, h.tmax) - + sum(log(ground_intensity(pp, t, h)) for t in h.times) + + g_poly = ForwardDiff.gradient( + p -> nll(InhomogeneousPoissonProcess(PolynomialIntensity(p)), h), [1.0, 0.3] + ) + @test all(isfinite, g_poly) + + g_exp = ForwardDiff.gradient( + p -> nll(InhomogeneousPoissonProcess(ExponentialIntensity(p[1], p[2])), h), + [1.0, 0.1], + ) + @test all(isfinite, g_exp) + + g_sin = ForwardDiff.gradient( + p -> nll( + InhomogeneousPoissonProcess(SinusoidalIntensity(p[1], p[2], 2π, 0.0)), h + ), + [3.0, 1.0], + ) + @test all(isfinite, g_sin) + end +end diff --git a/test/poisson_process.jl b/test/poisson_process.jl index 23ea80b..53fcdb7 100644 --- a/test/poisson_process.jl +++ b/test/poisson_process.jl @@ -82,3 +82,16 @@ end l_est = logdensityof(pp_est1, h1) @test l_est > l end + +@testset "Type promotion in ground_intensity_bound" begin + pp32 = PoissonProcess(1.0f0) + h32 = History(Float32[], 0.0f0, 1.0f0) + tup = ground_intensity_bound(pp32, 0.0f0, h32) + @test typeof(tup[1]) === typeof(tup[2]) === Float32 # Float32 preserved + @test tup[2] === typemax(Float32) + + # Mixed: Float32 λ + Float64 t → both Float64 + h64 = History(Float64[], 0.0, 1.0) + tup = ground_intensity_bound(pp32, 0.0, h64) + @test typeof(tup[1]) === typeof(tup[2]) === Float64 +end From 50c9256f90813b46190705623c165716a82a6688 Mon Sep 17 00:00:00 2001 From: Ryan Senne <50930199+rsenne@users.noreply.github.com> Date: Mon, 18 May 2026 08:29:32 -0400 Subject: [PATCH 2/4] Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/univariate/hawkes/hawkes_process.jl | 4 +--- test/inhomogeneous_poisson_process.jl | 18 ++++++++++++------ 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/src/univariate/hawkes/hawkes_process.jl b/src/univariate/hawkes/hawkes_process.jl index ebfaacd..57ffad6 100644 --- a/src/univariate/hawkes/hawkes_process.jl +++ b/src/univariate/hawkes/hawkes_process.jl @@ -186,9 +186,7 @@ function ground_intensity(hp::HawkesProcess, h::History, t) return hp.μ + (hp.α * activation / exp(hp.ω * t)) end -function integrated_ground_intensity( - hp::HawkesProcess{T}, h::History, tmin, tmax -) where {T} +function integrated_ground_intensity(hp::HawkesProcess{T}, h::History, tmin, tmax) where {T} U = promote_type(T, typeof(tmin), typeof(tmax)) times = event_times(h, h.tmin, tmax) integral = zero(U) diff --git a/test/inhomogeneous_poisson_process.jl b/test/inhomogeneous_poisson_process.jl index 92c77b3..a374039 100644 --- a/test/inhomogeneous_poisson_process.jl +++ b/test/inhomogeneous_poisson_process.jl @@ -1478,14 +1478,20 @@ end f_s32 = SinusoidalIntensity(3.0f0, 1.0f0, Float32(2π), 0.0f0) f_pw32 = PiecewiseConstantIntensity(Float32[0, 1, 2], Float32[1, 2]) - @test typeof(PointProcesses.integrated_intensity(f_p32, 0.0f0, 2.0f0, config)) === Float32 - @test typeof(PointProcesses.integrated_intensity(f_e32, 0.0f0, 2.0f0, config)) === Float32 - @test typeof(PointProcesses.integrated_intensity(f_s32, 0.0f0, 2.0f0, config)) === Float32 - @test typeof(PointProcesses.integrated_intensity(f_pw32, 0.0f0, 2.0f0, config)) === Float32 + @test typeof(PointProcesses.integrated_intensity(f_p32, 0.0f0, 2.0f0, config)) === + Float32 + @test typeof(PointProcesses.integrated_intensity(f_e32, 0.0f0, 2.0f0, config)) === + Float32 + @test typeof(PointProcesses.integrated_intensity(f_s32, 0.0f0, 2.0f0, config)) === + Float32 + @test typeof(PointProcesses.integrated_intensity(f_pw32, 0.0f0, 2.0f0, config)) === + Float32 # Mixed types: Float32 params + Float64 endpoints → Float64 - @test typeof(PointProcesses.integrated_intensity(f_p32, 0.0, 2.0, config)) === Float64 - @test typeof(PointProcesses.integrated_intensity(f_e32, 0.0, 2.0, config)) === Float64 + @test typeof(PointProcesses.integrated_intensity(f_p32, 0.0, 2.0, config)) === + Float64 + @test typeof(PointProcesses.integrated_intensity(f_e32, 0.0, 2.0, config)) === + Float64 end @testset "intensity_bound: tuple components share type" begin From 65c1d4677056775d9ba97195df213251b5000ebb Mon Sep 17 00:00:00 2001 From: Ryan Senne <50930199+rsenne@users.noreply.github.com> Date: Mon, 18 May 2026 08:48:06 -0400 Subject: [PATCH 3/4] Add degenerate test case --- test/inhomogeneous_poisson_process.jl | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/test/inhomogeneous_poisson_process.jl b/test/inhomogeneous_poisson_process.jl index 92c77b3..febf013 100644 --- a/test/inhomogeneous_poisson_process.jl +++ b/test/inhomogeneous_poisson_process.jl @@ -945,6 +945,29 @@ end @test integral_cross ≈ 12.0 rtol = 1e-6 end + @testset "PiecewiseConstantIntensity - empty rates (degenerate)" begin + # Single-breakpoint intensity has zero rate-intervals. This is the only + # input that reaches the `return zero(U)` branch of integrated_intensity: + # with any non-empty rates the post-clamp indices land in [1, length(rates)] + # and the inner `start_idx <= length(f.rates)` is always true. + intensity_empty = PiecewiseConstantIntensity([0.0], Float64[]) + + # Same-type call returns zero of the promoted type + @test PointProcesses.integrated_intensity(intensity_empty, -1.0, 1.0, config) === 0.0 + @test PointProcesses.integrated_intensity(intensity_empty, 0.5, 0.5, config) === 0.0 + + # Mixed-type call: Float64 rates + Float32 endpoints → Float64 + @test PointProcesses.integrated_intensity( + intensity_empty, -1.0f0, 1.0f0, config + ) === 0.0 + + # Float32 rates path + intensity_empty32 = PiecewiseConstantIntensity(Float32[0], Float32[]) + @test PointProcesses.integrated_intensity( + intensity_empty32, -1.0f0, 1.0f0, config + ) === 0.0f0 + end + @testset "LinearCovariateIntensity numerical integration" begin # λ(t) = exp(1.0 + 0.5*t + 0.2*sin(t)) cov1 = t -> t From 7a03df22107a03f2e4db8ea5707f882231a9fa76 Mon Sep 17 00:00:00 2001 From: Ryan Senne <50930199+rsenne@users.noreply.github.com> Date: Mon, 18 May 2026 08:49:45 -0400 Subject: [PATCH 4/4] Update test/inhomogeneous_poisson_process.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- test/inhomogeneous_poisson_process.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/inhomogeneous_poisson_process.jl b/test/inhomogeneous_poisson_process.jl index 668a100..b6cefb4 100644 --- a/test/inhomogeneous_poisson_process.jl +++ b/test/inhomogeneous_poisson_process.jl @@ -953,7 +953,8 @@ end intensity_empty = PiecewiseConstantIntensity([0.0], Float64[]) # Same-type call returns zero of the promoted type - @test PointProcesses.integrated_intensity(intensity_empty, -1.0, 1.0, config) === 0.0 + @test PointProcesses.integrated_intensity(intensity_empty, -1.0, 1.0, config) === + 0.0 @test PointProcesses.integrated_intensity(intensity_empty, 0.5, 0.5, config) === 0.0 # Mixed-type call: Float64 rates + Float32 endpoints → Float64