finish internal remake. Permit non-Float64 ints

Author: TorkelE

src/spatial_reaction_systems/lattice_reaction_systems.jl

@@ -29,17 +29,17 @@ struct LatticeReactionSystem{Q,R,S,T} <: MT.AbstractTimeDependentSystem
     All parameters related to the lattice reaction system
     (both those whose values are tied to vertices and edges).
     """
-    parameters::Vector{BasicSymbolic{Real}}
+    parameters::Vector{Any}
     """
     Parameters which values are tied to vertices, 
     e.g. that possibly could have unique values at each vertex of the system.
     """
-    vertex_parameters::Vector{BasicSymbolic{Real}}
+    vertex_parameters::Vector{Any}
     """
     Parameters whose values are tied to edges (adjacencies), 
     e.g. that possibly could have unique values at each edge of the system.
     """
-    edge_parameters::Vector{BasicSymbolic{Real}}
+    edge_parameters::Vector{Any}
     """
     An iterator over all the lattice's edges. Currently, the format is always a Vector{Pair{Int64,Int64}}.
     However, in the future, different types could potentially be used for different types of lattice

src/spatial_reaction_systems/spatial_ODE_systems.jl

@@ -62,7 +62,7 @@ struct LatticeTransportODEFunction{P,Q,R,S,T}
 
     function LatticeTransportODEFunction(ofunc::P, ps::Vector{<:Pair}, 
             lrs::LatticeReactionSystem, transport_rates::Vector{Pair{Int64, SparseMatrixCSC{S, Int64}}},
-            jac_transport::Union{Nothing, Matrix{S}, SparseMatrixCSC{S, Int64}}, sparse::Bool) where {P,S}        
+            jac_transport::Union{Nothing, Matrix{S}, SparseMatrixCSC{S, Int64}}) where {P,S}        
 
         # Creates a vector with the heterogeneous vertex parameters' indexes in the full parameter vector.
         p_dict = Dict(ps)
@@ -118,7 +118,7 @@ function (lt_ofun::LatticeTransportODEFunction)(du::AbstractVector, u, p, t)
         for e in lt_ofun.edge_iterator   
             idx_src = get_index(e[1], s, lt_ofun.num_species)
             idx_dst = get_index(e[2], s, lt_ofun.num_species)     
-            du[idx_dst] += get_transport_rate(s_idx, lt_ofun, e) * u[idx_src]
+            du[idx_dst] += get_transport_rate(rates, e, lt_ofun.t_rate_idx_types[s_idx]) * u[idx_src]
         end
     end
 end
@@ -182,10 +182,10 @@ function DiffEqBase.ODEProblem(lrs::LatticeReactionSystem, u0_in, tspan,
 end
 
 # Builds an ODEFunction for a spatial ODEProblem.
-function build_odefunction(lrs::LatticeReactionSystem, vert_ps::Vector{Pair{BasicSymbolic{Real},Vector{T}}},
-                           edge_ps::Vector{Pair{BasicSymbolic{Real},SparseMatrixCSC{T, Int64}}}, 
+function build_odefunction(lrs::LatticeReactionSystem, vert_ps::Vector{Pair{R,Vector{T}}},
+                           edge_ps::Vector{Pair{S,SparseMatrixCSC{T, Int64}}}, 
                            jac::Bool, sparse::Bool, name, include_zero_odes, combinatoric_ratelaws, 
-                           remove_conserved, checks) where {T}
+                           remove_conserved, checks) where {R,S,T}
     # Error check.
     if remove_conserved 
         error("Removal of conserved quantities is currently not supported for `LatticeReactionSystem`s")
@@ -209,7 +209,7 @@ function build_odefunction(lrs::LatticeReactionSystem, vert_ps::Vector{Pair{Basi
     end
 
     # Creates the `LatticeTransportODEFunction` functor (if `jac`, sets it as the Jacobian as well).
-    f = LatticeTransportODEFunction(ofunc_dense, [vert_ps; edge_ps], lrs, transport_rates, jac_transport, sparse)
+    f = LatticeTransportODEFunction(ofunc_dense, [vert_ps; edge_ps], lrs, transport_rates, jac_transport)
     J = (jac ? f : nothing)
 
     # Extracts the `Symbol` form for species and parameters. Creates and returns the `ODEFunction`.
@@ -268,7 +268,7 @@ function build_jac_prototype(ns_jac_prototype::SparseMatrixCSC{Float64, Int64},
     end
 
     # Create a sparse Jacobian prototype with 0-valued entries.
-    jac_prototype = sparse(i_idxs, j_idxs, zeros(num_entries))
+    jac_prototype = sparse(i_idxs, j_idxs, zeros(T, num_entries))
 
     # Set element values.
     if set_nonzero

src/spatial_reaction_systems/utility.jl

@@ -15,7 +15,7 @@ end
 
 # From u0 input, extract their values and store them in the internal format.
 # Internal format: a vector on the form [spec 1 at vert 1, spec 2 at vert 1, ..., spec 1 at vert 2, ...]).
-function lattice_process_u0(u0_in, u0_syms::Vector{BasicSymbolic{Real}}, lrs::LatticeReactionSystem)
+function lattice_process_u0(u0_in, u0_syms::Vector, lrs::LatticeReactionSystem)
     # u0 values can be given in various forms. This converts it to a Vector{Pair{Symbolics,...}} form.
     # Top-level vector: Maps each species to its value(s).
     u0 = lattice_process_input(u0_in, u0_syms) 
@@ -32,8 +32,8 @@ end
 
 # From a parameter input, split it into vertex parameters and edge parameters.
 # Store these in the desired internal format.
-function lattice_process_p(ps_in, ps_vertex_syms::Vector{BasicSymbolic{Real}}, 
-                           ps_edge_syms::Vector{BasicSymbolic{Real}}, lrs::LatticeReactionSystem)
+function lattice_process_p(ps_in, ps_vertex_syms::Vector, 
+                           ps_edge_syms::Vector, lrs::LatticeReactionSystem)
     # p values can be given in various forms. This converts it to a Vector{Pair{Symbolics,...}} form.
     # Top-level vector: Maps each parameter to its value(s).
     # Second-level: Contains either a vector (vertex parameters) or a sparse matrix (edge parameters).
@@ -52,7 +52,7 @@ end
 # The input (parameters or initial conditions) may either be a dictionary (symbolics to value(s).)
 # or a map (in vector or tuple form) from symbolics to value(s). This converts the input to a
 # (Vector) map from symbolics to value(s), where the entries have the same order as `syms`.
-function lattice_process_input(input::Dict{BasicSymbolic{Real}, T}, syms::Vector{BasicSymbolic{Real}}) where {T}
+function lattice_process_input(input::Dict{<:Any, T}, syms::Vector) where {T}
     # Error checks
     if !isempty(setdiff(keys(input), syms))
         error("You have provided values for the following unrecognised parameters/initial conditions: $(setdiff(keys(input), syms)).")
@@ -63,16 +63,15 @@ function lattice_process_input(input::Dict{BasicSymbolic{Real}, T}, syms::Vector
 
     return [sym => input[sym] for sym in syms]
 end
-function lattice_process_input(input, syms::Vector{BasicSymbolic{Real}})
+function lattice_process_input(input, syms::Vector)
     if ((input isa Vector) || (input isa Tuple)) && all(entry isa Pair for entry in input)
         return lattice_process_input(Dict(input), syms)
     end
     error("Input parameters/initial conditions have the wrong format ($(typeof(input))). These should either be a Dictionary, or a Tuple or a Vector (where each entry is a Pair taking a parameter/species to its value).")
 end
 
 # Splits parameters into vertex and edge parameters.
-# function split_parameters(ps::Vector{<: Pair}, p_vertex_syms::Vector, p_edge_syms::Vector)  
-function split_parameters(ps, p_vertex_syms::Vector{BasicSymbolic{Real}}, p_edge_syms::Vector{BasicSymbolic{Real}})  
+function split_parameters(ps, p_vertex_syms::Vector, p_edge_syms::Vector)  
     vert_ps = [p for p in ps if any(isequal(p[1]), p_vertex_syms)]     
     edge_ps = [p for p in ps if any(isequal(p[1]), p_edge_syms)]   
     return vert_ps, edge_ps
@@ -86,7 +85,7 @@ function vertex_value_map(values, lrs::LatticeReactionSystem)
 end
 
 # Converts the values for an individual species/vertex parameter to its correct vector form.
-function vertex_value_form(values, lrs::LatticeReactionSystem, sym::BasicSymbolic{Real})
+function vertex_value_form(values, lrs::LatticeReactionSystem, sym::BasicSymbolic)
     # If the value is a scalar (i.e. uniform across the lattice), return it in vector form.
     (values isa AbstractArray) || (return [values])
 
@@ -112,7 +111,7 @@ end
 
 # Converts values to the correct vector form for a Cartesian grid lattice.
 function vertex_value_form(values::AbstractArray, num_verts::Int64, lattice::CartesianGridRej{N,T}, 
-                           sym::BasicSymbolic{Real}) where {N,T}
+                           sym::BasicSymbolic) where {N,T}
     if size(values) != lattice.dims
         error("The values for $sym did not have the same format as the lattice. Expected a $(lattice.dims) array, got one of size $(size(values))")
     end
@@ -124,7 +123,7 @@ end
 
 # Converts values to the correct vector form for a masked grid lattice.
 function vertex_value_form(values::AbstractArray, num_verts::Int64, lattice::Array{Bool,T}, 
-                           sym::BasicSymbolic{Real}) where {T}
+                           sym::BasicSymbolic) where {T}
     if size(values) != size(lattice)
         error("The values for $sym did not have the same format as the lattice. Expected a $(size(lattice)) array, got one of size $(size(values))")
     end
@@ -174,9 +173,9 @@ end
 # The species is represented by its index (in species(lrs). 
 # If the rate is uniform across all edges, the transportation rate will be a size (1,1) sparse matrix.
 # Else, the rate will be a size (num_verts,num_verts) sparse matrix.
-function make_sidxs_to_transrate_map(vert_ps::Vector{Pair{BasicSymbolic{Real},Vector{T}}}, 
-                                     edge_ps::Vector{Pair{BasicSymbolic{Real},SparseMatrixCSC{T, Int64}}},
-                                     lrs::LatticeReactionSystem) where {T}
+function make_sidxs_to_transrate_map(vert_ps::Vector{Pair{R,Vector{T}}}, 
+                                     edge_ps::Vector{Pair{S,SparseMatrixCSC{T, Int64}}},
+                                     lrs::LatticeReactionSystem) where {R,S,T}
     # Creates a dictionary with each parameter's value(s).
     p_val_dict = Dict(vcat(vert_ps, edge_ps))
 
@@ -203,7 +202,7 @@ end
 # and the values of all our parameters, compute the transport rate(s).
 # If all parameters that the rate depends on are uniform across all edges, this becomes a length-1 vector.
 # Else it becomes a vector where each value corresponds to the rate at one specific edge.
-function compute_transport_rates(s::BasicSymbolic{Real}, p_val_dict, lrs::LatticeReactionSystem)
+function compute_transport_rates(s::BasicSymbolic, p_val_dict, lrs::LatticeReactionSystem)
     # Find parameters involved in the rate and create a function evaluating the rate law.
     rate_law = get_transport_rate_law(s, lrs)
     relevant_ps = Symbolics.get_variables(rate_law)
@@ -228,7 +227,7 @@ end
 # For a species, retrieve the symbolic expression for its transportation rate
 # (likely only a single parameter, such as `D`, but could be e.g. L*D, where L and D are parameters).
 # If there are several transportation reactions for the species, their sum is used.
-function get_transport_rate_law(s::BasicSymbolic{Real}, lrs::LatticeReactionSystem)
+function get_transport_rate_law(s::BasicSymbolic, lrs::LatticeReactionSystem)
     rates = filter(sr -> isequal(s, sr.species), spatial_reactions(lrs))
     return sum(getfield.(rates, :rate))
 end

test/spatial_modelling/lattice_reaction_systems.jl

@@ -223,10 +223,11 @@ let
     tr_3 = TransportReaction(dZ, Z)
     tr_macro_1 = @transport_reaction $dX X
     tr_macro_2 = @transport_reaction $(rate2) Y
+    @test_broken false
     # tr_macro_3 = @transport_reaction dZ $species3 # Currently does not work, something with meta programming.
 
     @test isequal(tr_1, tr_macro_1)
-    @test isequal(tr_2, tr_macro_2) # Unsure why these fails, since for components equality hold: `isequal(tr_1.species, tr_macro_1.species)` and `isequal(tr_1.rate, tr_macro_1.rate)` are both true.
+    @test isequal(tr_2, tr_macro_2) 
     # @test isequal(tr_3, tr_macro_3)
 end
 

test/spatial_modelling/lattice_reaction_systems_ODEs.jl

@@ -85,24 +85,23 @@ end
 
 ### Tests Simulation Correctness ###
 
-# Checks that non-spatial brusselator simulation is identical to all on an unconnected lattice.
+# Tests with non-Float64 parameter values.
 let
-    lrs = LatticeReactionSystem(brusselator_system, brusselator_srs_1, unconnected_graph)
-    u0 = [:X => 2.0 + 2.0 * rand(rng), :Y => 10.0 * (1.0 * rand(rng))]
-    pV = brusselator_p
-    pE = [:dX => 0.2]
-    oprob_nonspatial = ODEProblem(brusselator_system, u0, (0.0, 100.0), pV)
-    oprob_spatial = ODEProblem(lrs, u0, (0.0, 100.0), [pV; pE])
-    sol_nonspatial = solve(oprob_nonspatial, QNDF(); abstol = 1e-12, reltol = 1e-12)
-    sol_spatial = solve(oprob_spatial, QNDF(); abstol = 1e-12, reltol = 1e-12)
-
-    for i in 1:nv(unconnected_graph)
-        @test all(isapprox.(sol_nonspatial.u[end],
-                            sol_spatial.u[end][((i - 1) * 2 + 1):((i - 1) * 2 + 2)]))
+    lrs = LatticeReactionSystem(SIR_system, SIR_srs_2, very_small_2d_cartesian_grid)
+    u0 = [:S => 990.0, :I => rand_v_vals(lrs), :R => 0.0]
+    ps_1 = [:α => 0.1, :β => 0.01, :dS => 0.01, :dI => 0.01, :dR => 0.01]
+    ps_2 = [:α => 1//10, :β => 1//100, :dS => 1//100, :dI => 1//100, :dR => 1//100]
+    ps_3 = [:α => 1//10, :β => 0.01, :dS => 0.01, :dI => 1//100, :dR => 0.01]
+    sol_base = solve(ODEProblem(lrs, u0, (0.0, 100.0), ps_1), Rosenbrock23(); saveat = 0.1)
+    for ps in [ps_1, ps_2, ps_3]
+        for jac in [true, false], sparse in [true, false]
+            oprob = ODEProblem(lrs, u0, (0.0, 100.0), ps; jac, sparse)
+            @test sol_base ≈ solve(oprob, Rosenbrock23(); saveat = 0.1)
+        end
     end
 end
 
-# Compares Jacobian and forcing functions of spatial system to analytically computed on.
+# Compares Jacobian and forcing functions of spatial system to analytically computed ones.
 let
     # Creates LatticeReactionNetwork ODEProblem.
     rs = @reaction_network begin
@@ -119,7 +118,7 @@ let
     D_vals[1,2] = 0.2; D_vals[2,1] = 0.2; 
     D_vals[2,3] = 0.3; D_vals[3,2] = 0.3; 
     u0 = [:X => [1.0, 2.0, 3.0], :Y => 1.0]
-    ps = [:pX => [2.0, 2.5, 3.0], :pY => 0.5, :d => 0.1, :D => D_vals]
+    ps = [:pX => [2.0, 2.5, 3.0], :d => 0.1, :pY => 0.5, :D => D_vals]
     oprob = ODEProblem(lrs, u0, (0.0, 0.0), ps; jac=true, sparse=true)
 
     # Creates manual f and jac functions.
@@ -172,7 +171,7 @@ let
 
     # Sets test input values.
     u = rand(rng, 6)
-    p = [rand(rng, 3), rand(rng, 1), rand(rng, 1)]
+    p = [rand(rng, 3), ps[2][2], ps[3][2]]
 
     # Tests forcing function.
     du1 = fill(0.0, 6)
@@ -614,6 +613,22 @@ let
     end
 end
 
+# Tests with non-Int64 parameter values.
+let
+    lrs = LatticeReactionSystem(SIR_system, SIR_srs_2, very_small_2d_cartesian_grid)
+    u0 = [:S => 990.0, :I => rand_v_vals(lrs), :R => 0.0]
+    ps_1 = [:α => 0.1, :β => 0.01, :dS => 0.01, :dI => 0.01, :dR => 0.01]
+    ps_2 = [:α => Float32(0.1), :β => Float32(0.01), :dS => Float32(0.01), :dI => Float32(0.01), :dR => Float32(0.01)]
+    ps_3 = [:α => 1//10, :β => 0.01, :dS => 0.01, :dI => 1//100, :dR => Float32(0.01)]
+    sol_base = solve(ODEProblem(lrs, u0, (0.0, 100.0), ps_1), Rosenbrock23(); savetat = 0.1)
+    for ps in [ps_1, ps_2, ps_3]
+        for jac in [true, false], sparse in [true, false]
+            oprob = ODEProblem(lrs, u0, (0.0, 100.0), ps; jac, sparse)
+            @test sol_base ≈ solve(oprob, Rosenbrock23(); savetat = 0.1)
+        end
+    end
+end
+
 # Tests various types of numbers for initial conditions/parameters (e.g. Real numbers, Float32, etc.).
 let
     # Declare u0 versions.

test/spatial_modelling/lattice_reaction_systems_jumps.jl

@@ -64,7 +64,9 @@ let
 
     # Prepare various (diffusion) parameter input types.
     pE_1 = [:dI => 0.02, :dS => 0.01, :dR => 0.03]
-    pE_2 = [:dI => 0.02, :dS => uniform_e_vals(lrs, 0.01), :dR => 0.03]
+    dS_vals = spzeros(num_verts(lrs), num_verts(lrs))
+    foreach(e -> (dS_vals[e[1], e[2]] = 0.01), edge_iterator(lrs))
+    pE_2 = [:dI => 0.02, :dS => dS_vals, :dR => 0.03]
 
     # Checks hopping rates and u0 are correct.
     true_u0 = [fill(1.0, 1, 25); fill(2.0, 1, 25); fill(3.0, 1, 25)]