Avoid array copies in the calc_correction accumulation loop

[achubaty]

Problem

calc_correction builds null_current_total by adding the same array into overlapping regions, once for each cell of the block:

for i in 1:arguments["block_size"]
    for j in 1:arguments["block_size"]
        null_current_total[i:(...), j:(...)] += null_current
    end
end

Indexing with ranges on the left of += does not write in place. Each iteration copies the target region out, allocates the sum, and writes it back — two arrays of (2*radius + 1)² elements per iteration, block_size² times.

This runs before the moving window loop begins, on a single thread, so adding threads does not help.

Change

-            null_current_total[i:(i + arguments["radius"] * 2 + arguments["buffer"] * 2),
-                               j:(j + arguments["radius"] * 2 + arguments["buffer"] * 2)] += null_current
+            @views null_current_total[i:(i + arguments["radius"] * 2 + arguments["buffer"] * 2),
+                                      j:(j + arguments["radius"] * 2 + arguments["buffer"] * 2)] .+= null_current

The same additions happen in the same order; only the temporaries go. null_current_total and null_current are separate arrays, so there is no aliasing to worry about.

Results

Measured with the loop lifted out on its own, on Julia 1.11.7, buffer = 0. identical compares the raw bit patterns of the two result arrays, not ≈.

radius	block_size	identical	before	after	allocated before	after
50	11	yes	0.004 s	0.0004 s	19 MiB	0.09 MiB
150	15	yes	0.048 s	0.006 s	312 MiB	0.76 MiB
300	31	yes	1.95 s	0.100 s	5,300 MiB	3.04 MiB
503	51	yes	10.9 s	1.36 s	39.3 GiB	8.52 MiB

The allocation figures are exact and repeatable. The times were taken on a shared machine that was busy at the time, so read them as roughly an order of magnitude rather than as precise ratios.

At the sizes used in the test suite (radius = 5) the difference is not measurable either way.

Test suite on Julia 1.11.7 with Circuitscape 5.15.0: Internals 10/10, run_omniscape() 31/31.

Script that produces the table

Standalone — no Omniscape state, nothing outside the standard library.

using Printf

## as written today
function current(null_current, bs, r)
    tot = zeros(Float64, 2r + bs, 2r + bs)
    for i in 1:bs, j in 1:bs
        tot[i:(i + 2r), j:(j + 2r)] += null_current
    end
    tot
end

## proposed
function withviews(null_current, bs, r)
    tot = zeros(Float64, 2r + bs, 2r + bs)
    for i in 1:bs, j in 1:bs
        @views tot[i:(i + 2r), j:(j + 2r)] .+= null_current
    end
    tot
end

bits(a) = reinterpret(UInt64, vec(a))
best(f, args...) = minimum(@elapsed(f(args...)) for _ in 1:5)

@printf("%-8s %-11s %-10s %10s %10s %8s %12s %10s\n",
        "radius", "block_size", "identical", "current", "@views", "faster", "alloc before", "after")
for (r, bs) in ((5, 3), (50, 11), (150, 15), (300, 31), (503, 51))
    nc = rand(Float64, 2r + 1, 2r + 1)
    a, b = current(nc, bs, r), withviews(nc, bs, r)
    ta, tb = best(current, nc, bs, r), best(withviews, nc, bs, r)
    aa, ab = @allocated(current(nc, bs, r)), @allocated(withviews(nc, bs, r))
    @printf("%-8d %-11d %-10s %9.4fs %9.4fs %7.1fx %10.1f MiB %6.2f MiB\n",
            r, bs, bits(a) == bits(b), ta, tb, ta / tb, aa / 2^20, ab / 2^20)
end

A larger change, if you want it

The loop computes a two-dimensional convolution of null_current with a block_size × block_size box. A box filter is separable, so two passes of a sliding-window sum give the same result in far less work — about 650 times faster than the current code at radius = 503, against about 8 times for the change proposed here.

I have not included it. It changes the order of summation, so results differ by around 1e-13 relative, and the stored reference rasters would need to be looked at. Happy to open it separately if that is of interest.

🤖 Generated with Claude Code

[ViralBShah]

@achubaty I'll be happy to give you commit access here if you are comfortable with the codebase. While I am familiar with Circuitscape.jl, I am not familiar with Omniscape.jl - and it is great to see you opening PRs.

Let me know if you are comfortable with commit here.

[vlandau]

@ViralBShah would commit access to non-main branches as a start make things a bit easier, or is there not much point in that? I think we still want reviews and to keep main protected. I know I'm saying this as someone who hasn't touched the project in a long while, so I don't want to get in the way either! I'm happy to provide some general support on PR reviews, though I do want to make a disclaimer that I'm quite out of practice with Julia these days.

[ViralBShah]

@vlandau That makes sense. I was hoping to see if someone could claude their way into updating the package to 1.12 and fix some of the tolerance issues cropping up.

[vlandau]

@ViralBShah yeah that may be the likeliest path to an update (finally). Will just want to make sure the code is reviewed when AI-contributed. The very nice thing here though is that we have canonical correctness tests, so we can feel confident on changes if tests pass. Guards us a bit against AI errors (but not necessarily overly verbose and over-engineered AI slop that you can see when larger tasks are dispatched to a model like opus 5 (in my experience)). I will think about this. Perhaps I will be able to do a short hackathon to update the package over the weekend some time.

[achubaty]

@ViralBShah Thank you for the offer of commit access. I agree with @vlandau that non-main level access is preferred, so I can make small changes and have them reviewed as I'm still getting my head around things here (and using Claude).

I will need to update this PR to rebase using the updated julia@v3 action. THis is forthcoming

[achubaty]

The red checks here are the same infrastructure problem I described on #168, not anything to do with this change. The run died in dependency resolution, before a single test executed.

From the ubuntu job's log:

LinearSolve has a malformed Project.toml, the extension package SparseArrays is not listed in [weakdeps]
##[error]Process completed with exit code 1.

Briefly: this run used julia-actions/setup-julia@v2, where version: min against julia = "~1.11" resolves to Julia 1.11.0. That is the only patch in the entire 1.11 series whose bundled Pkg requires an extension trigger to appear in [weakdeps]; from 1.11.1 onward a trigger in [deps] is accepted. LinearSolve v3.87.0 declares SparseArrays in [deps] and uses it to trigger 15 extensions, so resolution fails outright. #164 bumped the action to @v3 — where min means the latest patch of the minimum minor, i.e. 1.11.9 — but it merged after this run had already started, so the run was computed against the older workflow.

I have rebased onto current main to pick that up. A plain "Re-run all jobs" would not help; it replays the recorded merge commit, still on @v2.

One thing worth flagging: unlike #168, this rebase will not turn CI green on its own. This branch carries no Circuitscape pin, so it resolves 5.17.1 and will now get past resolution only to hit the three failures main currently has (runtests.jl:102, :103, and the CHOLMOD residual at :95), plus the docs example. Those are the subject of #168. So the expected outcome here is red-for-main's-reason rather than red-for-a-stale-workflow — which at least means this change finally gets exercised, since calc_correction sits directly on the path those tests exercise. It should go green once #168 lands, and I am happy to rebase again then.