1+ # bracketing intervals for RSSA-based aggregators
2+ # see: "On the rejection-based algorithm for simulation and analysis of
3+ # large-scale reaction networks", Thanh et al, J. Chem. Phys., 2015
4+ # note, expects the type of the bracketing variables [ulow,uhigh] to be the
5+ # same as the fluct_rate and ushift.
16"""
27 BracketData(fluctrate, threshold, Δu)
38 BracketData{T1, T2}()
@@ -31,9 +36,9 @@ bd.fluctrate == 0.1
3136```
3237"""
3338struct BracketData{T1, T2}
34- fluctrate:: T1
35- threshold:: T2
36- Δu:: T2
39+ fluctrate:: T1 # interval should be [1-fluctrate,1+fluctrate] * u
40+ threshold:: T2 # for u below threshold interval is:
41+ Δu:: T2 # [max(u-Δu,0),u+Δu]
3742end
3843
3944# # suggested defaults
@@ -53,7 +58,7 @@ BracketData{T1, T2}() where {T1, T2} = BracketData(T1(0.1), T2(25), T2(4))
5358@inline getΔu (bd:: BracketData{T1, T2} , i) where {T1, T2 <: Number } = bd. Δu
5459
5560@inline function delta_bracket (u:: Integer , δ)
56- return (trunc (typeof (u), (one (δ) - δ) * u), trunc (typeof (u), (one (δ) + δ) * u))
61+ (trunc (typeof (u), (one (δ) - δ) * u), trunc (typeof (u), (one (δ) + δ) * u))
5762end
5863
5964@inline delta_bracket (u, δ) = ((one (δ) - δ) * u), ((one (δ) + δ) * u)
7580
7681# Get propensity brackets of massaction jump k.
7782@inline function get_majump_brackets (ulow, uhigh, k, majumps)
78- return evalrxrate (ulow, k, majumps), evalrxrate (uhigh, k, majumps)
83+ evalrxrate (ulow, k, majumps), evalrxrate (uhigh, k, majumps)
7984end
8085
8186# for constant rate jumps we must check the ordering of the bracket values
@@ -95,10 +100,8 @@ get brackets for the rate of reaction rx by first checking if the reaction is a
95100 if rx <= num_majumps
96101 return get_majump_brackets (p. ulow, p. uhigh, rx, ma_jumps)
97102 else
98- @inbounds return get_cjump_brackets (
99- p. ulow, p. uhigh, p. rates[rx - num_majumps],
100- params, t
101- )
103+ @inbounds return get_cjump_brackets (p. ulow, p. uhigh, p. rates[rx - num_majumps],
104+ params, t)
102105 end
103106end
104107
108111 @inbounds for (i, uval) in enumerate (u)
109112 ulow[i], uhigh[i] = get_spec_brackets (p. bracket_data, i, uval)
110113 end
111- return nothing
114+ nothing
112115end
113116
114117@inline function update_u_brackets! (p:: AbstractSSAJumpAggregator , u:: SVector )
@@ -117,12 +120,12 @@ end
117120 p. ulow = setindex (p. ulow, ulow, i)
118121 p. uhigh = setindex (p. uhigh, uhigh, i)
119122 end
120- return nothing
123+ nothing
121124end
122125
123126# For ExtendedJumpArray, only iterate over the species portion (u.u), not the jump tracking portion
124127@inline function update_u_brackets! (p:: AbstractSSAJumpAggregator , u:: ExtendedJumpArray )
125- return update_u_brackets! (p, u. u)
128+ update_u_brackets! (p, u. u)
126129end
127130
128131# Set up bracketing. The aggregator must have fields
@@ -151,5 +154,5 @@ function set_bracketing!(p::AbstractSSAJumpAggregator, u, params, t)
151154 end
152155 p. sum_rate = sum_rate
153156
154- return nothing
157+ nothing
155158end
0 commit comments