Skip to content

fix MassActionJump double scaling for non-unit substrate species in callbacks/remake #559

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
isaacsas merged 2 commits into from
Feb 27, 2026
Merged

fix MassActionJump double scaling for non-unit substrate species in callbacks/remake #559

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ecca871
Add rescale_rates_on_update field to MassActionJump to prevent double…
isaacsas Feb 26, 2026
44f5889
Fix PreScaledMapper test to actually catch double-scaling regression
isaacsas Feb 26, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
94 changes: 56 additions & 38 deletions src/jumps.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -333,10 +333,13 @@ struct MassActionJump{T, S, U, V} <: AbstractMassActionJump
net_stoch::U
"""Parameter mapping functor to identify reaction rate constants with parameters in `p` vectors."""
param_mapper::V
"""Whether `update_parameters!` should apply stoichiometric scaling to rates."""
rescale_rates_on_update::Bool

function MassActionJump{T, S, U, V}(rates::T, rs_in::S, ns::U, pmapper::V,
scale_rates::Bool, useiszero::Bool,
nocopy::Bool) where {T <: AbstractVector, S, U, V}
nocopy::Bool,
rescale_rates_on_update::Bool = scale_rates) where {T <: AbstractVector, S, U, V}
sr = nocopy ? rates : copy(rates)
rs = nocopy ? rs_in : copy(rs_in)
for i in eachindex(rs)
Expand All @@ -348,61 +351,67 @@ struct MassActionJump{T, S, U, V} <: AbstractMassActionJump
if scale_rates && !isempty(sr)
scalerates!(sr, rs)
end
new(sr, rs, ns, pmapper)
new(sr, rs, ns, pmapper, rescale_rates_on_update)
end
function MassActionJump{Nothing, Vector{S},
Vector{U}, V}(::Nothing, rs_in::Vector{S},
ns::Vector{U}, pmapper::V,
scale_rates::Bool,
useiszero::Bool,
nocopy::Bool) where {S <: AbstractVector,
nocopy::Bool,
rescale_rates_on_update::Bool = scale_rates) where {S <: AbstractVector,
U <: AbstractVector, V}
rs = nocopy ? rs_in : copy(rs_in)
for i in eachindex(rs)
if useiszero && (length(rs[i]) == 1) && iszero(rs[i][1][1])
rs[i] = typeof(rs[i])()
end
end
new(nothing, rs, ns, pmapper)
new(nothing, rs, ns, pmapper, rescale_rates_on_update)
end
function MassActionJump{T, S, U, V}(rate::T, rs_in::S, ns::U, pmapper::V,
scale_rates::Bool, useiszero::Bool,
nocopy::Bool) where {T <: Number, S, U, V}
nocopy::Bool,
rescale_rates_on_update::Bool = scale_rates) where {T <: Number, S, U, V}
rs = rs_in
if useiszero && (length(rs) == 1) && iszero(rs[1][1])
rs = typeof(rs)()
end
sr = scale_rates ? scalerate(rate, rs) : rate
new(sr, rs, ns, pmapper)
new(sr, rs, ns, pmapper, rescale_rates_on_update)
end
function MassActionJump{Nothing, S, U, V}(::Nothing, rs_in::S, ns::U, pmapper::V,
scale_rates::Bool, useiszero::Bool,
nocopy::Bool) where {S, U, V}
nocopy::Bool,
rescale_rates_on_update::Bool = scale_rates) where {S, U, V}
rs = rs_in
if useiszero && (length(rs) == 1) && iszero(rs[1][1])
rs = typeof(rs)()
end
new(nothing, rs, ns, pmapper)
new(nothing, rs, ns, pmapper, rescale_rates_on_update)
end
end
function MassActionJump(usr::T, rs::S, ns::U, pmapper::V; scale_rates = true,
useiszero = true, nocopy = false) where {T, S, U, V}
MassActionJump{T, S, U, V}(usr, rs, ns, pmapper, scale_rates, useiszero, nocopy)
useiszero = true, nocopy = false,
rescale_rates_on_update = scale_rates) where {T, S, U, V}
MassActionJump{T, S, U, V}(usr, rs, ns, pmapper, scale_rates, useiszero, nocopy,
rescale_rates_on_update)
end
function MassActionJump(usr::T, rs, ns; scale_rates = true, useiszero = true,
nocopy = false) where {T <: AbstractVector}
MassActionJump(usr, rs, ns, nothing; scale_rates = scale_rates, useiszero = useiszero,
nocopy = nocopy)
nocopy = false, rescale_rates_on_update = scale_rates) where {T <: AbstractVector}
MassActionJump(usr, rs, ns, nothing; scale_rates, useiszero, nocopy,
rescale_rates_on_update)
end
function MassActionJump(usr::T, rs, ns; scale_rates = true, useiszero = true,
nocopy = false) where {T <: Number}
MassActionJump(usr, rs, ns, nothing; scale_rates = scale_rates, useiszero = useiszero,
nocopy = nocopy)
nocopy = false, rescale_rates_on_update = scale_rates) where {T <: Number}
MassActionJump(usr, rs, ns, nothing; scale_rates, useiszero, nocopy,
rescale_rates_on_update)
end

# with parameter indices or mapping, multiple jump case
function MassActionJump(rs, ns; param_idxs = nothing, param_mapper = nothing,
scale_rates = true, useiszero = true, nocopy = false)
scale_rates = true, useiszero = true, nocopy = false,
rescale_rates_on_update = scale_rates)
if param_mapper === nothing
(param_idxs === nothing) &&
error("If no parameter indices are given via param_idxs, an explicit parameter mapping must be passed in via param_mapper.")
Expand All @@ -413,8 +422,8 @@ function MassActionJump(rs, ns; param_idxs = nothing, param_mapper = nothing,
pmapper = param_mapper
end

MassActionJump(nothing, nocopy ? rs : copy(rs), ns, pmapper; scale_rates = scale_rates,
useiszero = useiszero, nocopy = true)
MassActionJump(nothing, nocopy ? rs : copy(rs), ns, pmapper; scale_rates, useiszero,
nocopy = true, rescale_rates_on_update)
end

using_params(maj::MassActionJump{T, S, U, Nothing}) where {T, S, U} = false
Expand Down Expand Up @@ -470,18 +479,19 @@ function Base.merge(pmap1::MassActionJumpParamMapper{Int},
end

"""
update_parameters!(maj::MassActionJump, newparams; scale_rates=true)
update_parameters!(maj::MassActionJump, newparams; scale_rates=maj.rescale_rates_on_update)

Updates the passed in MassActionJump with the parameter values in `newparams`.

Notes:

- Requires the jump to have been constructed with a user-passed `param_idxs` or `param_mapper`.
- `scale_rates=true` will scale the parameter representing the jump rate by an
appropriate combinatoric factor. i.e for 3A --> B at rate k it will scale
k --> k/3!.
- `scale_rates` defaults to `maj.rescale_rates_on_update`, which itself defaults to the
`scale_rates` value used when constructing the jump. When `true`, the parameter
representing the jump rate will be scaled by an appropriate combinatoric factor, i.e.
for 3A --> B at rate k it will scale k --> k/3!.
"""
function update_parameters!(maj::MassActionJump, newparams; scale_rates = true, kwargs...)
function update_parameters!(maj::MassActionJump, newparams; scale_rates = maj.rescale_rates_on_update, kwargs...)
(maj.param_mapper === nothing) &&
error("MassActionJumps must be constructed with param_idxs or a param_mapper to be updateable.")
maj.param_mapper(maj, newparams; scale_rates, kwargs)
Expand Down Expand Up @@ -556,8 +566,12 @@ function JumpSet(vjs, cjs, rj, majv::Vector{T}) where {T <: MassActionJump}
error("JumpSets do not accept empty mass action jump collections; use \"nothing\" instead.")
end

sr_val = majv[1].rescale_rates_on_update
if !all(m -> m.rescale_rates_on_update == sr_val, majv)
error("Cannot merge MassActionJumps with different rescale_rates_on_update settings.")
end
maj = setup_majump_to_merge(majv[1].scaled_rates, majv[1].reactant_stoch,
majv[1].net_stoch, majv[1].param_mapper)
majv[1].net_stoch, majv[1].param_mapper, sr_val)
for i in 2:length(majv)
massaction_jump_combine(maj, majv[i])
end
Expand Down Expand Up @@ -618,35 +632,35 @@ end
# functionality to merge two mass action jumps together
function check_majump_type(maj::MassActionJump{S, T, U, V}) where {S <: Number, T, U, V}
setup_majump_to_merge(maj.scaled_rates, maj.reactant_stoch, maj.net_stoch,
maj.param_mapper)
maj.param_mapper, maj.rescale_rates_on_update)
end
function check_majump_type(maj::MassActionJump{Nothing, T, U, V}) where {T, U, V}
setup_majump_to_merge(maj.scaled_rates, maj.reactant_stoch, maj.net_stoch,
maj.param_mapper)
maj.param_mapper, maj.rescale_rates_on_update)
end

# if given containers of rates and stoichiometry directly create a jump
function setup_majump_to_merge(sr::T, rs::AbstractVector{S}, ns::AbstractVector{U},
pmapper) where {T <: AbstractVector, S <: AbstractArray,
function setup_majump_to_merge(sr::T, rs::AbstractVector{S}, ns::AbstractVector{U}, pmapper,
rescale_rates_on_update::Bool) where {T <: AbstractVector, S <: AbstractArray,
U <: AbstractArray}
MassActionJump(sr, rs, ns, pmapper; scale_rates = false)
MassActionJump(sr, rs, ns, pmapper; scale_rates = false, rescale_rates_on_update)
end

# if just given the data for one jump (and not in a container) wrap in a vector
function setup_majump_to_merge(sr::S, rs::T, ns::U,
pmapper) where {S <: Number, T <: AbstractArray,
function setup_majump_to_merge(sr::S, rs::T, ns::U, pmapper,
rescale_rates_on_update::Bool) where {S <: Number, T <: AbstractArray,
U <: AbstractArray}
MassActionJump([sr], [rs], [ns],
(pmapper === nothing) ? pmapper : to_collection(pmapper);
scale_rates = false)
scale_rates = false, rescale_rates_on_update)
end

# if no rate field setup yet
function setup_majump_to_merge(::Nothing, rs::T, ns::U,
pmapper) where {T <: AbstractArray, U <: AbstractArray}
function setup_majump_to_merge(::Nothing, rs::T, ns::U, pmapper,
rescale_rates_on_update::Bool) where {T <: AbstractArray, U <: AbstractArray}
MassActionJump(nothing, [rs], [ns],
(pmapper === nothing) ? pmapper : to_collection(pmapper);
scale_rates = false)
scale_rates = false, rescale_rates_on_update)
end

# when given a collection of reactions to add to maj
Expand Down Expand Up @@ -697,17 +711,21 @@ function majump_merge!(maj::MassActionJump{T, S, U, V}, sr::T, rs::S, ns::U,
(param_mapper === nothing) ||
error("Error, trying to merge a MassActionJump with a parameter mapping to one without a parameter mapping.")
return MassActionJump(rates, [maj.reactant_stoch, rs], [maj.net_stoch, ns],
param_mapper; scale_rates = false)
param_mapper; scale_rates = false,
rescale_rates_on_update = maj.rescale_rates_on_update)
else
return MassActionJump(rates, [maj.reactant_stoch, rs], [maj.net_stoch, ns],
merge(maj.param_mapper, param_mapper); scale_rates = false)
merge(maj.param_mapper, param_mapper); scale_rates = false,
rescale_rates_on_update = maj.rescale_rates_on_update)
end
end

massaction_jump_combine(maj1::MassActionJump, maj2::Nothing) = maj1
massaction_jump_combine(maj1::Nothing, maj2::MassActionJump) = maj2
massaction_jump_combine(maj1::Nothing, maj2::Nothing) = maj1
function massaction_jump_combine(maj1::MassActionJump, maj2::MassActionJump)
(maj1.rescale_rates_on_update == maj2.rescale_rates_on_update) ||
error("Cannot merge MassActionJumps with different rescale_rates_on_update settings.")
majump_merge!(maj1, maj2.scaled_rates, maj2.reactant_stoch, maj2.net_stoch,
maj2.param_mapper)
end
Expand Down
4 changes: 3 additions & 1 deletion src/spatial/spatial_massaction_jump.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,9 @@ function SpatialMassActionJump(urates::A, rs, ns; scale_rates = true, useiszero
useiszero = useiszero, nocopy = nocopy)
end

function SpatialMassActionJump(ma_jumps::MassActionJump{T, S, U, V}; scale_rates = true,
# scale_rates defaults to false since ma_jumps.scaled_rates are already scaled;
# passing true would double-scale.
function SpatialMassActionJump(ma_jumps::MassActionJump{T, S, U, V}; scale_rates = false,
useiszero = true, nocopy = false) where {T, S, U, V}
SpatialMassActionJump(ma_jumps.scaled_rates, ma_jumps.reactant_stoch,
ma_jumps.net_stoch, ma_jumps.param_mapper;
Expand Down
1 change: 1 addition & 0 deletions test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@ end
@time @safetestset "Mass Action Jump Tests; Nonlinear Rx Model" begin include("bimolerx_test.jl") end
@time @safetestset "Mass Action Jump Tests; Special Cases" begin include("degenerate_rx_cases.jl") end
@time @safetestset "Mass Action Jump Tests; Floating Point Inputs" begin include("fp_unknowns.jl") end
@time @safetestset "scale_rates Field Tests" begin include("scale_rates_field_test.jl") end
@time @safetestset "Direct allocations test" begin include("allocations.jl") end
@time @safetestset "Bracketing Tests" begin include("bracketing.jl") end
@time @safetestset "Composition-Rejection Table Tests" begin include("table_test.jl") end
Expand Down
Loading
Loading