Merge pull request #32 from cesaraustralia/modeltuple

Modeltuple

Project.toml

@@ -1,7 +1,7 @@
 name = "GrowthMaps"
 uuid = "a252ced4-3752-11e9-0615-256422162c4d"
 authors = ["Rafael Schouten", "James Maino"]
-version = "0.2.0"
+version = "0.2.1"
 
 [deps]
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
@@ -18,9 +18,9 @@ UnitfulRecipes = "42071c24-d89e-48dd-8a24-8a12d9b8861f"
 [compat]
 ConstructionBase = "1"
 GeoData = "0.2"
-InteractModels = "0.2"
+InteractModels = "0.3"
 LsqFit = "0.10, 0.11"
-ModelParameters = "0.2"
+ModelParameters = "0.3"
 Plots = "1"
 Reexport = "0.2"
 Unitful = "1.0"

src/framework.jl

@@ -6,8 +6,10 @@ Combine growth rates accross layers and subperiods for all required periods.
 
 ## Arguments
 
-- `model`: A `Model` or `Tuple` of `Layer` components. 
-  `Layer`s can also be passed in as separate arguments.
+- `model`: A `Model`, NamedTuple of `Model`,  or `Tuple` of `Layer` components. 
+
+`Layer`s can also be passed in as separate arguments. A `NamedTuple` of models
+will all be run from the same data source, reducing load time.
 
 ## Keyword Arguments
 
@@ -18,18 +20,20 @@ Combine growth rates accross layers and subperiods for all required periods.
 The output is a GeoArray with the same dimensions as the passed in stack layers, and a Time
 dimension with a length of `nperiods`.
 """
-mapgrowth(model::Model; kwargs...) = mapgrowth(parent(model); kwargs...)
-mapgrowth(layers...; kwargs...) = mapgrowth(layers; kwargs...)
-function mapgrowth(layers::Tuple; 
+mapgrowth(layers::Layer...; kw...) = mapgrowth(layers; kw...)
+mapgrowth(layers::Tuple{<:Layer,Vararg}; kw...) = mapgrowth(Model(layers); kw...)
+mapgrowth(model::Model; kw...) = mapgrowth((model,); kw...)[1]
+mapgrowth(m1::Model, ms::Model...; kw...) = mapgrowth((m1, ms...); kw...)
+function mapgrowth(models::Union{Tuple{<:Model,Vararg},NamedTuple{<:Any,Tuple{<:Model,Vararg}}}; 
     series::AbstractGeoSeries, tspan::AbstractRange, arraytype=Array
 )
-    layers = stripparams(layers)
+    models = map(stripparams ∘ parent, models)
     period = step(tspan); nperiods = length(tspan)
     startdate, enddate = first(tspan), last(tspan)
-    required_keys = Tuple(union(keys(layers)))
+    required_keys = Tuple(union(map(keys ∘ _astuple, models)...))
 
     # Copy only the required keys to a memory-backed stack
-    stack = GeoStack(deepcopy(first(series)); keys=required_keys)
+    stack = deepcopy(GeoStack(first(series); keys=required_keys))
 
     A = first(values(stack));
     missingval = eltype(A)(NaN)
@@ -42,15 +46,17 @@ function mapgrowth(layers::Tuple;
     ti = Ti(tspan; mode=Sampled(Ordered(), Regular(period), Intervals(Start())))
     outdims = (dims(A)..., ti)
     outA = arraytype(zeros(eltype(A), size(A)..., nperiods))
-    output = GeoArray(outA, outdims; name=:growthrate, missingval=missingval)
+    outputs = map(models) do m
+        GeoArray(deepcopy(outA), outdims; name=:growthrate, missingval=missingval)
+    end
 
-    runperiods!(output, stackbuffer, series, mask, layers, tspan)
+    runperiods!(outputs, stackbuffer, series, mask, models, tspan)
 
     # Return a GeoArray wrapping a regular Array, not arraytype
-    GeoData.modify(Array, output)
+    map(o -> GeoData.modify(Array, o), outputs)
 end
 
-function runperiods!(output, stackbuffer, series, mask, layers, tspan)
+function runperiods!(outputs, stackbuffer, series, mask, models, tspan)
     period = step(tspan); nperiods = length(tspan)
     println("Running for $(1:nperiods)")
     for p in 1:nperiods
@@ -67,15 +73,21 @@ function runperiods!(output, stackbuffer, series, mask, layers, tspan)
             println("    ", val(dims(subseries, Ti))[t])
             # Copy the arrays we need from disk to the buffer stack
             copy!(stackbuffer, subseries[t])
-            # For some reason now this is broken with DD getindex, view is a workaround
-            parent(view(output, Ti(p))) .+= combinelayers(layers, stackbuffer)
+            map(outputs, models, keys(outputs)) do output, model, key
+                # For some reason now this is broken with DD getindex, view is a workaround
+                parent(view(output, Ti(p))) .+= combinelayers(model, stackbuffer)
+            end
             n += 1
         end
         if n > 0
-            parent(view(output, Ti(p))) .*= mask ./ n
+            map(outputs) do output
+                parent(view(output, Ti(p))) .*= mask ./ n
+            end
         else
             @warn ("No files found for the $period period starting $periodstart")
-            parent(view(output, Ti(p))) .*= mask
+            map(outputs) do output
+                parent(view(output, Ti(p))) .*= mask
+            end
         end
     end
 end
@@ -100,3 +112,6 @@ maybetoK(val) = unit(val) == u"°C" ? val |> K : val
     conditionalrate(model(l), maybetoK(val))
 @inline conditionalrate(model::RateModel, val) =
     condition(model, val) ? rate(model, val) : zero(rate(model, val))
+
+_astuple(x::Tuple) = x
+_astuple(x) = (x,)