Interpolated stats

src/Output.jl

@@ -61,21 +61,25 @@ show(io::IO, o::MemoryOutput) = print(io, "MemoryOutput$(collect(keys(o.data)))"
 """
 function (o::MemoryOutput)(y, t, dt, yfun)
     save, ts = o.save_cond(y, t, dt, o.saved)
-    append_stats!(o, o.statsfun(y, t, dt))
     !haskey(o.data, o.yname) && initialise(o, y)
     while save
+        ys = yfun(ts)
+        append_stats!(o, o.statsfun(ys, ts, dt))
         s = size(o.data[o.yname])
         if s[end] < o.saved+1
-            o.data[o.yname] = fastcat(o.data[o.yname], yfun(ts))
+            o.data[o.yname] = fastcat(o.data[o.yname], ys)
             push!(o.data[o.tname], ts)
         else
             idcs = fill(:, ndims(y))
-            o.data[o.yname][idcs..., o.saved+1] = yfun(ts)
+            o.data[o.yname][idcs..., o.saved+1] = ys
             o.data[o.tname][o.saved+1] = ts
         end
         o.saved += 1
         save, ts = o.save_cond(y, t, dt, o.saved)
     end
+    if ts != t
+        append_stats!(o, o.statsfun(y, t, dt))
+    end
 end
 
 function append_stats!(o::MemoryOutput, d)
@@ -209,7 +213,7 @@ function HDF5Output(fpath::AbstractString)
     HDF5Output(fpath, 0, 0, 1; readonly=true)
 end
 
-function initialise(o::HDF5Output, y)
+function initialise(o::HDF5Output, y, statsnames)
     ydims = size(y)
     idims = init_dims(ydims, o.save_cond)
     cdims = collect(idims)
@@ -228,7 +232,6 @@ function initialise(o::HDF5Output, y)
         end
         HDF5.create_dataset(file, o.tname, HDF5.datatype(Float64), ((dims[end],), (-1,)),
                       chunk=(1,))
-        statsnames = sort(collect(keys(o.stats_tmp[end])))
         o.cachehash = hash((statsnames, size(y)))
         file["meta"]["cachehash"] = o.cachehash
         if o.cache
@@ -300,22 +303,26 @@ end
 function (o::HDF5Output)(y, t, dt, yfun)
     o.readonly && error("Cannot add data to read-only output!")
     save, ts = o.save_cond(y, t, dt, o.saved)
-    push!(o.stats_tmp, o.statsfun(y, t, dt))
     if save
         @hlock HDF5.h5open(o.fpath, "r+") do file
-            !HDF5.haskey(file, o.yname) && initialise(o, y)
-            statsnames = sort(collect(keys(o.stats_tmp[end])))
-            cachehash = hash((statsnames, size(y)))
-            cachehash == o.cachehash || error(
-                "the hash for this propagation does not agree with cache in file")
+            statsnames = sort(collect(keys(o.statsfun(y, t, dt))))
+            if  !HDF5.haskey(file, o.yname)
+                initialise(o, y, statsnames)
+            else    
+                cachehash = hash((statsnames, size(y)))
+                cachehash == o.cachehash || error(
+                    "the hash for this propagation does not agree with cache in file")
+            end
             while save
                 s = collect(size(file[o.yname]))
                 idcs = fill(:, length(s)-1)
                 if s[end] < o.saved+1
                     s[end] += 1
                     HDF5.set_extent_dims(file[o.yname], Tuple(s))
                 end
-                file[o.yname][idcs..., o.saved+1] = yfun(ts)
+                ys = yfun(ts)
+                push!(o.stats_tmp, o.statsfun(ys, ts, dt))
+                file[o.yname][idcs..., o.saved+1] = ys
                 s = collect(size(file[o.tname]))
                 if s[end] < o.saved+1
                     s[end] += 1
@@ -335,6 +342,9 @@ function (o::HDF5Output)(y, t, dt, yfun)
             end
         end
     end
+    if ts != t
+        push!(o.stats_tmp, o.statsfun(y, t, dt))
+    end
 end
 
 function append_stats!(parent, a::Array{Dict{String,Any},1})

test/test_output.jl

@@ -193,6 +193,11 @@ fpath_comp = joinpath(dirpath, "test_comp.h5")
     ω, Eω, zac = Processing.getEω(o, 5e-2)
     @test (ω, Eω[:, 1]) == Processing.getEω(grid, mem["Eω"][:, 51])
     @test zac[1] == 5e-2
+    # test stat save locations
+    @test length(unique(mem["stats"]["z"])) == length(mem["stats"]["z"])
+    @test length(unique(o["stats"]["z"])) == length(o["stats"]["z"])
+    @test all([z in o["stats"]["z"] for z in o["z"]])
+    @test all([z in mem["stats"]["z"] for z in mem["z"]])
 end
 rm(fpath)
 rm(fpath_comp)