-
Notifications
You must be signed in to change notification settings - Fork 7
/
simulationdata.jl
271 lines (234 loc) · 9.06 KB
/
simulationdata.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
"""
Simulation data specific to a single grid.
"""
abstract type GridData{T,N,I} <: AbstractArray{T,N} end
(::Type{T})(d::GridData) where T <: GridData =
T(init(d), mask(d), radius(d), overflow(d), source(d), dest(d),
sourcestatus(d), deststatus(d), localstatus(d))
# Common fields for GridData and WritableGridData, which are
# identical except for their indexing methods
@mix struct GridDataMixin{T,N,I<:AbstractArray{T,N},M,R,O,S,St,LSt}
init::I
mask::M
radius::R
overflow::O
source::S
dest::S
sourcestatus::St
deststatus::St
localstatus::LSt
end
GridDataOrReps = Union{GridData, Vector{<:GridData}}
# Array interface
Base.size(d::GridData) = size(source(d))
Base.axes(d::GridData) = axes(source(d))
Base.eltype(d::GridData) = eltype(source(d))
Base.firstindex(d::GridData) = firstindex(source(d))
Base.lastindex(d::GridData) = lastindex(source(d))
# Getters
init(d::GridData) = d.init
mask(d::GridData) = d.mask
radius(d::GridData) = d.radius
radius(d::Tuple{<:GridData,Vararg}) = map(radius, d)
overflow(d::GridData) = d.overflow
source(d::GridData) = d.source
dest(d::GridData) = d.dest
sourcestatus(d::GridData) = d.sourcestatus
deststatus(d::GridData) = d.deststatus
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
"""
ReadableGridData(griddata::GridData)
ReadableGridData(init::AbstractArray, mask, radius, overflow)
Simulation data and storage passed to rules for each timestep.
"""
@GridDataMixin struct ReadableGridData{} <: GridData{T,N,I} end
# Generate simulation data to match a ruleset and init array.
ReadableGridData(init::AbstractArray, mask, radius, overflow) = begin
r = radius
# We add one extra row and column of status blocks so
# we dont have to worry about special casing the last block
if r > 0
hoodsize = 2r + 1
blocksize = 2r
source = addpadding(init, r)
dest = addpadding(init, r)
nblocs = indtoblock.(size(source), blocksize) .+ 1
sourcestatus = zeros(Bool, nblocs)
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
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...)
"""
ReadableGridData(griddata::GridData)
Passed to rules `<: ManualRule`, and can be written to directly as
an array. This handles updates to SparseOpt() and writing to
the correct source/dest array.
"""
@GridDataMixin struct WritableGridData{} <: GridData{T,N,I} end
Base.@propagate_inbounds Base.setindex!(d::WritableGridData, x, I...) = begin
r = radius(d)
@inbounds dest(d)[I...] = x
if deststatus(d) isa AbstractArray
@inbounds deststatus(d)[indtoblock.(I .+ r, 2r)...] = true
end
end
Base.parent(d::WritableGridData) = parent(dest(d))
Base.@propagate_inbounds Base.getindex(d::WritableGridData, I...) =
getindex(dest(d), I...)
abstract type AbstractSimData end
"""
SimData(extent::Extent, ruleset::Ruleset)
Simulation dataset to hold all intermediate arrays, timesteps
and frame numbers for the current frame of the simulation.
A simdata object is accessable in [`applyrule`](@ref) as the first parameter.
Multiple grids can be indexed into using their key. In single grid
simulations `SimData` can be indexed directly as if it is a `Matrix`.
"""
struct SimData{G<:NamedTuple,E,Ru,F} <: AbstractSimData
grids::G
extent::E
ruleset::Ru
currentframe::F
end
# 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
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;
SimData(griddata, extent, ruleset, currentframe)
end
# Getters
extent(d::SimData) = d.extent
ruleset(d::SimData) = d.ruleset
grids(d::SimData) = d.grids
init(d::SimData) = init(extent(d))
mask(d::SimData) = mask(first(d))
aux(d::SimData) = aux(extent(d))
tspan(d::SimData) = tspan(extent(d))
timestep(d::SimData) = step(tspan(d))
currentframe(d::SimData) = d.currentframe
currenttime(d::SimData) = tspan(d)[currentframe(d)]
currenttime(d::Vector{<:SimData}) = currenttime(d[1])
# Getters forwarded to data
Base.getindex(d::SimData, i::Symbol) =
getindex(grids(d), i)
# For resolving method ambiguity
Base.getindex(d::SimData{<:NamedTuple{(:_default_,)}}, i::Symbol) =
getindex(grids(d), i)
Base.getindex(d::SimData{<:NamedTuple{(:_default_,)}}, I...) =
getindex(first(grids(d)), 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))
Base.first(d::SimData) = first(grids(d))
Base.last(d::SimData) = last(grids(d))
gridsize(d::SimData) = gridsize(first(d))
rules(d::SimData) = rules(ruleset(d))
overflow(d::SimData) = overflow(ruleset(d))
opt(d::SimData) = opt(ruleset(d))
# Get the actual current timestep, e.g. seconds instead of variable periods like Month
currenttimestep(d::SimData) = currenttime(d) + timestep(d) - currenttime(d)
# Swap source and dest arrays. Allways returns regular SimData.
swapsource(d::Tuple) = map(swapsource, d)
swapsource(grid::GridData) = begin
src = grid.source
dst = grid.dest
@set! grid.dest = src
@set! grid.source = dst
srcstatus = grid.sourcestatus
dststatus = grid.deststatus
@set! grid.deststatus = srcstatus
@set grid.sourcestatus = dststatus
end
# Uptate timestamp
updatetime(simdata::SimData, f::Integer) = begin
@set! simdata.currentframe = f
end
updatetime(simdata::AbstractVector{<:SimData}, f) = updatetime.(simdata, f)
#=
Find the maximum radius required by all rules
Add padding around the original init array, offset into the negative
So that the first real cell is still 1, 1
=#
addpadding(init::AbstractArray{T,N}, r) where {T,N} = begin
sze = size(init)
paddedsize = sze .+ 2r
paddedindices = -r + 1:sze[1] + r, -r + 1:sze[2] + r
sourceparent = fill!(similar(init, paddedsize), zero(T))
source = OffsetArray(sourceparent, paddedindices...)
# Copy the init array to the middle section of the source array
for j in 1:size(init, 2), i in 1:size(init, 1)
@inbounds source[i, j] = init[i, j]
end
source
end
#=
Initialise the block status array.
This tracks whether anything has to be done in an area of the main array.
=#
updatestatus!(grid::GridData) =
updatestatus!(parent(source(grid)), sourcestatus(grid), deststatus(grid), radius(grid))
updatestatus!(source, sourcestatus::Bool, deststatus::Bool, radius) = nothing
updatestatus!(source, sourcestatus, deststatus, radius) = begin
blocksize = 2radius
source = parent(source)
for i in CartesianIndices(source)
# Mark the status block if there is a non-zero value
if source[i] != 0
bi = indtoblock.(Tuple(i), blocksize)
@inbounds sourcestatus[bi...] = true
@inbounds deststatus[bi...] = true
end
end
end
# When replicates are an Integer, construct a vector of SimData
initdata!(::Nothing, extent, ruleset::Ruleset, nreplicates::Integer) =
[SimData(extent, ruleset) for r in 1:nreplicates]
# When simdata is a Vector, the existing SimData arrays are re-initialised
initdata!(simdata::AbstractVector{<:AbstractSimData}, extent, ruleset, nreplicates::Integer) =
map(d -> initdata!(d, extent, ruleset, nothing), simdata)
# When no simdata is passed in, create new SimData
initdata!(::Nothing, extent, ruleset::Ruleset, nreplicates::Nothing) =
SimData(extent, ruleset)
# Initialise a SimData object with a new `Extent` and `Ruleset`.
initdata!(simdata::AbstractSimData, extent::Extent, ruleset::Ruleset, nreplicates::Nothing) = begin
map(values(simdata), values(init(extent))) do grid, init
for j in 1:gridsize(grid)[2], i in 1:gridsize(grid)[1]
@inbounds source(grid)[i, j] = dest(grid)[i, j] = init[i, j]
end
updatestatus!(grid)
end
@set! simdata.extent = extent
@set! simdata.ruleset = ruleset
simdata
end
# Convert regular index to block index
indtoblock(x, blocksize) = (x - 1) ÷ blocksize + 1
# Convert block index to regular index
blocktoind(x, blocksize) = (x - 1) * blocksize + 1