Skip to content

Fix minor bugs #1

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
merged 2 commits into from
Jun 29, 2021
Merged

Fix minor bugs #1

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
323586d
fix `convert`
augustinas1 Jun 24, 2021
bf8d34e
fix ReducingIndexHandler
augustinas1 Jun 24, 2021
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
50 changes: 25 additions & 25 deletions src/build_rhs.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,8 +10,8 @@ See also: [`build_rhs_header`](@ref), [`build_rhs`](@ref)
"""
function unpackparams(sys::FSPSystem, psym::Symbol)::Expr
param_names = Expr(:tuple, map(par -> par.name, params(sys.rs))...)
quote

quote
$(param_names) = ps
end
end
Expand All @@ -26,7 +26,7 @@ just unpacks parameters from `p`.
See also: [`unpackparams`](@ref), [`build_rhs`](@ref)
"""
function build_rhs_header(sys::FSPSystem)::Expr
quote
quote
ps::AbstractVector{Float64} = p
$(unpackparams(sys, :ps))
end
Expand All @@ -38,18 +38,18 @@ end
build_rhs_firstpass(sys::FSPSystem, rfs)::Expr

Return code for the first pass of the RHS function, for the time-dependent
FSP. Goes through all reactions and computes the negative part of the CME
(probability flowing out of states). This is a simple array traversal and
FSP. Goes through all reactions and computes the negative part of the CME
(probability flowing out of states). This is a simple array traversal and
can be done in one go for all reactions.

See also: [`build_rhs`](@ref)
"""
function build_rhs_firstpass(sys::FSPSystem)
isempty(sys.rfs) && return quote end
isempty(sys.rfs) && return quote end

first_line = :(du[idx_in] = -u[idx_in] * $(sys.rfs[1].body))
other_lines = (:(du[idx_in] -= u[idx_in] * $(rf.body)) for rf in sys.rfs[2:end])

quote
for idx_in in singleindices($(sys.ih), u)
$first_line
Expand All @@ -66,55 +66,55 @@ end
Return code for the second pass of the RHS function. Goes through
all reactions and computes the positive part of the CME (probability
flowing into states). This requires accessing `du` and `u` at different
locations depending on the net stoichiometries. In order to reduce
locations depending on the net stoichiometries. In order to reduce
random memory access reactions are processed one by one.

See also: [`build_rhs`](@ref)
"""
function build_rhs_secondpass(sys::FSPSystem)::Expr
isempty(sys.rfs) && return quote end

S = netstoichmat(sys.rs)
ret = Expr(:block)

for (i, rf) in enumerate(sys.rfs)
ex = quote
for (idx_in, idx_out) in pairedindices($(sys.ih), u, $(CartesianIndex(S[i,:]...)))
du[idx_out] += u[idx_in] * $(rf.body)
end
end

append!(ret.args, ex.args)
end

return ret
end

##

function build_rhs_ex(sys::FSPSystem; striplines::Bool=true)
function build_rhs_ex(sys::FSPSystem; striplines::Bool=false)
header = build_rhs_header(sys)
first_pass = build_rhs_firstpass(sys)
second_pass = build_rhs_secondpass(sys)

body = Expr(:block, header, first_pass, second_pass)
ex = :((du, u, p, t) -> $(body))

ex = :((du, u, p, t) -> $(body))

striplines && (ex = MacroTools.striplines(ex))

ex = ex |> MacroTools.flatten |> MacroTools.prettify

ex
end

"""
build_rhs(sys::FSPSystem)

Builds the function `f(du,u,p,t)` that defines the right-hand side of the CME
for use with DifferentialEquations.jl.
for use with DifferentialEquations.jl.
"""
function build_rhs(sys::FSPSystem)
function build_rhs(sys::FSPSystem)
@RuntimeGeneratedFunction(build_rhs_ex(sys; striplines=false))
end

Expand All @@ -126,14 +126,14 @@ end

Return an `ODEFunction` defining the right-hand side of the CME.

Creates an `ODEFunction` for use with `DifferentialEquations`. This is where most of
the work in the package happens; for best performance it is suggested to build an `ODEFunction`
Creates an `ODEFunction` for use with `DifferentialEquations`. This is where most of
the work in the package happens; for best performance it is suggested to build an `ODEFunction`
once for a given reaction system and reuse it instead of directly converting
a reaction system to an `ODEProblem` (which implicitly calls this function).
"""
function Base.convert(::Type{ODEFunction}, sys::FSPSystem;
combinatoric_ratelaw::Bool=true)::ODEFunction
rhs = build_rhs(sys; striplines=false)
rhs = build_rhs(sys)
ODEFunction{true}(rhs)
end

Expand Down
74 changes: 37 additions & 37 deletions src/indexhandlers.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,13 +1,13 @@
"""
"""
singleindices(idxhandler::AbstractIndexHandler, arr)

Returns all indices `I` in `arr`. Defaults to CartesianIndices, but can
be overloaded for arbitrary index handlers.
be overloaded for arbitrary index handlers.
"""
singleindices(::AbstractIndexHandler, arr::AbstractArray) = CartesianIndices(arr)
singleindices(::AbstractIndexHandler, arr::Tuple) = CartesianIndices(arr)

"""
"""
pairedindices(idxhandler::AbstractIndexHandler, arr, shift::CartesianIndex)

Returns all pairs of indices `(I .- shift, I)` in `arr`.
Expand All @@ -17,16 +17,16 @@ function pairedindices end
"""
getsubstitutions(idxhandler::AbstractIndexHandler, rs::ReactionSystem; state_sym::Symbol)

Returns a dict of the form `S_i => f_i(state_sym)`, where each `f_i` is an expression
Returns a dict of the form `S_i => f_i(state_sym)`, where each `f_i` is an expression
for the abundance of species `S_i` in terms of the state variable `state_sym`.
"""
function getsubstitutions end

"""
vec(idxhandler::AbstractIndexHandler, arr)

Converts the right-hand side defining the solution of the CME into a
one-dimensional vector to which a matrix can be applied.
Converts the right-hand side defining the solution of the CME into a
one-dimensional vector to which a matrix can be applied.

See also: [`LinearIndices`](@ref Base.LinearIndices)
"""
Expand All @@ -50,15 +50,15 @@ function LinearIndices end
##


"""
"""
struct NaiveIndexHandler <: AbstractIndexHandler
offset::Int
end

Basic index handler that stores the state of a system with
`s` species in an `s`-dimensional array. The `offset` parameter
denotes the offset by which the array is indexed (defaults to 1
in Julia).
in Julia).

This is the simplest index handler, but it will not be optimal
if some states cannot be reached from the initial state, e.g.
Expand All @@ -79,24 +79,24 @@ NaiveIndexHandler() = NaiveIndexHandler(1)
Base.vec(::NaiveIndexHandler, arr) = vec(arr)
Base.LinearIndices(::NaiveIndexHandler, arr) = LinearIndices(arr)

function pairedindices(ih::NaiveIndexHandler, arr::AbstractArray{T,N},
function pairedindices(ih::NaiveIndexHandler, arr::AbstractArray{T,N},
shift::CartesianIndex{N}) where {T,N}
pairedindices(ih, axes(arr), shift)
end

function pairedindices(ih::NaiveIndexHandler, dims::NTuple{N,T},
function pairedindices(ih::NaiveIndexHandler, dims::NTuple{N,T},
shift::CartesianIndex{N}) where {N,T<:Number}
pairedindices(ih, Base.OneTo.(dims), shift)
end

function pairedindices(::NaiveIndexHandler, dims::NTuple{N,T},
function pairedindices(::NaiveIndexHandler, dims::NTuple{N,T},
shift::CartesianIndex{N}) where {N,T<:AbstractVector}
ranges = tuple((UnitRange(max(first(ax), first(ax)+shift[i]),
min(last(ax), last(ax)+shift[i]))
ranges = tuple((UnitRange(max(first(ax), first(ax)+shift[i]),
min(last(ax), last(ax)+shift[i]))
for (i, ax) in enumerate(dims))...)

ranges_shifted = tuple((rng .- shift[i] for (i, rng) in enumerate(ranges))...)

zip(CartesianIndices(ranges_shifted), CartesianIndices(ranges))
end

Expand All @@ -111,13 +111,13 @@ end

##

"""
"""
struct ReducingIndexHandler <: AbstractIndexHandler

More efficient index handler that improves upon [`NaiveIndexHandler`](@ref)
by eliminating variables whose abundances can be computed from other variables
using conservation laws. Describes the system using a subset of the original
species which can be obtained via [`reducedspecies`](@ref). Reduces the
species which can be obtained via [`reducedspecies`](@ref). Reduces the
dimensionality of the FSP by the number of conservation laws in the system.

Note: This feature is currently being moved into Catalyst.jl.
Expand All @@ -135,12 +135,12 @@ struct ReducingIndexHandler <: AbstractIndexHandler
end

function ReducingIndexHandler(rs::ReactionSystem, offset::Int=1)
cons_laws = conservationlaws(sys)
cons_laws = conservationlaws(rs)
specs_elided = elidedspecies(cons_laws)
specs_red = [ i for i in 1:length(species(sys.rs)) if !(i in specs_elided) ]
specs_red = [ i for i in 1:length(species(rs)) if !(i in specs_elided) ]

cons_syms = [ gensym("c$i") for i in 1:size(cons_laws, 1) ]

return ReducingIndexHandler(cons_laws, offset, specs_red, specs_elided, cons_syms)
end

Expand Down Expand Up @@ -170,7 +170,7 @@ elidedspecies(idxhandler::ReducingIndexHandler) = idxhandler.specs_elided

##

"""
"""
elidedspecies(cons_laws::AbstractMatrix{Int})::Vector

Returns a list of species ``[ s_1, ... ]`` which can be removed from the reaction system
Expand All @@ -186,9 +186,9 @@ function elidedspecies(cons_laws::AbstractMatrix{Int})::Vector{Int}
idx = possible_specs[end]
ret[i] = idx
end

@assert length(ret) == length(unique(ret))

return ret
end

Expand All @@ -204,19 +204,19 @@ See also: [`getsubstitutions`](@ref)
function elisions(idxhandler::ReducingIndexHandler, rs::ReactionSystem)
ret = Dict()
spec_syms = species(rs)

for (i, spec) in enumerate(elidedspecies(idxhandler))
sym = spec_syms[spec]
cons_law = idxhandler.cons_laws[i,:]

# Does this always work? What if some of the species on the RHS
# also end up getting elided at some point?
rhs = (Variable(idxhandler.cons_syms[i]) - sum(cons_law[j] * spec_syms[j] for j in 1:length(spec_syms) if j != spec))
rhs /= cons_law[i]

ret[sym] = rhs
end

ret
end

Expand All @@ -230,28 +230,28 @@ from the conserved quantities and the reduced species.
"""
function getsubstitutions(idxhandler::ReducingIndexHandler, rs::ReactionSystem; state_sym::Symbol)
symbols = states(rs)
ret = Dict{Any,Any}(symbols[spec] => Term(Base.getindex, (state_sym, i)) - idxhandler.offset

ret = Dict{Any,Any}(symbols[spec] => Term(Base.getindex, (state_sym, i)) - idxhandler.offset
for (i, spec) in enumerate(reducedspecies(idxhandler)))

elision_dict = elisions(idxhandler, rs)
for (spec, ex) in pairs(elision_dict)
ret[spec] = substitute(ex, ret)
end

ret
end

"""
build_rhs_header(idxhandler::ReducingIndexHandler, sys::FSPSystem)::Expr

Assumes `p` is of the form `(params, cons::AbstractVector{Int})` where `params`
Assumes `p` is of the form `(params, cons::AbstractVector{Int})` where `params`
are the system parameters and `cons` the conserved quantities.
"""
function build_rhs_header(sys::FSPSystem{ReducingIndexHandler})::Expr
cons_names = Expr(:tuple, sys.ih.cons_syms...)
quote

quote
(ps, $(cons_names)) = p
$(unpackparams(sys, :ps))
end
Expand All @@ -263,13 +263,13 @@ end
Similar to its `NaiveIndexHandler` variant, but converts the indices into indices into
the reduced state space array.
"""
function pairedindices(idxhandler::ReducingIndexHandler, arr::AbstractArray{T,M},
function pairedindices(idxhandler::ReducingIndexHandler, arr::AbstractArray{T,M},
shift::CartesianIndex{N}) where {T,M,N}
shift_red = CartesianIndex{M}(convert(Tuple, shift)[reducedspecies(idxhandler)]...)
pairedindices(NaiveIndexHandler(idxhandler.offset), arr, shift_red)
end

function pairedindices(idxhandler::ReducingIndexHandler, dims::NTuple{M},
function pairedindices(idxhandler::ReducingIndexHandler, dims::NTuple{M},
shift::CartesianIndex{N}) where {M,N}
shift_red = CartesianIndex{M}(convert(Tuple, shift)[reducedspecies(idxhandler)]...)
pairedindices(NaiveIndexHandler(idxhandler.offset), dims, shift_red)
Expand Down