Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
3 changes: 0 additions & 3 deletions .JuliaFormatter.toml

This file was deleted.

19 changes: 19 additions & 0 deletions .github/workflows/FormatCheck.yml
Original file line number Diff line number Diff line change
@@ -0,0 +1,19 @@
name: format-check

on:
push:
branches:
- 'master'
- 'main'
- 'release-'
tags: '*'
pull_request:

jobs:
runic:
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v4
- uses: fredrikekre/runic-action@v1
with:
version: '1'
2 changes: 1 addition & 1 deletion demo/birth_death.jl
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ prob = ODEProblem(sys, u0, 10.0, ps)

##

sol = solve(prob, Vern7(), dense = false, save_everystep = false, abstol = 1e-6)
sol = solve(prob, Vern7(), dense = false, save_everystep = false, abstol = 1.0e-6)

##

Expand Down
2 changes: 1 addition & 1 deletion demo/telegraph.jl
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ prob = ODEProblem(sys, u0, 10.0, ps)

##

sol = solve(prob, Vern7(), dense = false, save_everystep = false, abstol = 1e-6)
sol = solve(prob, Vern7(), dense = false, save_everystep = false, abstol = 1.0e-6)

##

Expand Down
14 changes: 9 additions & 5 deletions docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,8 @@ using Documenter
using FiniteStateProjection
using SparseArrays

makedocs(sitename = "FiniteStateProjection.jl",
makedocs(
sitename = "FiniteStateProjection.jl",
modules = [FiniteStateProjection],
format = Documenter.HTML(prettyurls = false),
pages = [
Expand All @@ -13,9 +14,12 @@ makedocs(sitename = "FiniteStateProjection.jl",
"Internal API" => "internal.md",
"Tips & Tricks" => "tips.md",
"Troubleshooting" => "troubleshoot.md",
"Examples" => "examples.md"
])
"Examples" => "examples.md",
]
)

deploydocs(repo = "github.com/SciML/FiniteStateProjection.jl.git",
deploydocs(
repo = "github.com/SciML/FiniteStateProjection.jl.git",
devbranch = "main",
push_preview = true)
push_preview = true
)
12 changes: 6 additions & 6 deletions src/build_rhs.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ See also: [`build_rhs_header`](@ref), [`build_rhs`](@ref)
function unpackparams(sys::FSPSystem, psym::Symbol)
param_names = Expr(:tuple, map(par -> par.name, Catalyst.parameters(sys.rs))...)

quote
return quote
$(param_names) = $(psym)
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)
quote
return quote
ps = p
$(unpackparams(sys, :ps))
end
Expand All @@ -50,7 +50,7 @@ function build_rhs_firstpass(sys::FSPSystem)
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
return quote
for idx_in in singleindices($(sys.ih), u)
$first_line
$(other_lines...)
Expand Down Expand Up @@ -106,7 +106,7 @@ function build_rhs_ex(sys::FSPSystem; striplines::Bool = false)

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

ex
return ex
end

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

##
Expand All @@ -141,5 +141,5 @@ calls `convert(ODEFunction, sys)`. It is usually more efficient to create an `OD
first and then use that to create `ODEProblem`s.
"""
function DiffEqBase.ODEProblem(sys::FSPSystem, u0, tint, pmap = NullParameters())
ODEProblem(convert(ODEFunction, sys), u0, tint, pmap_to_p(sys, pmap))
return ODEProblem(convert(ODEFunction, sys), u0, tint, pmap_to_p(sys, pmap))
end
8 changes: 4 additions & 4 deletions src/build_rhs_ss.jl
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ function build_rhs_ex_ss(sys::FSPSystem; striplines::Bool = false)

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

ex
return ex
end

"""
Expand All @@ -55,7 +55,7 @@ Builds the function `f(du,u,p,t)` that defines the right-hand side of the CME
for use with `SteadyStateProblem`s.
"""
function build_rhs_ss(sys::FSPSystem)
@RuntimeGeneratedFunction(build_rhs_ex_ss(sys; striplines = false))
return @RuntimeGeneratedFunction(build_rhs_ex_ss(sys; striplines = false))
end

##
Expand All @@ -67,7 +67,7 @@ Return an `ODEFunction` defining the right-hand side of the CME, for use
with `SteadyStateProblem`s.
"""
function Base.convert(::Type{ODEFunction}, sys::FSPSystem, ::SteadyState)
ODEFunction{true}(build_rhs_ss(sys))
return ODEFunction{true}(build_rhs_ss(sys))
end

"""
Expand All @@ -76,5 +76,5 @@ end
Return a `SteadyStateProblem` for use in `DifferentialEquations.
"""
function DiffEqBase.SteadyStateProblem(sys::FSPSystem, u0, pmap = NullParameters())
SteadyStateProblem(convert(ODEFunction, sys, SteadyState()), u0, pmap_to_p(sys, pmap))
return SteadyStateProblem(convert(ODEFunction, sys, SteadyState()), u0, pmap_to_p(sys, pmap))
end
30 changes: 19 additions & 11 deletions src/fspsystem.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,20 +9,22 @@ struct FSPSystem{IHT <: AbstractIndexHandler, RT}
rfs::RT
end

function FSPSystem(rs::ReactionSystem,
function FSPSystem(
rs::ReactionSystem,
ih::AbstractIndexHandler = DefaultIndexHandler{length(Catalyst.species(rs))}();
combinatoric_ratelaw::Bool = true)
combinatoric_ratelaw::Bool = true
)
isempty(Catalyst.get_systems(rs)) ||
error("Supported Catalyst models can not contain subsystems. Please use `rs = Catalyst.flatten(rs::ReactionSystem)` to generate a single system with no subsystems from your Catalyst model.")
any(eq -> !(eq isa Reaction), equations(rs)) &&
error("Catalyst models that include constraint ODEs or algebraic equations are not supported.")

rfs = create_ratefuncs(rs, ih; combinatoric_ratelaw = combinatoric_ratelaw)
FSPSystem(rs, ih, rfs)
return FSPSystem(rs, ih, rfs)
end

function FSPSystem(rs::ReactionSystem, order::AbstractVector{Symbol}; kwargs...)
FSPSystem(rs, PermutingIndexHandler(rs, order); kwargs...)
return FSPSystem(rs, PermutingIndexHandler(rs, order); kwargs...)
end

"""
Expand All @@ -33,8 +35,10 @@ Return the rate functions converted to Julia expressions in the state variable

See also: [`getsubstitutions`](@ref), [`build_rhs`](@ref)
"""
function build_ratefuncs(rs::ReactionSystem, ih::AbstractIndexHandler;
state_sym::Symbol, combinatoric_ratelaw::Bool = true)
function build_ratefuncs(
rs::ReactionSystem, ih::AbstractIndexHandler;
state_sym::Symbol, combinatoric_ratelaw::Bool = true
)
substitutions = getsubstitutions(ih, rs, state_sym = state_sym)

return map(Catalyst.reactions(rs)) do reac
Expand All @@ -47,18 +51,22 @@ end
function create_ratefuncs(rs::ReactionSystem, ih::AbstractIndexHandler; combinatoric_ratelaw::Bool = true)
paramsyms = Symbol.(Catalyst.parameters(rs))

return tuple(map(ex -> compile_ratefunc(ex, paramsyms),
build_ratefuncs(rs, ih; state_sym = :idx_in, combinatoric_ratelaw))...)
return tuple(
map(
ex -> compile_ratefunc(ex, paramsyms),
build_ratefuncs(rs, ih; state_sym = :idx_in, combinatoric_ratelaw)
)...
)
end

function compile_ratefunc(ex_rf, params)
# Make this nicer in the future
ex = :((idx_in, t, $(params...)) -> $(ex_rf)) |> MacroTools.flatten
@RuntimeGeneratedFunction(ex)
return @RuntimeGeneratedFunction(ex)
end

function pmap_to_p(sys::FSPSystem, pmap)
pmap_conv = eltype(pmap) <: Pair{Symbol} ? Catalyst.symmap_to_varmap(sys.rs, pmap) :
pmap
ModelingToolkit.varmap_to_vars(pmap_conv, Catalyst.parameters(sys.rs))
pmap
return ModelingToolkit.varmap_to_vars(pmap_conv, Catalyst.parameters(sys.rs))
end
50 changes: 32 additions & 18 deletions src/indexhandlers.jl
Original file line number Diff line number Diff line change
Expand Up @@ -79,31 +79,45 @@ DefaultIndexHandler{N}() where {N} = DefaultIndexHandler{N}(1, Tuple(1:N))
Base.vec(::DefaultIndexHandler, arr) = vec(arr)
Base.LinearIndices(::DefaultIndexHandler, arr) = LinearIndices(arr)

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

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

# Important: the species in `shift` are ordered according to `Catalyst.species`!
function pairedindices(ih::DefaultIndexHandler{N}, dims::NTuple{N, T},
shift::CartesianIndex{N}) where {N, T <: AbstractVector}
ranges = tuple((UnitRange(max(first(ax), first(ax)+shift[ih.perm[i]]),
min(last(ax), last(ax)+shift[ih.perm[i]]))
for (i, ax) in enumerate(dims))...)
function pairedindices(
ih::DefaultIndexHandler{N}, dims::NTuple{N, T},
shift::CartesianIndex{N}
) where {N, T <: AbstractVector}
ranges = tuple(
(
UnitRange(
max(first(ax), first(ax) + shift[ih.perm[i]]),
min(last(ax), last(ax) + shift[ih.perm[i]])
)
for (i, ax) in enumerate(dims)
)...
)

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

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

function pairedindices(::DefaultIndexHandler, dims::NTuple{N, T},
shift::CartesianIndex{M}) where {N, M, T <: AbstractVector}
@error "Dimension of state space ($(length(dims))) does not match number of species ($(length(shift)))"
function pairedindices(
::DefaultIndexHandler, dims::NTuple{N, T},
shift::CartesianIndex{M}
) where {N, M, T <: AbstractVector}
return @error "Dimension of state space ($(length(dims))) does not match number of species ($(length(shift)))"
end

"""
Expand All @@ -118,7 +132,7 @@ function getsubstitutions(ih::DefaultIndexHandler, rs::ReactionSystem; state_sym
species_orig = species(rs)
species_perm = [species_orig[ih.perm[i]] for i in 1:nspecs]

Dict(symbol => state_sym_vec[i] - ih.offset for (i, symbol) in enumerate(species_perm))
return Dict(symbol => state_sym_vec[i] - ih.offset for (i, symbol) in enumerate(species_perm))
end

#"""
Expand All @@ -128,7 +142,7 @@ end
#defined by the vector `order`.
#"""
function PermutingIndexHandler(rs::ReactionSystem, order::AbstractVector{Symbol})
PermutingIndexHandler(rs, map(sym -> Catalyst._symbol_to_var(rs, sym), order))
return PermutingIndexHandler(rs, map(sym -> Catalyst._symbol_to_var(rs, sym), order))
end

function PermutingIndexHandler(rs::ReactionSystem, order::AbstractVector)
Expand Down Expand Up @@ -158,5 +172,5 @@ function PermutingIndexHandler(rs::ReactionSystem, order::AbstractVector)

@assert count == ones(Int, nspec)

DefaultIndexHandler(1, Tuple(perm))
return DefaultIndexHandler(1, Tuple(perm))
end
8 changes: 4 additions & 4 deletions src/matrix.jl
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ function create_sparsematrix(sys::FSPSystem, dims::NTuple, ps, t)
end
end

sparse(I, J, V)
return sparse(I, J, V)
end

function create_sparsematrix_ss(sys::FSPSystem, dims::NTuple, ps)
Expand Down Expand Up @@ -72,7 +72,7 @@ function create_sparsematrix_ss(sys::FSPSystem, dims::NTuple, ps)
end
end

sparse(I, J, V)
return sparse(I, J, V)
end

"""
Expand All @@ -84,7 +84,7 @@ Chemical Master Equation. `dims` is a tuple denoting the dimensions of the FSP a
of the state obtained using `vec`.
"""
function SparseArrays.SparseMatrixCSC(sys::FSPSystem, dims::NTuple, pmap, t::Real)
create_sparsematrix(sys, dims, pmap_to_p(sys, pmap), t)
return create_sparsematrix(sys, dims, pmap_to_p(sys, pmap), t)
end

"""
Expand All @@ -94,5 +94,5 @@ Convert the reaction system into a sparse matrix defining the right-hand side of
Chemical Master Equation, steady-state version.
"""
function SparseArrays.SparseMatrixCSC(sys::FSPSystem, dims::NTuple, pmap, ::SteadyState)
create_sparsematrix_ss(sys, dims, pmap_to_p(sys, pmap))
return create_sparsematrix_ss(sys, dims, pmap_to_p(sys, pmap))
end
Loading
Loading