Merge remote-tracking branch 'refs/remotes/upstream/master'

Conflicts:
	src/spatial/topology.jl

Author: Vilin97

.github/workflows/Documentation.yml

@@ -22,4 +22,8 @@ jobs:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token
           DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key
           GKSwstype: "100" # https://discourse.julialang.org/t/generation-of-documentation-fails-qt-qpa-xcb-could-not-connect-to-display/60988
-        run: julia --project=docs/ docs/make.jl
+        run: julia --project=docs/ --code-coverage=user docs/make.jl
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v1
+        with:
+          file: lcov.info

Project.toml

@@ -1,7 +1,7 @@
 name = "JumpProcesses"
 uuid = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "9.1.0"
+version = "9.2.0"
 
 [deps]
 ArrayInterfaceCore = "30b0a656-2188-435a-8636-2ec0e6a096e2"
@@ -33,7 +33,7 @@ PoissonRandom = "0.4"
 RandomNumbers = "1.3"
 RecursiveArrayTools = "2"
 Reexport = "0.2, 1.0"
-SciMLBase = "1.35.1"
+SciMLBase = "1.51"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"
 TreeViews = "0.3"
 UnPack = "1.0.2"

src/coupling.jl

@@ -20,9 +20,13 @@ end
 function cat_problems(prob::DiffEqBase.AbstractODEProblem,
                       prob_control::DiffEqBase.AbstractODEProblem)
     l = length(prob.u0) # add l_c = length(prob_control.u0)
+
+    _f = SciMLBase.unwrapped_f(prob.f)
+    _f_control = SciMLBase.unwrapped_f(prob_control.f)
+
     new_f = function (du, u, p, t)
-        prob.f(@view(du[1:l]), u.u, p, t)
-        prob_control.f(@view(du[(l + 1):(2 * l)]), u.u_control, p, t)
+        _f(@view(du[1:l]), u.u, p, t)
+        _f_control(@view(du[(l + 1):(2 * l)]), u.u_control, p, t)
     end
     u0_coupled = CoupledArray(prob.u0, prob_control.u0, true)
     ODEProblem(new_f, u0_coupled, prob.tspan, prob.p)
@@ -33,9 +37,13 @@ function cat_problems(prob::DiscreteProblem, prob_control::DiffEqBase.AbstractOD
     if !(typeof(prob.f) <: typeof(DiffEqBase.DISCRETE_INPLACE_DEFAULT))
         @warn("Coupling to DiscreteProblem with nontrivial f. Note that, unless scale_by_time=true, the meaning of f will change when using an ODE/SDE/DDE/DAE solver.")
     end
+
+    _f = SciMLBase.unwrapped_f(prob.f)
+    _f_control = SciMLBase.unwrapped_f(prob_control.f)
+
     new_f = function (du, u, p, t)
-        prob.f(@view(du[1:l]), u.u, p, t)
-        prob_control.f(@view(du[(l + 1):(2 * l)]), u.u_control, p, t)
+        _f(@view(du[1:l]), u.u, p, t)
+        _f_control(@view(du[(l + 1):(2 * l)]), u.u_control, p, t)
     end
     u0_coupled = CoupledArray(prob.u0, prob_control.u0, true)
     ODEProblem(new_f, u0_coupled, prob.tspan, prob.p)
@@ -59,9 +67,13 @@ end
 function cat_problems(prob::DiffEqBase.AbstractSDEProblem,
                       prob_control::DiffEqBase.AbstractODEProblem)
     l = length(prob.u0)
+
+    _f = SciMLBase.unwrapped_f(prob.f)
+    _f_control = SciMLBase.unwrapped_f(prob_control.f)
+
     new_f = function (du, u, p, t)
-        prob.f(@view(du[1:l]), u.u, p, t)
-        prob_control.f(@view(du[(l + 1):(2 * l)]), u.u_control, p, t)
+        _f(@view(du[1:l]), u.u, p, t)
+        _f_control(@view(du[(l + 1):(2 * l)]), u.u_control, p, t)
     end
     new_g = function (du, u, p, t)
         prob.g(@view(du[1:l]), u.u, p, t)

src/problem.jl

@@ -227,9 +227,13 @@ function extend_problem(prob::DiffEqBase.AbstractDiscreteProblem, jumps; rng = D
 end
 
 function extend_problem(prob::DiffEqBase.AbstractODEProblem, jumps; rng = DEFAULT_RNG)
-    function jump_f(du::ExtendedJumpArray, u::ExtendedJumpArray, p, t)
-        prob.f(du.u, u.u, p, t)
-        update_jumps!(du, u, p, t, length(u.u), jumps.variable_jumps...)
+    _f = SciMLBase.unwrapped_f(prob.f)
+
+    jump_f = let _f = _f
+        function jump_f(du::ExtendedJumpArray, u::ExtendedJumpArray, p, t)
+            _f(du.u, u.u, p, t)
+            update_jumps!(du, u, p, t, length(u.u), jumps.variable_jumps...)
+        end
     end
     ttype = eltype(prob.tspan)
     u0 = ExtendedJumpArray(prob.u0,

src/spatial/topology.jl

@@ -109,7 +109,7 @@ struct CartesianGridRej{N, T}
     dims::NTuple{N, Int}
 
     "number of neighbor for each site"
-    nums_neighbors::Vector{Int}
+    nums_neighbors::Vector{Int8}
     CI::CartesianIndices{N, T}
     LI::LinearIndices{N, T}
 
@@ -127,7 +127,7 @@ function CartesianGridRej(dims::Tuple)
     CI = CartesianIndices(dims)
     LI = LinearIndices(dims)
     offsets = potential_offsets(dim)
-    nums_neighbors = [count(x -> x + CI[site] in CI, offsets) for site in 1:prod(dims)]
+    nums_neighbors = Int8[count(x -> x + CI[site] in CI, offsets) for site in 1:prod(dims)]
     CartesianGridRej(dims, nums_neighbors, CI, LI, offsets)
 end
 CartesianGridRej(dims) = CartesianGridRej(Tuple(dims))
@@ -144,6 +144,30 @@ function rand_nbr(rng, grid::CartesianGridRej, site::Int)
     end
 end
 
+# neighbor sampling is iterator-based
+struct CartesianGridIter{N, T}
+    dims::NTuple{N, Int}
+    nums_neighbors::Vector{Int8}
+    CI::CartesianIndices{N, T}
+    LI::LinearIndices{N, T}
+    offsets::Vector{CartesianIndex{N}}
+end
+function CartesianGridIter(dims::Tuple)
+    dim = length(dims)
+    CI = CartesianIndices(dims)
+    LI = LinearIndices(dims)
+    offsets = potential_offsets(dim)
+    nums_neighbors = Int8[count(x -> x + CI[site] in CI, offsets) for site in 1:prod(dims)]
+    CartesianGridIter(dims, nums_neighbors, CI, LI, offsets)
+end
+CartesianGridIter(dims) = CartesianGridIter(Tuple(dims))
+function CartesianGridIter(dimension, linear_size::Int)
+    CartesianGridIter([linear_size for i in 1:dimension])
+end
+function rand_nbr(rng, grid::CartesianGridIter, site::Int)
+    nth_nbr(grid, site, rand(rng, 1:outdegree(grid, site)))
+end
+
 function Base.show(io::IO, ::MIME"text/plain",
                    grid::CartesianGridRej)
     println(io, "A Cartesian grid with dimensions $(grid.dims)")

test/spatial/diffusion.jl

@@ -55,8 +55,8 @@ Bt = B'
 analytic_solution(t) = B * diagm(ℯ .^ (t * evals)) * Bt * reshape(prob.u0, num_nodes, 1)
 
 num_time_points = 10
-Nsims = 10000
-rel_tol = 0.01
+Nsims = 50000
+rel_tol = 0.02
 times = 0.0:(tf / num_time_points):tf
 
 algs = [NSM(), DirectCRDirect()]
@@ -143,7 +143,7 @@ for (j, spatial_jump_prob) in enumerate(jump_problems)
     end
 end
 
-# testing non-uniform hopping rates
+# testing non-uniform hopping rates. Only hopping up or left.
 dims = (2, 2)
 num_nodes = prod(dims)
 grid = Graphs.grid(dims)
@@ -153,7 +153,7 @@ for ci in CartesianIndices(hopping_constants)
     hopping_constants[species, site] = zeros(outdegree(grid, site))
     for (n, nb) in enumerate(neighbors(grid, site))
         if nb < site
-            hopping_constants[species, site][n] = 1.0
+            hopping_constants[species, site][n] = 10.0
         end
     end
 end