Merge pull request #1 from augustinas1/main

Fix minor bugs

Author: kaandocal

src/build_rhs.jl

@@ -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
@@ -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
@@ -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
@@ -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
 
@@ -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
 

src/indexhandlers.jl

@@ -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`.
@@ -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)
 """
@@ -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.
@@ -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
 
@@ -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.
@@ -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
 
@@ -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
@@ -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
 
@@ -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
 
@@ -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
@@ -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)