DirectCR-RSSA

[Vilin97]

A new spatial SSA called DirectCRRSSA uses RSSA-style sampling to sample the site+jump. It is about 20% faster than DirectCRDirect, at least on the problem in the test/spatial/ABC.jl (5 sites in a line, 3 species, 2 reactions).

using BenchmarkTools
@btime solve(jump_problems[1], SSAStepper()) # NSM with a `CartesianGrid`, 52ms
@btime solve(jump_problems[2], SSAStepper()) # NSM with a `Graph`, 49ms
@btime solve(jump_problems[3], SSAStepper()) # DirectCRDirect with  a `CartesianGrid`, 51ms
@btime solve(jump_problems[4], SSAStepper()) # DirectCRRSSA with a `CartesianGrid`, 42ms

[github-actions]

Pull Request Test Coverage Report for Build 6141985002

• 172 of 181 (95.03%) changed or added relevant lines in 8 files are covered.
• 12 unchanged lines in 5 files lost coverage.
• Overall coverage decreased (-0.1%) to 89.759%

---

Changes Missing Coverage	Covered Lines	Changed/Added Lines	%
src/spatial/bracketing.jl	35	38	92.11%
src/spatial/rssacrdirect.jl	105	111	94.59%

Files with Coverage Reduction	New Missed Lines	%
src/aggregators/prioritytable.jl	1	86.86%
src/extended_jump_array.jl	1	80.52%
src/coupling.jl	2	78.18%
src/spatial/spatial_massaction_jump.jl	3	85.45%
src/spatial/directcrdirect.jl	5	87.65%

Totals
Change from base Build 6133295964:	-0.1%
Covered Lines:	2419
Relevant Lines:	2695

---

💛 - Coveralls

[Vilin97]

There was a very hard-to-find bug in the ordering of arguments in generate_jumps!. I fixed it here and in NSM.jl.

[coveralls]

Pull Request Test Coverage Report for Build 6180696485

• 155 of 161 (96.27%) changed or added relevant lines in 8 files are covered.
• No unchanged relevant lines lost coverage.
• Overall coverage increased (+0.3%) to 90.175%

---

Changes Missing Coverage	Covered Lines	Changed/Added Lines	%
src/spatial/bracketing.jl	23	24	95.83%
src/spatial/directcrrssa.jl	106	111	95.5%

Totals
Change from base Build 6167199515:	0.3%
Covered Lines:	2423
Relevant Lines:	2687

---

💛 - Coveralls

[isaacsas]

This needs rebasing too it appears.

[isaacsas]

Does switching to using a branch have performance implications? We create these structures a lot right in the scalar case right? I'm not sure if Julia will remove it during compilation.

[Vilin97]

Even if it's not optimized away, it shouldn't be expensive, right? I am not good at reading low level code but it seems that it does get optimized away in the end.

# Old struct without branching
struct LowHighOld{T}
    low::T
    high::T

    function LowHighOld(low::T, high::T) where {T}
        new{T}(deepcopy(low), deepcopy(high))
    end
end

# New struct with branching
struct LowHighNew{T}
    low::T
    high::T

    function LowHighNew(low::T, high::T; do_copy=true) where {T}
        if do_copy
            return new{T}(deepcopy(low), deepcopy(high))
        else
            return new{T}(low, high)
        end
    end
end

# Test values
low_value = 10
high_value = 20

# Code introspection to see if the branch is removed in scalar cases
@code_llvm debuginfo=:none LowHighOld(low_value, high_value)

@code_llvm debuginfo=:none LowHighNew(low_value, high_value; do_copy=true)

The output of the code introspection is as follows:

define void @julia_LowHighOld_2364([2 x i64]* noalias nocapture noundef nonnull sret([2 x i64]) 
align 8 dereferenceable(16) %0, {}* noundef nonnull readonly %1, i64 signext %2, i64 signext %3) #0 {
top:
  %newstruct.sroa.0.0..sroa_idx = getelementptr inbounds [2 x i64], [2 x i64]* %0, i64 0, i64 0 
  store i64 %2, i64* %newstruct.sroa.0.0..sroa_idx, align 8
  %newstruct.sroa.2.0..sroa_idx1 = getelementptr inbounds [2 x i64], [2 x i64]* %0, i64 0, i64 1  store i64 %3, i64* %newstruct.sroa.2.0..sroa_idx1, align 8
  ret void
}

define void @julia_LowHighNew_2366([2 x i64]* noalias nocapture noundef nonnull sret([2 x i64]) 
align 8 dereferenceable(16) %0, [1 x i8]* nocapture noundef nonnull readonly align 1 dereferenceable(1) %1, {}* noundef nonnull readonly %2, i64 signext %3, i64 signext %4) #0 {
top:
  %.sroa.025.0..sroa_idx = getelementptr inbounds [2 x i64], [2 x i64]* %0, i64 0, i64 0        
  store i64 %3, i64* %.sroa.025.0..sroa_idx, align 8
  %.sroa.2.0..sroa_idx26 = getelementptr inbounds [2 x i64], [2 x i64]* %0, i64 0, i64 1
  store i64 %4, i64* %.sroa.2.0..sroa_idx26, align 8
  ret void
}

[isaacsas]

Suggested change

	push!(jump_problems, JumpProblem(prob, alg, majumps, hopping_constants = hopping_constants, spatial_system = grids[1], save_positions = (false, false), rng = rng))
	push!(jump_problems, JumpProblem(prob, alg, majumps; hopping_constants, spatial_system = grids[1], save_positions = (false, false), rng))

[isaacsas]

Does this handle if there is no next jump?

[Vilin97]

To be honest, I do not rembmer any of the logic by now. Do we have a test for what you are asking? If not, we should, for all SSAs.

[Vilin97]

Addressed comments.

[Vilin97]

What do I need to do to get rid of the format and the spell check CI failures? @isaacsas

[isaacsas]

For spell check fix the listed spelling errors I guess. Haven’t ever used it or seen it before.

For formatting use JuliaFormatter to format the repo and add to the PR.

[isaacsas]

Just FYI, I'm traveling for the next week, but will try to review once I am back in Boston.

[isaacsas]

@Vilin97 sorry again for the long delays with this. I’ll try to put a test for the case I was worried about into the general tests in a followup. Assuming tests now pass I’ll merge.

I do think there is a chance the removal of the dispatch for HighLow could have performance implications, as I’m not sure Julia can always remove the branch in this situation (what it can do from the REPL may not be the same as what it can do here). But in any case we can merge and that can be examined later as a performance optimization.

Thanks again for your efforts on this — I look forward to updating the benchmarks to use it to see how it does.

[isaacsas]

Ahh, I didn’t realize that the updates to address my earlier comments actually never passed tests.