Systematize intensity types

src/univariate/hawkes/hawkes_process.jl

@@ -186,9 +186,10 @@ 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].

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)

src/univariate/poisson/inhomogeneous/intensity_methods.jl

@@ -11,28 +11,27 @@ 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;
             abstol=config.abstol,
             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(

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
 

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