Merge pull request #4 from DanielVandH/extrapolation

Allow for extrapolation to be enabled

Author: DanielVandH

Project.toml

@@ -1,7 +1,7 @@
 name = "EpithelialDynamics1D"
 uuid = "ace8a2d7-7779-48a6-a8a4-cf6831a7e55b"
 authors = ["Daniel VandenHeuvel <danj.vandenheuvel@gmail.com>"]
-version = "1.7.2"
+version = "1.8.0"
 
 [deps]
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"

src/statistics.jl

@@ -127,7 +127,8 @@ end
         interp_fnc=(u, t) -> LinearInterpolation{true}(u, t),
         smooth_boundary=true)
 
-Computes summary statistics for the node densities from an `EnsembleSolution` to a [`CellProblem`](@ref).
+Computes summary statistics for the node densities from an `EnsembleSolution` to a [`CellProblem`](@ref). Any 
+negative values are set to zero.
 
 # Arguments 
 - `sol::EnsembleSolution`: The ensemble solution to a `CellProblem`.
@@ -140,6 +141,7 @@ Computes summary statistics for the node densities from an `EnsembleSolution` to
 - `alpha::Float64 = 0.05`: The significance level for the confidence intervals.
 - `interp_fnc = (u, t) -> LinearInterpolation{true}(u, t)`: The function to use for constructing the interpolant.
 - `smooth_boundary::Bool = true`: Whether to use the smooth boundary node densities.
+- `extrapolate = false`: Whether to allow for extrapolation when considering evaluating outside the range of cells for a given time and a given simulation.
 
 # Outputs 
 - `q::Vector{Vector{Vector{Float64}}}`: The node densities for each cell simulation.
@@ -156,7 +158,8 @@ function node_densities(sol::EnsembleSolution;
     knots=get_knots(sol, num_knots; indices, stat),
     alpha=0.05,
     interp_fnc=(u, t) -> LinearInterpolation{true}(u, t),
-    smooth_boundary=true)
+    smooth_boundary=true,
+    extrapolate=false)
     q = Vector{Vector{Vector{Float64}}}(undef, length(indices))
     r = Vector{Vector{Vector{Float64}}}(undef, length(indices))
     Base.Threads.@threads for i in eachindex(indices)
@@ -175,7 +178,7 @@ function node_densities(sol::EnsembleSolution;
             cell_positions = r[k][j]
             interp = interp_fnc(densities, cell_positions)
             for i in eachindex(knots[j])
-                if knots[j][i] > r[k][j][end]
+                if !extrapolate && knots[j][i] > r[k][j][end] 
                     q_splines[i, j, k] = 0.0
                 else
                     q_splines[i, j, k] = max(0.0, interp(knots[j][i]))

test/step_function.jl

@@ -613,4 +613,42 @@ end
             @test quantile(all_q, 0.975) ≈ uppers[j][i] rtol = 1e-3
         end
     end
+
+    # Using the average leading edge and allowing extrapolation
+    L = leading_edges(sol).L
+    _L = stack(L)
+    _indices = rand(eachindex(sol), 20)
+    _L = _L[:, _indices]
+    _mL = mean.(eachrow(_L))
+    q, r, means, lowers, uppers, knots = node_densities(sol; indices=_indices, stat=mean, extrapolate=true)
+    @inferred node_densities(sol; indices=_indices, stat=mean,extrapolate=true)
+    for j in eachindex(knots)
+        a = mean(sol[k][j][begin] for k in _indices)
+        b = mean(sol[k][j][end] for k in _indices)
+        @test knots[j] ≈ LinRange(a, b, 500)
+        @test knots[j][end] ≈ _mL[j]
+    end
+    for (enum_k, k) in enumerate(_indices)
+        for j in rand(1:length(sol[k]), 40)
+            for i in 1:length(sol[k][j])
+                if i == 1
+                    @test q[enum_k][j][1] ≈ LinearInterpolation([1 / (r[enum_k][j][2] - r[enum_k][j][1]), 2 / (r[enum_k][j][3] - r[enum_k][j][1])], [(r[enum_k][j][1] + r[enum_k][j][2]) / 2, r[enum_k][j][2]])(r[enum_k][j][1])
+                elseif i == length(sol[k][j])
+                    n = length(sol[k][j])
+                    @test q[enum_k][j][n] ≈ LinearInterpolation([2 / (r[enum_k][j][end] - r[enum_k][j][end-2]), 1 / (r[enum_k][j][end] - r[enum_k][j][end-1])], [r[enum_k][j][end-1], (r[enum_k][j][end-1] + r[enum_k][j][end]) / 2])(r[enum_k][j][end])
+                else
+                    @test q[enum_k][j][i] ≈ 2 / (r[enum_k][j][i+1] - r[enum_k][j][i-1])
+                end
+                @test r[enum_k][j][i] == sol[k][j][i]
+            end
+        end
+    end
+    for j in rand(eachindex(mb_sol), 40)
+        for i in eachindex(knots[j])
+            all_q = max.(0.0, [LinearInterpolation(q[k][j], r[k][j])(knots[j][i])  for k in eachindex(_indices)])
+            @test mean(all_q) ≈ means[j][i] rtol = 1e-3
+            @test quantile(all_q, 0.025) ≈ lowers[j][i] rtol = 1e-3
+            @test quantile(all_q, 0.975) ≈ uppers[j][i] rtol = 1e-3
+        end
+    end
 end