Is there an easy way to speed this up?
sims = 10_000
T = 1
σ = .1
Δ = .000001
σrΔ = σ√Δ
steps = (Int64 ∘ round)(T/Δ)
L = .9
H = 1.1
x = ones( sims )
alive = trues( sims )
for t ∈ 1:steps
x .+= randn(sims)σrΔ
@. alive &= L < x < H
end
mean(alive)
around 2x faster
julia> function main()
sims = 10_000
T = 1
σ = .1
Δ = .000001
σrΔ = σ√Δ
steps = (Int64 ∘ round)(T/Δ)
L = .9
H = 1.1
x = ones( sims )
alive = trues( sims )
cache = similar(x)
for t ∈ 1:steps
x .+= Random.randn!(cache) .* σrΔ
@. alive &= L < x < H
end
mean(alive)end
main (generic function with 1 method)
julia> @time main()
18.912667 seconds (9 allocations: 157.811 KiB)
0.3624
you could also try to switch random generator to check.
(post deleted by author)
How fast you can go depends on what you want to record: if you only care about mean(alive), you can get very fast (if I didn’t make a mistake…).
sim_loop simply uses a loop instead of a vector of normals. You can then strip this further by using some problem-specific insights if you only care about mean(alive): 1) you can stop path-wise as soon as alive is false, 2) you don’t need to record alive as a vector, 3) instead of multiplying by a Float each time, you can also divide L and H by a Float once, which I would expect to be faster (but isn’t?)
using Random, BenchmarkTools
function sim_opt(sims::Int)
T = 1.0
σ = 0.1
Δ = 0.000001
σrΔ = σ * √Δ
steps = round(Int64, T / Δ)
L = 0.9
H = 1.1
x = ones(sims)
alive = trues(sims)
rnd = Vector{Float64}(undef, sims)
for t ∈ 1:steps
randn!(rnd)
@. x += rnd * σrΔ
@. alive &= (L < x < H)
end
mean(alive)
end
function sim_loop(sims::Int)
T = 1.0
σ = 0.1
Δ = 0.000001
σrΔ = σ * √Δ
steps = round(Int64, T / Δ)
L = 0.9
H = 1.1
x = ones(sims)
alive = trues(sims)
Threads.@threads for i in 1:sims
for _ in 1:steps
x[i] += randn() * σrΔ
alive[i] &= (L < x[i] < H)
end
end
mean(alive)
end
function sim_loop_only_alive(sims::Int)
T = 1.0
σ = 0.1
Δ = 0.000001
σrΔ = σ * √Δ
steps = round(Int64, T / Δ)
L = 0.9
H = 1.1
x = ones(sims)
alive = trues(sims)
Threads.@threads for i in 1:sims
for _ in 1:steps
x[i] += randn() * σrΔ
if !(L < x[i] < H)
alive[i] = false
break
end
end
end
mean(alive)
end
function sim_loop_only_alive_no_vecs(sims::Int)
T = 1.0
σ = 0.1
Δ = 0.000001
σrΔ = σ * √Δ
steps = round(Int64, T / Δ)
L = 0.9
H = 1.1
x₀ = 1.
n_alive = sims
Threads.@threads for i in 1:sims
x = x₀
for _ in 1:steps
x += randn() * σrΔ
if !(L < x < H)
n_alive -= 1
break
end
end
end
n_alive / sims
end
function sim_loop_only_alive_no_vecs_no_mul(sims::Int)
T = 1.0
σ = 0.1
Δ = 0.000001
σrΔ = σ * √Δ
steps = round(Int64, T / Δ)
L_sc = 0.9 / σrΔ
H_sc = 1.1 / σrΔ
x₀ = 1. / σrΔ
n_alive = sims
Threads.@threads for i in 1:sims
x = x₀
for _ in 1:steps
x += randn()
if !(L_sc < x < H_sc)
n_alive -= 1
break
end
end
end
n_alive / sims
end
function tests(N)
@time sim_opt(N);
@time sim_loop(N);
@time sim_loop_only_alive(N);
@time sim_loop_only_alive_no_vecs(N);
@time sim_loop_only_alive_no_vecs_no_mul(N);
end
tests(1000)
gives
3.063867 seconds (1.00 M allocations: 4.053 GiB, 6.50% gc time)
1.090280 seconds (16.19 k allocations: 1.043 MiB, 6.25% compilation time)
0.722454 seconds (14.61 k allocations: 995.672 KiB, 10.36% compilation time)
0.699596 seconds (13.65 k allocations: 907.141 KiB, 10.81% compilation time)
0.731058 seconds (13.49 k allocations: 896.547 KiB, 10.69% compilation time)
0.38
Note that if you only care about alive, a loop-based version will get relatively faster when you add more and more sims: a vectorized approach requires all calculations to be done for all paths, also for the ones that are no longer alive.
much appreciated.
could @turbo be used? or does the sequential calculation make it ineffective
When dealing with unreasonably slow code, it’s best to consider the algorithm. The example takes 24 seconds on my well provisioned Mac.
For a Brownian motion with constant drift=0, σ, and absorbing barriers at L and H starting from x₀=1 (the midpoint), there are closed-form survival probabilities via eigenfunction expansion. That would replace the entire simulation. Your Δ=1e-6 is very small, likely to avoid missing barrier crossings between discrete steps. Brownian bridge correction lets you use much larger Δ (e.g., 1e-3) while still accounting for intra-step barrier hits, reducing steps from 1M to 1K:
function barrier_sim_bb(; sims=10_000, T=1, σ=0.1, Δ=1e-3, L=0.9, H=1.1)
σrΔ = σ * √Δ
steps = round(Int, T / Δ)
x = ones(sims)
alive = trues(sims)
for t ∈ 1:steps
@inbounds for i ∈ 1:sims
alive[i] || continue
x_old = x[i]
x[i] += randn() * σrΔ
x_new = x[i]
# direct exit
if !(L < x_new < H)
alive[i] = false
continue
end
# Brownian bridge: probability of crossing barrier between steps
p_low = exp(-2(x_old - L) * (x_new - L) / (σ^2 * Δ))
p_high = exp(-2(H - x_old) * (H - x_new) / (σ^2 * Δ))
if rand() < p_low + p_high
alive[i] = false
end
end
end
mean(alive)
end
@time barrier_sim_bb(sims=10_000, T=1, σ=0.1, Δ=1e-3, L=0.9, H=1.1)
quite a bit faster
julia> @time barrier_sim_bb(sims=10_000, T=1, σ=0.1, Δ=1e-3, L=0.9, H=1.1)
0.093072 seconds (5 allocations: 81.422 KiB)
0.3762
I was thinking the same but figured that OP might have more complicated dynamics in mind in the long run. If not, directly simulating hitting times of fixed barriers is also an approach: first simulate the first hitting time of either barrier; then simulate the state x at the first of {hitting time, T} using a conditional distribution. If you also need the terminal state and hitting occurred before time T, simulate x_T using one normal random variate.
thank you.
I’ll combine this with JADekker’s approach.
Combining, I get
function barrier_sim_bb(sims=10_000; T=1, σ=0.1, Δ=1e-3, L=0.9, H=1.1)
σrΔ = σ * √Δ
steps = round(Int, T / Δ)
x = ones(sims)
alive = trues(sims)
for t ∈ 1:steps
@inbounds for i ∈ 1:sims
alive[i] || continue
x_old = x[i]
x[i] += randn() * σrΔ
x_new = x[i]
# direct exit
if !(L < x_new < H)
alive[i] = false
continue
end
# Brownian bridge: probability of crossing barrier between steps
p_low = exp(-2(x_old - L) * (x_new - L) / (σ^2 * Δ))
p_high = exp(-2(H - x_old) * (H - x_new) / (σ^2 * Δ))
if rand() < p_low + p_high
alive[i] = false
end
end
end
mean(alive)
end
function barrier_sim_bb_quick(sims=10_000; T=1, σ=0.1, Δ=1e-3, L=0.9, H=1.1)
σrΔ = σ * √Δ
steps = round(Int, T / Δ)
x = 1.
n_alive = sims
fr_inv = 2 / (σ^2 * Δ)
for i ∈ 1:sims
x_old = x
x_new = x_old
for t ∈ 1:steps
x_old = x_new
x_new += randn() * σrΔ
# direct exit
if !(L < x_new < H)
n_alive -= 1
break
end
# Brownian bridge: probability of crossing barrier between steps
p_low = exp(-(x_old - L) * (x_new - L) * fr_inv)
p_high = exp(-(H - x_old) * (H - x_new) * fr_inv)
if rand() < p_low + p_high
n_alive -= 1
break
end
end
end
n_alive / sims
end
gives with 1000 simulations
0.010199 seconds (3 allocations: 8.219 KiB) # barrier_sim_bb
0.006918 seconds # barrier_sim_bb_quick
and 10.000 simulations
0.126836 seconds (4 allocations: 79.531 KiB)
0.070636 seconds