Skip to content
Open
Changes from 2 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
68 changes: 61 additions & 7 deletions src/OutputWriters/windowed_time_average.jl
Original file line number Diff line number Diff line change
Expand Up @@ -102,17 +102,70 @@ Base.copy(sch::AveragedTimeInterval) = AveragedTimeInterval(sch.interval, window

A schedule for averaging over windows that precede SpecifiedTimes.
"""
mutable struct AveragedSpecifiedTimes <: AbstractSchedule
mutable struct AveragedSpecifiedTimes{W <: Union{Float64, Vector{Float64}}} <: AbstractSchedule
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
mutable struct AveragedSpecifiedTimes{W <: Union{Float64, Vector{Float64}}} <: AbstractSchedule
mutable struct AveragedSpecifiedTimes{W} <: AbstractSchedule

specified_times :: SpecifiedTimes
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since SpecifiedTimes is not a concrete type, this should also be a type parameter

window :: Float64
window :: W
stride :: Int
collecting :: Bool
end

const VaryingWindowAveragedSpecifiedTimes = AveragedSpecifiedTimes{Vector{Float64}}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
const VaryingWindowAveragedSpecifiedTimes = AveragedSpecifiedTimes{Vector{Float64}}
const VaryingWindowAveragedSpecifiedTimes = AveragedSpecifiedTimes{<:Vector}

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think you want to limit to Float64. That won't work with dates.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should we relax AveragedTimeInterval as well?

mutable struct AveragedTimeInterval <: AbstractSchedule
interval :: Float64
window :: Float64
stride :: Int
first_actuation_time :: Float64
actuations :: Int
collecting :: Bool
end

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes! good catch. I can do that in another PR if you like

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

in the above, interval and window have a "time period" type, while first_actuation_time has a "time" type. These are different when using dates.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is there currently a TimePeriod type? I can try to define one.

Should we also make a SpecifiedWindows like SpecifiedTimes in the case for AveragedSpecifiedTimes?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We don't want a specific type, because sometimes a "period" is Float64 (the assumption here), but with dates, it may need to be Dates.Second or Dates.Day (for example). However, a "time type" can be Float64 (the assumption here) or DateTime.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should we relax it and make it parametric? or constrain it so that interval :: Union{Float64, Dates.Period}

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Right, we need to use type parameters, and there need to be two of them.

We don't use the type union approach because it is not concrete; with that approach types can no longer be inferred.


AveragedSpecifiedTimes(specified_times::SpecifiedTimes; window, stride=1) =
AveragedSpecifiedTimes(specified_times, window, stride, false)

AveragedSpecifiedTimes(times; kw...) = AveragedSpecifiedTimes(SpecifiedTimes(times); kw...)
AveragedSpecifiedTimes(times; window, kw...) = AveragedSpecifiedTimes(times, window; kw...)
# AveragedSpecifiedTimes(times; window::Float64; kw...) = AveragedSpecifiedTimes(SpecifiedTimes(times); window, kw...)

function AveragedSpecifiedTimes(times, window::Vector{Float64}; kw...)
length(window) == length(times) || throw(ArgumentError("When providing a vector of windows, its length $(length(window)) must match the number of specified times $(length(times))."))
perm = sortperm(times)
sorted_times = times[perm]
sorted_window = window[perm]
time_diff = diff(vcat(0, sorted_times))

@info "timediff", time_diff
@info "sortedwindow", sorted_window

any(time_diff .- sorted_window .< -eps(eltype(window))) && throw(ArgumentError("Averaging windows overlap. Ensure that for each specified time tᵢ, tᵢ - windowᵢ ≥ tᵢ₋₁."))

return AveragedSpecifiedTimes(SpecifiedTimes(times); window, kw...)
end

function AveragedSpecifiedTimes(times, window::Float64; kw...)
perm = sortperm(times)
sorted_times = times[perm]
time_diff = diff(vcat(0, sorted_times))

any(time_diff .- window .< -eps(typeof(window))) && throw(ArgumentError("Averaging window $window is too large and causes overlapping windows. Ensure that for each specified time tᵢ, tᵢ - window ≥ tᵢ₋₁."))

return AveragedSpecifiedTimes(SpecifiedTimes(times); window, kw...)
end

# function AveragedSpecifiedTimes(times; window, kw...)
# perm = sortperm(times)
# sorted_times = times[perm]
# time_diff = diff(vcat(0, sorted_times))

# if window isa Vector{Float64}
# length(window) == length(times) || throw(ArgumentError("When providing a vector of windows, its length $(length(window)) must match the number of specified times $(length(times))."))

# sorted_window = window[perm]
# @info "timediff", time_diff
# @info "sortedwindow", sorted_window

# any(time_diff .- sorted_window .< -eps(eltype(window))) && throw(ArgumentError("Averaging windows overlap. Ensure that for each specified time tᵢ, tᵢ - windowᵢ ≥ tᵢ₋₁."))
# elseif window isa Number
# any(time_diff .- window .< -eps(typeof(window))) && throw(ArgumentError("Averaging window $window is too large and causes overlapping windows. Ensure that for each specified time tᵢ, tᵢ - window ≥ tᵢ₋₁."))
# else
# throw(ArgumentError("window must be a Float64 or a Vector{Float64}, got $(typeof(window))"))
# end

# return AveragedSpecifiedTimes(SpecifiedTimes(times); window=window, kw...)
# end

get_next_window(schedule::VaryingWindowAveragedSpecifiedTimes) = schedule.window[schedule.specified_times.previous_actuation + 1]
get_next_window(schedule::AveragedSpecifiedTimes) = schedule.window

function (schedule::AveragedSpecifiedTimes)(model)
time = model.clock.time
Expand All @@ -121,7 +174,7 @@ function (schedule::AveragedSpecifiedTimes)(model)
next > length(schedule.specified_times.times) && return false

next_time = schedule.specified_times.times[next]
window = schedule.window
window = get_next_window(schedule)

schedule.collecting || time >= next_time - window
end
Expand All @@ -132,7 +185,8 @@ function outside_window(schedule::AveragedSpecifiedTimes, clock)
next = schedule.specified_times.previous_actuation + 1
next > length(schedule.specified_times.times) && return true
next_time = schedule.specified_times.times[next]
return clock.time < next_time - schedule.window
window = get_next_window(schedule)
return clock.time < next_time - window
end

function end_of_window(schedule::AveragedSpecifiedTimes, clock)
Expand Down Expand Up @@ -171,7 +225,7 @@ stride(wta::SpecifiedWindowedTimeAverage) = wta.schedule.stride
WindowedTimeAverage(operand, model=nothing; schedule)

Returns an object for computing running averages of `operand` over `schedule.window` and
recurring on `schedule.interval`, where `schedule` is an `AveragedTimeInterval`.
recurring on `schedule.interval`, where `schedule` is an `AveragedTimeInterval` or `AveragedSpecifiedTimes`.
During the collection period, averages are computed every `schedule.stride` iteration.

`operand` may be a `Oceananigans.Field` or a function that returns an array or scalar.
Expand Down Expand Up @@ -275,7 +329,7 @@ function advance_time_average!(wta::SpecifiedWindowedTimeAverage, model)
# Begin collecting window-averaged increments
wta.schedule.collecting = true

wta.window_start_time = next_actuation_time(wta.schedule) - wta.schedule.window
wta.window_start_time = next_actuation_time(wta.schedule) - get_next_window(wta.schedule)
wta.previous_collection_time = wta.window_start_time
wta.window_start_iteration = model.clock.iteration - 1
# @info "t $(prettytime(model.clock.time)), next actuation time: $(prettytime(next_actuation_time(wta.schedule))), window $(prettytime(wta.schedule.window))"
Expand Down
Loading