Merge pull request #55 from cesaraustralia/perf

Perf

Project.toml

@@ -1,7 +1,7 @@
 name = "DynamicGrids"
 uuid = "a5dba43e-3abc-5203-bfc5-584ca68d3f5b"
 authors = ["Rafael Schouten <[email protected]>"]
-version = "0.10.4"
+version = "0.10.5"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"

src/extent.jl

@@ -16,6 +16,17 @@ mutable struct Extent{I,M,A}
     aux::A
     tspan::AbstractRange
 end
+Extent(init::I, mask::M, aux::A, tspan::T) where {I,M,A,T} = begin
+    # Check grid sizes match
+    gridsize = if init isa NamedTuple
+        size_ = size(first(init_))
+        all(map(i -> size(i) == size_, init)) || throw(ArgumentError("`init` grid sizes do not match"))
+    else
+        size_ = size(init_)
+    end
+    (mask !== nothing) && (size(mask) != size_) && throw(ArgumentError("`mask` size do not match `init`"))
+    Extent{I,M,A,T}(init, mask, aux, tspan)
+end
 Extent(; init, mask=nothing, aux=nothing, tspan, kwargs...) =
     Extent(init, mask, aux, tspan)
 

src/framework.jl

@@ -70,6 +70,7 @@ end
 
 Run a simulation passing in rules without defining a `Ruleset`.
 """
+sim!(output::Output, rules::Tuple; kwargs...) = sim!(output::Output, rules...; kwargs...)
 sim!(output::Output, rules::Rule...; tspan=tspan(output), kwargs...) = begin
     ruleset = Ruleset(rules...; timestep=step(tspan), kwargs...)
     sim!(output::Output, ruleset; tspan=tspan, kwargs...)
@@ -105,7 +106,6 @@ function resume!(output::GraphicOutput, ruleset::Ruleset=ruleset(output);
         nreplicates=nothing)
     initialise(output)
     # Check status and arguments
-    length(output) > 0 || error("There is no simulation to resume. Run `sim!` first")
     isrunning(output) && error("A simulation is already running in this output")
     setrunning!(output, true) || error("Could not start the simulation with this output")
 

src/neighborhoods.jl

@@ -1,7 +1,7 @@
 
 """
-Neighborhoods define the pattern of surrounding cells in the "neighborhood" 
-of the current cell. The `neighbors` function returns the surrounding 
+Neighborhoods define the pattern of surrounding cells in the "neighborhood"
+of the current cell. The `neighbors` function returns the surrounding
 cells as an iterable.
 
 The main kinds of neighborhood are demonstrated below:
@@ -53,7 +53,7 @@ end
 """
     Moore(radius::Int=1)
 
-Moore neighborhoods define the neighborhood as all cells within a horizontal or 
+Moore neighborhoods define the neighborhood as all cells within a horizontal or
 vertical distance of the central cell. The central cell is omitted.
 
 The `buffer` argument may be required for performance
@@ -65,9 +65,9 @@ end
 # Buffer is updated later during the simulation.
 # but can be passed in now to avoid the allocation.
 # This might be bad design. SimData could instead hold a list of
-# ruledata for the rule that holds this buffer, with 
+# ruledata for the rule that holds this buffer, with
 # the neighborhood. So you can do neighbors(data)
-Moore(radius::Int=1, buffer=nothing) = 
+Moore(radius::Int=1, buffer=nothing) =
     Moore{radius}(buffer)
 Moore{R}(buffer=nothing) where R =
     Moore{R,typeof(buffer)}(buffer)
@@ -99,8 +99,8 @@ end
 
 """
 Neighborhoods are tuples or vectors of custom coordinates tuples
-that are specified in relation to the central point of the current cell. 
-They can be any arbitrary shape or size, but should be listed in column-major 
+that are specified in relation to the central point of the current cell.
+They can be any arbitrary shape or size, but should be listed in column-major
 order for performance.
 """
 abstract type AbstractPositional{R,B} <: Neighborhood{R,B} end
@@ -113,8 +113,8 @@ const CustomCoords = Union{AbstractArray{<:CustomCoord},Tuple{Vararg{<:CustomCoo
     Positional(coords::Tuple{Tuple{Vararg{Int}}}, [buffer=nothing])
     Positional{R}(coords::Tuple, buffer)
 
-Neighborhoods that can take arbitrary shapes by specifying each coordinate, 
-as `Tuple{Int,Int}` of the row/column distance (positive and negative) 
+Neighborhoods that can take arbitrary shapes by specifying each coordinate,
+as `Tuple{Int,Int}` of the row/column distance (positive and negative)
 from the central point.
 
 The neighborhood radius is calculated from the most distance coordinate.
@@ -192,7 +192,7 @@ end
 """
     neighbors(hood::Positional)
 
-Returns a tuple of iterators over each `Positional` neighborhood 
+Returns a tuple of iterators over each `Positional` neighborhood
 layer, for the cells around the current index.
 """
 @inline neighbors(hood::LayeredPositional) =
@@ -206,7 +206,7 @@ layer, for the cells around the current index.
 """
     VonNeumann(radius=1)
 
-A convenience wrapper to build Von-Neumann neighborhoods as 
+A convenience wrapper to build Von-Neumann neighborhoods as
 a [`Positional`](@ref) neighborhood.
 """
 VonNeumann(radius=1, buffer=nothing) = begin
@@ -240,9 +240,9 @@ radius(rule::Rule, args...) = 0
 
 # Build rules and neighborhoods for each buffer, so they
 # don't have to be constructed in the loop
-spreadbuffers(chain::Chain, init::AbstractArray) = begin
+spreadbuffers(chain::Chain{R,W}, init::AbstractArray) where {R,W} = begin
     buffers, bufrules = spreadbuffers(rules(chain)[1], init)
-    buffers, map(r -> Chain(r, tail(rules(chain))...), bufrules)
+    buffers, map(r -> Chain{R,W}((r, tail(rules(chain))...)), bufrules)
 end
 spreadbuffers(rule::Rule, init::AbstractArray) =
     spreadbuffers(rule, neighborhood(rule), buffer(neighborhood(rule)), init)
@@ -255,14 +255,14 @@ spreadbuffers(rule::NeighborhoodRule, hood::Neighborhood, buffers::Tuple, init::
     allocbuffers(init::AbstractArray, hood::Neighborhood)
     allocbuffers(init::AbstractArray, radius::Int)
 
-Allocate buffers for the Neighborhood. The `init` array should 
+Allocate buffers for the Neighborhood. The `init` array should
 be of the same type as the grid the neighborhood runs on.
 """
-allocbuffers(init::AbstractArray, hood::Neighborhood) = allocbuffers(init, radius(hood))
-allocbuffers(init::AbstractArray, r::Int) = Tuple(allocbuffer(init, r) for i in 1:2r)
+@inline allocbuffers(init::AbstractArray, hood::Neighborhood{R}) where R = allocbuffers(init, R)
+@inline allocbuffers(init::AbstractArray, r::Int) = ntuple(i -> allocbuffer(init, r), 2r)
 
-allocbuffer(init::AbstractArray, hood::Neighborhood) = allocbuffer(init, radius(hood))
-allocbuffer(init::AbstractArray, r::Int) = zeros(eltype(init), 2r+1, 2r+1)
+@inline allocbuffer(init::AbstractArray, hood::Neighborhood{R}) where R = allocbuffer(init, R)
+@inline allocbuffer(init::AbstractArray, r::Int) = zeros(eltype(init), 2r+1, 2r+1)
 
 """
     hoodsize(radius)

src/outputs/graphic.jl

@@ -72,7 +72,7 @@ isshowable(o::GraphicOutput, f) = true # TODO working max fps. o.timestamp + (t
 
 storeframe!(o::GraphicOutput, data::AbstractSimData) = begin
     f = frameindex(o, data)
-    if isstored(o)
+    if f > length(o)
         _pushgrid!(eltype(o), o)
     end
     storeframe!(eltype(o), o, data, f)

src/outputs/output.jl

@@ -29,9 +29,9 @@ Base.@propagate_inbounds Base.setindex!(o::Output, x, i::Union{Int,AbstractVecto
 Base.push!(o::Output, x) = push!(frames(o), x)
 Base.step(o::Output) = step(tspan(o))
 
-DimensionalData.DimensionalArray(o::Output{<:NamedTuple}; key=first(keys(o[1]))) = 
+DimensionalData.DimensionalArray(o::Output{<:NamedTuple}; key=first(keys(o[1]))) =
     cat(map(f -> f[key], frames(o)...); dims=timedim(o))
-DimensionalData.DimensionalArray(o::Output{<:DimensionalArray}) = 
+DimensionalData.DimensionalArray(o::Output{<:DimensionalArray}) =
     cat(frames(o)...; dims=timedim(o))
 DimensionalData.dims(o::Output) = begin
     ts = tspan(o)
@@ -141,7 +141,7 @@ storeframe!(::Type{<:AbstractArray}, output::Output, simdata::AbstractSimData, f
     _storeloop(outgrid, first(grids(simdata)))
 end
 _storeloop(outgrid, grid) =
-    for j in axes(outgrid, 2), i in axes(outgrid, 1) 
+    for j in axes(outgrid, 2), i in axes(outgrid, 1)
         @inbounds outgrid[i, j] = grid[i, j]
     end
 # Replicated frames

src/overflow.jl

@@ -97,12 +97,14 @@ handleoverflow!(griddata::GridData, ::WrapOverflow) = begin
 end
 
 function wrapstatus!(status)
-    status[end, :] .|= status[1, :]
-    status[:, end] .|= status[:, 1]
-    status[end-1, :] .|= status[1, :] .|= status[2, :]
-    status[:, end-1] .|= status[:, 1] .|= status[:, 2]
-    status[1, :] .|= status[end-1, :] .|= status[end, :]
-    status[:, 1] .|= status[:, end-1] .|= status[:, end]
-    status[end-1, :] .= true
+    # This could be further optimised.
+    status[end-1, :] .|= status[1, :]
+    status[:, end-1] .|= status[:, 1]
+    #status[end-2, :] .|= status[1, :] .|= status[2, :]
+    status[end-2, :] .= true
+    status[:, end-2] .|= status[:, 1] .|= status[:, 2]
+    #status[1, :] .|= status[end-2, :] .|= status[end-1, :]
     status[1, :] .= true
+    status[:, 1] .|= status[:, end-2] .|= status[:, end-1]
+    status .= true
 end

src/precalc.jl

@@ -10,8 +10,8 @@ end
     precalcrules(rule::Rule, simdata::SimData)
 
 Precalculates rule at each timestep, if there are any fields that need
-to be updated over time. Rules are usually immutable (it's faster), so 
-return a whole new rule object with changes you need applied. 
+to be updated over time. Rules are usually immutable (it's faster), so
+return a whole new rule object with changes you need applied.
 They will be discarded, and `rule` will always be the original object passed in.
 
 Setfield.jl and Flatten.jl may help for this.
@@ -23,7 +23,5 @@ precalcrules(rule, simdata) = rule
 precalcrules(rules::Tuple, simdata) =
     (precalcrules(rules[1], simdata), precalcrules(tail(rules), simdata)...)
 precalcrules(rules::Tuple{}, simdata) = ()
-precalcrules(chain::Chain{R,W}, simdata) where {R,W} = begin
-    ch = precalcrules(rules(chain), simdata)
-    Chain{R,W,typeof(ch)}(ch)
-end
+precalcrules(chain::Chain{R,W}, simdata) where {R,W} =
+    Chain{R,W}(precalcrules(rules(chain), simdata))

src/rules.jl

@@ -194,9 +194,9 @@ Cell(f; read=:_default_, write=read) = Cell{read,write}(f)
 
 const Map = Cell
 
-astuple(rule::Rule, read) = astuple(readkeys(rule), read)
-astuple(::Tuple, read) = read
-astuple(::Symbol, read) = (read,)
+astuple(rule::Rule, state) = astuple(readkeys(rule), state)
+astuple(::Tuple, state) = state
+astuple(::Symbol, state) = (state,)
 
 """
     Neighbors(f, neighborhood)

src/rulesets.jl

@@ -57,11 +57,12 @@ struct NoOpt <: PerformanceOpt end
 abstract type AbstractRuleset end
 
 # Getters
+ruleset(rs::AbstractRuleset) = rs
+rules(rs::AbstractRuleset) = rs.rules
 overflow(rs::AbstractRuleset) = rs.overflow
 opt(rs::AbstractRuleset) = rs.opt
 cellsize(rs::AbstractRuleset) = rs.cellsize
 timestep(rs::AbstractRuleset) = rs.timestep
-ruleset(rs::AbstractRuleset) = rs
 
 Base.step(rs::AbstractRuleset) = timestep(rs)
 
@@ -83,7 +84,7 @@ Rules will be run in the order they are passed, ie. `Ruleset(rule1, rule2, rule3
 """
 @default_kw @flattenable mutable struct Ruleset{O<:Overflow,Op<:PerformanceOpt,C,T} <: AbstractRuleset
     # Rules are intentionally not type stable. This allows `precalc` and Interact.jl 
-    # updates to change the rule type. Function barriers remove any performance overheads.
+    # updates to change the rule type. Function barriers remove most performance overheads.
     rules::Tuple{Vararg{<:Rule}} | ()               | true
     overflow::O                  | RemoveOverflow() | false
     opt::Op                      | NoOpt()          | false
@@ -92,5 +93,4 @@ Rules will be run in the order they are passed, ie. `Ruleset(rule1, rule2, rule3
 end
 Ruleset(rules::Vararg{<:Rule}; kwargs...) = Ruleset(; rules=rules, kwargs...)
 Ruleset(rules::Tuple; kwargs...) = Ruleset(; rules=rules, kwargs...)
-
-rules(rs::Ruleset) = rs.rules
+Ruleset(rs::Ruleset) = Ruleset(rules(rs), overflow(rs), opt(rs), cellsize(rs), timestep(rs))

src/simulationdata.jl

@@ -45,7 +45,9 @@ localstatus(d::GridData) = d.localstatus
 gridsize(d::GridData) = size(init(d))
 gridsize(A::AbstractArray) = size(A)
 gridsize(nt::NamedTuple) = gridsize(first(nt))
+gridsize(nt::NamedTuple{(),Tuple{}}) = 0, 0
 gridsize(t::Tuple) = gridsize(first(t))
+gridsize(t::Tuple{}) = 0, 0
 
 
 """
@@ -65,31 +67,32 @@ ReadableGridData(init::AbstractArray, mask, radius, overflow) = begin
         hoodsize = 2r + 1
         blocksize = 2r
         source = addpadding(init, r)
-        nblocs = indtoblock.(size(source), blocksize)
+        dest = addpadding(init, r)
+        nblocs = indtoblock.(size(source), blocksize) .+ 1
         sourcestatus = zeros(Bool, nblocs)
-        deststatus = deepcopy(sourcestatus)
+        deststatus = zeros(Bool, nblocs)
         updatestatus!(source, sourcestatus, deststatus, r)
         localstatus = zeros(Bool, 2, 2)
     else
         source = deepcopy(init)
+        dest = deepcopy(init)
         sourcestatus = deststatus = true
         localstatus = nothing
     end
-    dest = deepcopy(source)
 
     ReadableGridData(init, mask, radius, overflow, source, dest,
                      sourcestatus, deststatus, localstatus)
 end
 
 Base.parent(d::ReadableGridData) = parent(source(d))
-# Base.@propagate_inbounds
-Base.getindex(d::ReadableGridData, I...) = getindex(source(d), I...)
+Base.@propagate_inbounds Base.getindex(d::ReadableGridData, I...) = getindex(source(d), I...)
+
 
 """
     ReadableGridData(griddata::GridData)
 
 Passed to rules `<: ManualRule`, and can be written to directly as
-an array. This handles updates to SparseOpt() and writing to 
+an array. This handles updates to SparseOpt() and writing to
 the correct source/dest array.
 """
 @GridDataMixin struct WritableGridData{} <: GridData{T,N,I} end
@@ -103,7 +106,7 @@ Base.@propagate_inbounds Base.setindex!(d::WritableGridData, x, I...) = begin
 end
 
 Base.parent(d::WritableGridData) = parent(dest(d))
-Base.@propagate_inbounds Base.getindex(d::WritableGridData, I...) = 
+Base.@propagate_inbounds Base.getindex(d::WritableGridData, I...) =
     getindex(dest(d), I...)
 
 
@@ -127,19 +130,20 @@ struct SimData{G<:NamedTuple,E,Ru,F} <: AbstractSimData
     ruleset::Ru
     currentframe::F
 end
-SimData(extent, ruleset::Ruleset) = begin
-    extent = asnamedtuple(extent)
+# Convert grids in extent to NamedTuple
+SimData(extent, ruleset::Ruleset) =
+    SimData(asnamedtuple(extent), ruleset::Ruleset)
+SimData(extent::Extent{<:NamedTuple{Keys}}, ruleset::Ruleset) where Keys = begin
     # Calculate the neighborhood radus (and grid padding) for each grid
-    keys_ = keys(init(extent))
-    radii = NamedTuple{keys_}(get(radius(ruleset), key, 0) for key in keys_)
+    radii = NamedTuple{Keys}(get(radius(ruleset), key, 0) for key in Keys)
     # Construct the SimData for each grid
     griddata = map(init(extent), radii) do in, ra
         ReadableGridData(in, mask(extent), ra, overflow(ruleset))
     end
     SimData(griddata, extent, ruleset)
 end
 SimData(griddata::NamedTuple, extent, ruleset::Ruleset) = begin
-    currentframe = 1; 
+    currentframe = 1;
     SimData(griddata, extent, ruleset, currentframe)
 end
 
@@ -158,14 +162,14 @@ currenttime(d::SimData) = tspan(d)[currentframe(d)]
 currenttime(d::Vector{<:SimData}) = currenttime(d[1])
 
 # Getters forwarded to data
-Base.getindex(d::SimData, i::Symbol) = 
+Base.getindex(d::SimData, i::Symbol) =
     getindex(grids(d), i)
 # For resolving method ambiguity
-Base.getindex(d::SimData{<:NamedTuple{(:_default_,)}}, i::Symbol) = 
+Base.getindex(d::SimData{<:NamedTuple{(:_default_,)}}, i::Symbol) =
     getindex(grids(d), i)
-Base.getindex(d::SimData{<:NamedTuple{(:_default_,)}}, I...) = 
+Base.getindex(d::SimData{<:NamedTuple{(:_default_,)}}, I...) =
     getindex(first(grids(d)), I...)
-Base.setindex!(d::SimData{<:NamedTuple{(:_default_,)}}, x, I...) = 
+Base.setindex!(d::SimData{<:NamedTuple{(:_default_,)}}, x, I...) =
     setindex!(first(grids(d)), x, I...)
 Base.keys(d::SimData) = keys(grids(d))
 Base.values(d::SimData) = values(grids(d))