update

src/allee.jl

@@ -23,3 +23,10 @@ AlleeExtinction{R,W}(; minfounders=Param(5.0, bounds=(1.0, 200.0))) where {R,W}
 
     return N >= f ? N : zero(N)
 end
+@inline function applyrule(data, rule::AlleeExtinction, Ns::AbstractArray, I)
+    fs = get(data, rule.minfounders, I...)
+
+    return broadcast(Ns, fs) do N, f
+        N >= f ? N : zero(N)
+    end
+end

src/growth.jl

@@ -53,7 +53,7 @@ modifyrule(rule::ExponentialGrowth, data) = precalc_nsteps(rule, data)
     N > zero(N) || return zero(N)
     rt = get(data, rule.rate, I...) * rule.nsteps
 
-    return @fastmath N * exp(rt)
+    return N * exp(rt)
 end
 
 """
@@ -106,12 +106,32 @@ modifyrule(rule::LogisticGrowth, data) = precalc_nsteps(rule, data)
     rt = get(data, rule.rate, I...) * rule.nsteps
     k = get(data, rule.carrycap, I...)
 
-    if rt > zero(rt)
-        return @fastmath (N * k) / (N + (k - N) * exp(-rt))
+    new_N = if rt > zero(rt)
+        (N * k) / (N + (k - N) * exp(-rt))
     else
-        return @fastmath N * exp(rt)
+        N * exp(rt)
     end
+    isnan(new_NS) && _nan_pop_error(new_N, Ns, k)
+    return min(max(zero(new_N), new_N), k)
 end
+@inline function applyrule(data, rule::LogisticGrowth, Ns::AbstractArray, I)
+    rts = get(data, rule.rate, I...) .* rule.nsteps
+    ks = get(data, rule.carrycap, I...)
+
+    new_Ns = map(Ns, rts, ks) do N, rt, k
+        if rt > zero(rt)
+            (N * k) / (N + (k - N) * exp(-rt))
+        else
+            N * exp(rt)
+        end
+    end
+    any(isnan, new_Ns) && return zero(Ns)#_nan_pop_error(new_Ns, Ns, ks)
+    any(isinf, new_Ns) && return ks #_inf_pop_error(new_Ns, Ns, ks)
+    return min.(max.(zero(eltype(new_Ns)), new_Ns), ks)
+end
+
+@noinline _nan_pop_error(new_N, N, k) = error("NaN population found: $new_N, from original $N and carrycap $k")
+@noinline _inf_pop_error(new_N, N, k) = error("Inf population found: $new_N, from original $N and carrycap $k")
 
 """
     ThresholdGrowth <: CellRule

src/kernel/distancemethods.jl

@@ -33,7 +33,8 @@ Joseph D. Chipperfield et al 2011"
 """
 struct CentroidToCentroid <: DistanceMethod end
 
-dispersalprob(f, ::CentroidToCentroid, x, y, cellsize) = f(sqrt(x^2 + y^2) * cellsize)
+dispersalprob(f, ::CentroidToCentroid, x, y, cellsize::T) where T =
+    f(T(sqrt(x^2 + y^2)) * cellsize)
 
 """
     AreaToCentroid <: DistanceMethod
@@ -52,15 +53,15 @@ Base.@kwdef struct AreaToCentroid{SS<:Number} <: DistanceMethod
     subsampling::SS = Param(10.0; bounds=(2.0, 40.0))
 end
 
-@inline function dispersalprob(f, dm::AreaToCentroid, x, y, cellsize)
-    prob = zero(cellsize)
-    centerfirst = 1 / subsampling(dm) / 2 - 0.5
+@inline function dispersalprob(f, dm::AreaToCentroid, x, y, cellsize::T) where T
+    prob = zero(T)
+    centerfirst = 1 / subsampling(dm) / 2 - oneunit(T) / 2
     centerlast = centerfirst * -1
     range = LinRange(centerfirst, centerlast, round(Int, subsampling(dm)))
     for j in range, i in range
         prob += sqrt((x + i)^2 + (y + j)^2) * cellsize |> f
     end
-    prob / subsampling(dm)^2
+    return T(prob / subsampling(dm)^2)
 end
 
 """
@@ -79,16 +80,16 @@ Base.@kwdef struct AreaToArea{SS<:Number} <: DistanceMethod
     subsampling::SS = Param(10.0; bounds=(2.0, 40.0))
 end
 
-@inline function dispersalprob(f, dm::AreaToArea, x, y, cellsize)
-    prob = zero(cellsize)
+@inline function dispersalprob(f, dm::AreaToArea, x, y, cellsize::T) where T
+    prob = zero(T)
     # Get the center point of the first cell (for both dimensions)
-    centerfirst = 1 / subsampling(dm) / 2 - 0.5
+    centerfirst = 1 / subsampling(dm) / 2 - oneunit(T) / 2
     centerlast = centerfirst * -1
     range = LinRange(centerfirst, centerlast, round(Int, subsampling(dm)))
     for i in range, j in range
         for a in range, b in range
-            prob += sqrt((x + i + a)^2 + (y + j + b)^2) * cellsize |> f
+            prob += T(sqrt((x + i + a)^2 + (y + j + b)^2) * cellsize) |> f
         end
     end
-    prob / subsampling(dm)^4
+    return T(prob / subsampling(dm)^4)
 end

src/kernel/kernel.jl

@@ -64,6 +64,7 @@ function DG.Stencils.rebuild(n::DispersalKernel{R,N,L,<:Any,<:Any,K,F,C,D}, buff
         newstencil, kernel(n), formulation(n), cellsize(n), distancemethod(n)
     )
 end
+DG.Stencils.neighbors(hood::DispersalKernel) = neighbors(hood.stencil)
 
 cellsize(stencil::DispersalKernel) = stencil.cellsize
 distancemethod(stencil::DispersalKernel) = stencil.distancemethod