text:
##
# Sparse BroadcastStyle
##

for Typ in (:Diagonal, :SymTridiagonal, :Tridiagonal, :Bidiagonal)
    @eval begin
        BroadcastStyle(::StructuredMatrixStyle{<:$Typ}, ::BandedStyle) =
            BandedStyle()
        BroadcastStyle(::BandedStyle, ::StructuredMatrixStyle{<:$Typ}) =
            BandedStyle()
    end
end

# Here we implement the banded matrix interface for some key examples
isbanded(::Zeros) = true
bandwidths(::Zeros) = (-40320,-40320) # 40320 == prod(1:8), used for special cases involving gcd
inbands_getindex(::Zeros{T}, k::Integer, j::Integer) where T = zero(T)

isbanded(::Eye) = true
bandwidths(::Eye) = (0,0)
inbands_getindex(::Eye{T}, k::Integer, j::Integer) where T = one(T)

isbanded(::Diagonal) = true
bandwidths(::Diagonal) = (0,0)
inbands_getindex(D::Diagonal, k::Integer, j::Integer) = D.diag[k]
inbands_setindex!(D::Diagonal, v, k::Integer, j::Integer) = (D.diag[k] = v)
bandeddata(D::Diagonal) = permutedims(D.diag)

# treat subinds as banded
sublayout(::DiagonalLayout{L}, inds::Type) where L = sublayout(bandedcolumns(L()), inds)
sublayout(::DiagonalLayout{L}, inds::Type{<:NTuple{2,AbstractUnitRange{Int}}}) where L = sublayout(bandedcolumns(L()), inds)

# bandeddata(V::SubArray{<:Any,2,<:Diagonal}) = view(bandeddata(parent(V)), :, parentindices(V)[2])

isbanded(::Bidiagonal) = true
bandwidths(A::Bidiagonal) = A.uplo == 'U' ? (0,1) : (1,0)

sublayout(::BidiagonalLayout, ::Type{<:Tuple{AbstractUnitRange{Int},AbstractUnitRange{Int}}}) =
    BandedLayout()

isbanded(::SymTridiagonal) = true
bandwidths(::SymTridiagonal) = (1,1)
inbands_getindex(J::SymTridiagonal, k::Integer, j::Integer) =
    k == j ? J.dv[k] : J.ev[min(k,j)]
inbands_setindex!(J::SymTridiagonal, v, k::Integer, j::Integer) =
    k == j ? (J.dv[k] = v) : (J.ev[min(k,j)] = v)

isbanded(::Tridiagonal) = true
bandwidths(::Tridiagonal) = (1,1)

sublayout(::AbstractTridiagonalLayout, ::Type{<:Tuple{AbstractUnitRange{Int},AbstractUnitRange{Int}}}) =
    BandedLayout()

###
# rot180
###

function rot180(A::AbstractBandedMatrix)
    m,n = size(A)
    sh = m-n
    l,u = bandwidths(A)
    _BandedMatrix(bandeddata(A)[end:-1:1,end:-1:1], m, u+sh,l-sh)
end

function getindex(D::Diagonal{T,V}, b::Band) where {T,V}
    iszero(b.i) && return copy(D.diag)
    convert(V, Zeros{T}(size(D,1)-abs(b.i)))
end

function getindex(D::Tridiagonal{T,V}, b::Band) where {T,V}
    b.i == -1 && return copy(D.dl)
    iszero(b.i) && return copy(D.d)
    b.i == 1 && return copy(D.du)
    convert(V, Zeros{T}(size(D,1)-abs(b.i)))
end

function getindex(D::SymTridiagonal{T,V}, b::Band) where {T,V}
    iszero(b.i) && return copy(D.dv)
    abs(b.i) == 1 && return copy(D.ev)
    convert(V, Zeros{T}(size(D,1)-abs(b.i)))
end

function getindex(D::Bidiagonal{T,V}, b::Band) where {T,V}
    iszero(b.i) && return copy(D.dv)
    D.uplo == 'L' && b.i == -1 && return copy(D.ev)
    D.uplo == 'U' && b.i == 1 && return copy(D.ev)
    convert(V, Zeros{T}(size(D,1)-abs(b.i)))
end


text:
format_defs = Dict()

format_defs[:heading] = Dict(
	"bold" => true,
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "center",
	"border" => true
)

format_defs[:text] = Dict(
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "left",
	"border" => true
)

format_defs[:heading_right] = Dict(
	"bold" => true,
	"font" => "Arial",
	"size" => 10,
	"valign" => "vcenter",
	"align" => "right",
	"left" => true, "bottom" => true, "top" => true
)

format_defs[:heading_left] = Dict(
	"bold" => true,
	"font" => "Arial",
	"size" => 10,
	"valign" => "vcenter",
	"align" => "left",
	"right" => true, "bottom" => true, "top" => true
)

format_defs[:model_name] = Dict(
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "left",
	"border" => true
)

format_defs[:varname_1indent] = Dict(
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "left",
	"border" => true,
	"indent" => 1
)

format_defs[:varname_2indent] = Dict(
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "left",
	"border" => true,
	"indent" => 2
)

format_defs[:n_fmt_right] = Dict(
	"num_format" => "#,##0",
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "right",
	"left" => true, "bottom" => true, "top" => true
)

format_defs[:n_fmt_left_parens] = Dict(
	"num_format" => "(#,##0)",
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "right",
	"right" => true, "bottom" => true, "top" => true
)

format_defs[:f_date] = Dict(
	"num_format" => "mm/dd/yyyy",
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "right",
	"border" => true
)

format_defs[:f_datetime] = Dict(
	"num_format" => "mm/dd/yyyy hh:mm:ss",
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "right",
	"border" => true
)

format_defs[:f_fmt_right] = Dict(
	"num_format" => "#,##0.00",
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "right",
	"left" => true, "bottom" => true, "top" => true
)

format_defs[:f_fmt_center] = Dict(
	"num_format" => "#,##0.0000",
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "center",
	"border" => true
)

format_defs[:f_fmt] = Dict(
	"num_format" => "#,##0.00",
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "right",
	"border" => true
)

format_defs[:n_fmt] = Dict(
	"num_format" => "#,##0",
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "right",
	"border" => true
)

format_defs[:n_fmt_center] = Dict(
	"num_format" => "#,##0",
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "center",
	"border" => true
)

format_defs[:f_fmt_left_parens] = Dict(
	"num_format" => "(#,##0.00)",
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "left",
	"right" => true, "bottom" => true, "top" => true
)

format_defs[:pct_fmt_parens] = Dict(
	"num_format" => "(0.00%)",
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "left",
	"right" => true, "bottom" => true, "top" => true
)

format_defs[:pct_fmt] = Dict(
	"num_format" => "0.00%",
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "right",
	"border" => true
)

format_defs[:or_fmt] = Dict(
	"num_format" => "0.000",
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "right",
	"left" => true, "bottom" => true, "top" => true
)

format_defs[:cilb_fmt] = Dict(
	"num_format" => "(0.000 -",
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "right",
	"bottom" => true, "top" => true
)

format_defs[:ciub_fmt] = Dict(
	"num_format" => "0.000)",
	"font" => "Arial",
	"font_size" => 9,
	"valign" => "vcenter",
	"align" => "left",
	"right" => true, "bottom" => true, "top" => true
)

format_defs[:or_fmt_red] = Dict(
	"num_format" => "0.000",
	"font" => "Arial",
	"font_size" => 9,
	"font_color" => "red",
	"valign" => "vcenter",
	"align" => "right",
	"left" => true, "bottom" => true, "top" => true
)

format_defs[:cilb_fmt_red] = Dict(
	"num_format" => "(0.000 -",
	"font" => "Arial",
	"font_size" => 9,
	"font_color" => "red",
	"valign" => "vcenter",
	"align" => "right",
	"bottom" => true, "top" => true
)

format_defs[:ciub_fmt_red] = Dict(
	"num_format" => "0.000)",
	"font" => "Arial",
	"font_size" => 9,
	"font_color" => "red",
	"valign" => "vcenter",
	"align" => "left",
	"right" => true, "bottom" => true, "top" => true
)

format_defs[:or_fmt_bold] = Dict(
	"num_format" => "0.000",
	"font" => "Arial",
	"font_size" => 9,
	"bold" => true,
	"valign" => "vcenter",
	"align" => "right",
	"left" => true, "bottom" => true, "top" => true
)

format_defs[:cilb_fmt_bold] = Dict(
	"num_format" => "(0.000 -",
	"font" => "Arial",
	"font_size" => 9,
	"bold" => true,
	"valign" => "vcenter",
	"align" => "right",
	"bottom" => true, "top" => true
)

format_defs[:ciub_fmt_bold] = Dict(
	"num_format" => "0.000)",
	"font" => "Arial",
	"font_size" => 9,
	"bold" => true,
	"valign" => "vcenter",
	"align" => "left",
	"right" => true, "bottom" => true, "top" => true
)

format_defs[:p_fmt] = Dict(
	"num_format" => "0.000",
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "right",
	"border" => true
)

format_defs[:p_fmt_center] = Dict(
	"num_format" => "0.000",
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "center",
	"border" => true
)

format_defs[:p_fmt2] = Dict(
	#"num_format" => "0.000",
	"font" => "Arial",
	"size" => 9,
	"valign" => "vcenter",
	"align" => "right",
	"border" => true
)

format_defs[:p_fmt_red] = Dict(
	"num_format" => "0.000",
	"font" => "Arial",
	"font_size" => 9,
	"font_color" => "red",
	"valign" => "vcenter",
	"align" => "right",
	"border" => true
)

format_defs[:p_fmt2_red] = Dict(
	#"num_format" => "0.000",
	"font" => "Arial",
	"font_size" => 9,
	"font_color" => "red",
	"valign" => "vcenter",
	"align" => "right",
	"border" => true
)

format_defs[:p_fmt_bold] = Dict(
	"num_format" => "0.000",
	"font" => "Arial",
	"font_size" => 9,
	"bold" => true,
	"valign" => "vcenter",
	"align" => "right",
	"border" => true
)

format_defs[:p_fmt2_bold] = Dict(
	#"num_format" => "0.000",
	"font" => "Arial",
	"font_size" => 9,
	"bold" => true,
	"valign" => "vcenter",
	"align" => "right",
	"border" => true
)

format_defs[:empty_border] = Dict(
	#"num_format" => "0.000",
	"valign" => "vcenter",
	"border" => true
)

format_defs[:empty_left] = Dict(
	#"num_format" => "0.000",
	"valign" => "vcenter",
	"right" => true, "bottom" => true, "top" => true
)

format_defs[:empty_right] = Dict(
	#"num_format" => "0.000",
	"valign" => "vcenter",
	"left" => true, "bottom" => true, "top" => true
)

format_defs[:empty_both] = Dict(
	#"num_format" => "0.000",
	"valign" => "vcenter",
	"bottom" => true, "top" => true
)

function attach_formats(workbook;formats::Dict = format_defs)
	newfmts = Dict()
	for key in keys(formats)
		newfmts[key] = workbook.add_format(formats[key])
	end
	return newfmts
end

ext: jl


text:
module Events

import ..Dagger
import ..Dagger: Event, Chunk
import ..Dagger: init_similar

import ..Tables

"""
    CoreMetrics

Tracks the timestamp, category, and kind of the `Event` object generated by log
events.
"""
struct CoreMetrics end

(::CoreMetrics)(ev::Event{T}) where T =
    (;timestamp=ev.timestamp,
      category=ev.category,
      kind=T)

"""
    IDMetrics

Tracks the ID of `Event` objects generated by log events.
"""
struct IDMetrics end

(::IDMetrics)(ev::Event{T}) where T = ev.id

"""
    TimelineMetrics

Tracks the timeline of `Event` objects generated by log events.
"""
struct TimelineMetrics end

(::TimelineMetrics)(ev::Event{T}) where T = ev.timeline

"""
    FullMetrics

Tracks the full `Event` object generated by log events.
"""
struct FullMetrics end

(::FullMetrics)(ev) = ev

"""
    BytesAllocd

Tracks memory allocated for `Chunk`s.
"""
mutable struct BytesAllocd
    allocd::Int
end
BytesAllocd() = BytesAllocd(0)
init_similar(::BytesAllocd) = BytesAllocd()

function (ba::BytesAllocd)(ev::Event{:start})
    if ev.category in (:move, :evict) && ev.timeline.data isa Chunk
        sz = Int(ev.timeline.data.handle.size)
        if ev.category == :move && !haskey(Dagger.Sch.CHUNK_CACHE, ev.timeline.data)
            ba.allocd += sz
        elseif ev.category == :evict && haskey(Dagger.Sch.CHUNK_CACHE, ev.timeline.data)
            ba.allocd -= sz
        end
    end
    ba.allocd
end
(ba::BytesAllocd)(ev::Event{:finish}) = ba.allocd

"""
    CPULoadAverages

Monitors the CPU load averages.
"""
struct CPULoadAverages end

(::CPULoadAverages)(ev::Event) = Sys.loadavg()[1]

"""
    MemoryFree

Monitors the percentage of free system memory.
"""
struct MemoryFree end

(::MemoryFree)(ev::Event) = Sys.free_memory() / Sys.total_memory()

"""
    EventSaturation

Tracks the compute saturation (running tasks) per-processor.
"""
mutable struct EventSaturation
    saturation::Dict{Symbol,Int}
end
EventSaturation() = EventSaturation(Dict{Symbol,Int}())
init_similar(::EventSaturation) = EventSaturation()

function (es::EventSaturation)(ev::Event{:start})
    old = get(es.saturation, ev.category, 0)
    es.saturation[ev.category] = old + 1
    NamedTuple(filter(x->x[2]>0, es.saturation))
end
function (es::EventSaturation)(ev::Event{:finish})
    old = get(es.saturation, ev.category, 0)
    es.saturation[ev.category] = old - 1
    NamedTuple(filter(x->x[2]>0, es.saturation))
end

"""
    WorkerSaturation

Tracks the compute saturation (running tasks).
"""
mutable struct WorkerSaturation
    saturation::Int
end
WorkerSaturation() = WorkerSaturation(0)
init_similar(::WorkerSaturation) = WorkerSaturation()

function (ws::WorkerSaturation)(ev::Event{:start})
    if ev.category == :compute
        ws.saturation += 1
    end
    ws.saturation
end
function (ws::WorkerSaturation)(ev::Event{:finish})
    if ev.category == :compute
        ws.saturation -= 1
    end
    ws.saturation
end

"""
    ProcessorSaturation

Tracks the compute saturation (running tasks) per-processor.
"""
mutable struct ProcessorSaturation
    saturation::Dict{Dagger.Processor,Int}
end
ProcessorSaturation() = ProcessorSaturation(Dict{Dagger.Processor,Int}())
init_similar(::ProcessorSaturation) = ProcessorSaturation()

function (ps::ProcessorSaturation)(ev::Event{:start})
    if ev.category == :compute
        proc = ev.timeline.to_proc
        old = get(ps.saturation, proc, 0)
        ps.saturation[proc] = old + 1
    end
    filter(x->x[2]>0, ps.saturation)
end
function (ps::ProcessorSaturation)(ev::Event{:finish})
    if ev.category == :compute
        proc = ev.timeline.to_proc
        old = get(ps.saturation, proc, 0)
        ps.saturation[proc] = old - 1
    end
    filter(x->x[2]>0, ps.saturation)
end

"Debugging metric, used to log event start/finish via `@debug`."
struct DebugMetrics
    sat::EventSaturation
end
DebugMetrics() = DebugMetrics(EventSaturation())
function (dm::DebugMetrics)(ev::Event{T}) where T
    dm.sat(ev)
    sat_string = "$(join(["$cat => $(dm.sat.saturation[cat])" for cat in keys(dm.sat.saturation) if dm.sat.saturation[cat] > 0 ], ", "))"
    @debug "Event: $(ev.category) $T ($sat_string)" _line=nothing _file=nothing _module=nothing
    nothing
end
init_similar(dm::DebugMetrics) = DebugMetrics()

"""
    LogWindow

Aggregator that prunes events to within a given time window.
"""
mutable struct LogWindow
    window_length::UInt64
    core_name::Symbol
    creation_handlers::Vector{Any}
    deletion_handlers::Vector{Any}
end
LogWindow(window_length, core_name) = LogWindow(window_length, core_name, [], [])

function (lw::LogWindow)(logs::Dict)
    core_log = logs[lw.core_name]

    # Inform creation hooks
    if length(core_log) > 0
        log = Dict{Symbol,Any}()
        for key in keys(logs)
            log[key] = last(logs[key])
        end
        for obj in lw.creation_handlers
            creation_hook(obj, log)
        end
    end

    # Find entries outside the window
    window_start = time_ns() - lw.window_length
    idx = something(findfirst(ev->ev.timestamp >= window_start, core_log),
                    length(core_log))

    # Return if all events are in the window
    if length(core_log) == 0 || (core_log[1].timestamp >= window_start)
        return
    end

    # Inform deletion hooks
    for obj in lw.deletion_handlers
        deletion_hook(obj, idx)
    end

    # Remove entries outside the window
    for name in keys(logs)
        deleteat!(logs[name], 1:idx)
    end
end

creation_hook(x, log) = nothing
deletion_hook(x, idx) = nothing

"""
    TableStorage

LogWindow-compatible aggregator which stores logs in a Tables.jl-compatible
sink.

Using a `TableStorage` is reasonably simple:

```julia
ml = Dagger.MultiEventLog()

... # Add some events

lw = Dagger.Events.LogWindow(5*10^9, :core)

# Create a DataFrame with one Any[] for each event
df = DataFrame([key=>[] for key in keys(ml.consumers)]...)

# Create the TableStorage and register its creation handler
ts = Dagger.Events.TableStorage(df)
push!(lw.creation_handlers, ts)

ml.aggregators[:lw] = lw

# Logs will now be saved into `df` automatically, and packages like
# DaggerWebDash.jl will automatically use it to retrieve subsets of the logs.
```
"""
struct TableStorage{T}
    sink::T
    function TableStorage(sink::T) where T
        @assert Tables.istable(sink)
        new{T}(sink)
    end
end
init_similar(ts::TableStorage) = TableStorage(similar(ts.sink, 0))
function creation_hook(ts::TableStorage, log)
    try
        push!(ts.sink, NamedTuple(log))
    catch err
        rethrow(err)
    end
end
function Base.getindex(ts::TableStorage, ts_range::UnitRange)
    ts_low, ts_high = ts_range.start, ts_range.stop
    return filter(row->ts_low <= row.core.timestamp <= ts_high,
                  Tables.rows(ts.sink)) |> Tables.materializer(ts.sink)
end

end # module Events


text:
using Compat

# Glyph has to be a Bokehjs type as it's defined directly in Bokeh's JSON, but these methods are
# more porcelain than plumbing, so are defined here.
function Bokehjs.Glyph(glyphtype::Symbol,
               linecolor::NullString,
               linewidth::NullInt,
               linealpha::NullFloat,
               fillcolor::NullString,
               fillalpha::NullFloat,
               size::NullInt,
               dash::Union(Nothing, Vector{Int64}),
               fields::Union(Nothing, Dict{Symbol, Symbol}))
    props = @Compat.compat Dict(
        :linecolor => linecolor == nothing ? omit : @Compat.compat(Dict{Symbol, BkAny}(:value => linecolor)),
        :linewidth => linewidth == nothing ? omit : @Compat.compat(Dict{Symbol, BkAny}(:units => :data, :value =>linewidth)),
        :linealpha => linealpha == nothing ? omit : @Compat.compat(Dict{Symbol, BkAny}(:units => :data, :value =>linealpha)),
        :fillalpha => fillalpha == nothing ? omit : @Compat.compat(Dict{Symbol, BkAny}(:units => :data, :value =>fillalpha)),
        :size => size == nothing ? omit : @Compat.compat(Dict{Symbol, BkAny}(:units => :screen, :value =>size)),
        :fillcolor => fillcolor == nothing ? omit : @Compat.compat(Dict{Symbol, BkAny}(:value =>fillcolor)),
    )
    if fields != nothing
        for (field, val) in fields
            !haskey(props, field) && error("unknown field $(field) passed to Glyph")
            props[field] = @Compat.compat Dict{Symbol, BkAny}(:field => val, :units => :data)
        end
    end
    Glyph(Bokehjs.uuid4(),
          glyphtype,
          props[:linecolor],
          props[:linewidth],
          props[:linealpha],
          props[:fillcolor],
          props[:fillalpha],
          props[:size],
          dash == nothing ? omit : dash,
          @Compat.compat(Dict(:units =>:data, :field => :x)),
          @Compat.compat(Dict(:units =>:data, :field => :y)))
end

function Bokehjs.Glyph(;glyphtype=nothing,
               linecolor=nothing,
               linewidth=nothing,
               linealpha=nothing,
               fillcolor=nothing,
               fillalpha=nothing,
               size=nothing,
               dash=nothing,
               fields=nothing)
    glyphtype = glyphtype == nothing ? :Line : glyphtype
    Glyph(glyphtype, linecolor, linewidth, linealpha, fillcolor, fillalpha, size, dash, fields)
end

function Bokehjs.Glyph(glyphtype::Symbol; kwargs...)
    Glyph(glyphtype=glyphtype; kwargs...)
end

function Base.show(io::IO, g::Bokehjs.Glyph)
    names = Glyph.names
    features = String[]
    for name in Glyph.names
        showname = name == :_type_name ? :type : name
        g.(name) != nothing && push!(features, "$showname: $(g.(name))")
    end
    print(io, "Glyph(", join(features, ", "), ")")
end

type BokehDataSet
    data::Dict{Symbol, Vector}
    glyph::Glyph
    legend::NullString
    BokehDataSet(data::Dict{Symbol, Vector}, glyph::Glyph, legend::NullString=nothing) = new(data, glyph, legend)
end

function BokehDataSet(xdata::RealVect, ydata::RealVect, args...)
    data = @compat Dict{Symbol, Vector}(:x => xdata, :y => ydata)
    BokehDataSet(data, args...)
end

type Plot
    datacolumns::Array{BokehDataSet, 1}
    tools::Vector{Symbol}
    filename::String
    title::String
    width::Int
    height::Int
    x_axis_type::NullSymbol
    y_axis_type::NullSymbol
    legendsgo::NullSymbol
end

# both the string and the reversed string will be tried, eg. "ox" is equivilent to "xo"
STRINGTOKENS = @compat Dict(
    "--"=>@Compat.compat(Dict(:dash=>[4, 4])),
    "-."=>@Compat.compat(Dict(:dash=>[1, 4, 2])),
    "ox"=>@Compat.compat(Dict(:glyphtype=>:CircleX)),
    "o+"=>@Compat.compat(Dict(:glyphtype=>:CircleCross)),
    "sx"=>@Compat.compat(Dict(:glyphtype=>:SquareX)),
    "s+"=>@Compat.compat(Dict(:glyphtype=>:SquareCross)),
)

for (k, v) in Dict(STRINGTOKENS)
    STRINGTOKENS[reverse(k)] = v
end

# heavily borrowed from Winston, thanks Winston!
const CHARTOKENS = @compat Dict(
    '-'=>@Compat.compat(Dict(:dash=>nothing)),
    ':'=>@Compat.compat(Dict(:dash=>[1, 4])),
    ';'=>@Compat.compat(Dict(:dash=>[1, 4, 2])),
    '+'=>@Compat.compat(Dict(:glyphtype=>:Cross)),
    'o'=>@Compat.compat(Dict(:glyphtype=>:Circle)),
    '*'=>@Compat.compat(Dict(:glyphtype=>:Asterisk)),
    '.'=>@Compat.compat(Dict(:glyphtype=>:Circle, :size=>2)),
    'x'=>@Compat.compat(Dict(:glyphtype=>:X)),
    's'=>@Compat.compat(Dict(:glyphtype=>:Square)),
    'd'=>@Compat.compat(Dict(:glyphtype=>:Diamond)),
    '^'=>@Compat.compat(Dict(:glyphtype=>:Triangle)),
    'v'=>@Compat.compat(Dict(:glyphtype=>:InvertedTriangle)),
    'y'=>@Compat.compat(Dict(:linecolor=>"yellow")),
    'm'=>@Compat.compat(Dict(:linecolor=>"magenta")),
    'c'=>@Compat.compat(Dict(:linecolor=>"cyan")),
    'r'=>@Compat.compat(Dict(:linecolor=>"red")),
    'g'=>@Compat.compat(Dict(:linecolor=>"green")),
    'b'=>@Compat.compat(Dict(:linecolor=>"blue")),
    'w'=>@Compat.compat(Dict(:linecolor=>"white")),
    'k'=>@Compat.compat(Dict(:linecolor=>"black")),
)

Base.convert(::Type{Array{Glyph, 1}}, glyph::Glyph) = [glyph]

function Base.convert(::Type{Array{Glyph, 1}}, styles::String)
    map(style -> convert(Glyph, style), split(styles, '|'))
end

function Base.convert(::Type{Glyph}, style::String)
    styd = @compat Dict(:glyphtype=>:Line, :linecolor=>"blue", :linewidth=>1, :linealpha=>1.0)
    for key in keys(STRINGTOKENS)
        splitstyle = split(style, key)
        if length(splitstyle) > 1
            for (k, v) in STRINGTOKENS[key]
                styd[k] = v
            end
            style = join(splitstyle)
        end
    end

    for char in style
        if haskey(CHARTOKENS, char)
            for (k, v) in CHARTOKENS[char]
                styd[k] = v
            end
        else
            warn("unrecognized char '$char'")
        end
    end

    filledglyphs = [:Circle, :Square, :Diamond, :Triangle, :InvertedTriangle]
    if in(styd[:glyphtype], filledglyphs)
        styd[:fillcolor] = styd[:linecolor]
        styd[:fillalpha] = DEFAULT_FILL_ALPHA
        # this seems to be the best way of making plots look right, better ideas?
        styd[:linealpha] = DEFAULT_FILL_ALPHA
        if !haskey(styd, :size)
            styd[:size] = DEFAULT_SIZE
        end
    end
    emptyglyphs = [:CircleX, :CircleCross, :SquareX, :SquareCross]
    if in(styd[:glyphtype], emptyglyphs)
        styd[:fillcolor] = "transparent"
        styd[:size] = DEFAULT_SIZE
    end
    Glyph(;styd...)
end

max_stars_repo_name: Laeeth/dbokeh

max_issues_repo_name: Laeeth/dbokeh

max_forks_repo_name: Laeeth/dbokeh


text:
read_cstring(io) = String(readuntil(io, 0x00))

function read_elem(io::IO)
    T = bson_type[read(io, UInt8)]
    Element(read_cstring(io), read_value(io, T))
end

function read_elist(io::IO)
    head = position(io)
    len = read(io, Int32)
    @assert len >= 0
    tail = head + len - 1
    ret = Element[]
    while position(io) < tail
        push!(ret, read_elem(io))
    end
    @assert position(io) == tail
    @assert read(io, UInt8) == 0x00
    ret
end

read_value(io::IO, ::Type{Missing}) = missing
read_value(io::IO, ::Type{Nothing}) = nothing
read_value(io::IO, ::Type{Maxkey}) = maxkey
read_value(io::IO, ::Type{Minkey}) = minkey
read_value(io::IO, ::Type{Bool}) = Bool(read(io, UInt8))
read_value(io::IO, ::Type{Int32}) = read(io, Int32)
read_value(io::IO, ::Type{Int64}) = read(io, Int64)
read_value(io::IO, ::Type{Float64}) = read(io, Float64)
read_value(io::IO, ::Type{DateTime}) = DateTime(UTM(UNIXEPOCH + read(io, Int64)))
read_value(io::IO, ::Type{BSONSymbol}) = BSONSymbol(read_value(io, String))
read_value(io::IO, ::Type{Timestamp}) = Timestamp(read(io, UInt64))
read_value(io::IO, ::Type{ObjectId}) = ObjectId(read(io, 12))
read_value(io::IO, ::Type{Decimal128}) = Decimal128(read(io, UInt128))
read_value(io::IO, ::Type{DBPointer}) = DBPointer(read_value(io, String), read_value(io, ObjectId))
read_value(io::IO, ::Type{Code}) = Code(read_value(io, String))
read_value(io::IO, ::Type{BSONRegex}) = BSONRegex(read_cstring(io), read_cstring(io))
read_value(io::IO, ::Type{Document}) = Document(read_elist(io))
read_value(io::IO, ::Type{BSONArray})= BSONArray(read_elist(io))

function read_value(io, ::Type{String})
    len = read(io, Int32)
    @assert len >= 0
    ret = String(read(io, len - 1))
    @assert isvalid(ret)
    @assert read(io, UInt8) == 0x00
    ret
end

function read_value(io::IO, ::Type{Binary})
    len = read(io, Int32)
    @assert len >= 0
    subtype = read(io, UInt8)
    if subtype == 0x02
        len = read(io, Int32)
        @assert len >= 0
    end
    ret = Binary(subtype, read(io, len))
end

function read_value(io::IO, ::Type{CodeWithScope})
    head = position(io)
    len = read(io, Int32)
    @assert len >= 0
    tail = head + len
    code = read_value(io, String)
    doc = read_value(io, Document)
    @assert position(io) == tail
    CodeWithScope(code, doc)
end

# high level functions

function read_bson(io::IO; codec::Union{AbstractBSONCodec,Symbol} = :bson)
    codec = codec isa Symbol ? _codec[codec] : codec
    raw = read_value(io, Document)
    @assert(eof(io))
    decode(codec, raw)
end

read_bson(buf::AbstractVector{UInt8}; kwargs...) = read_bson(IOBuffer(buf, read = true, write = false); kwargs...)

read_bson(filename::AbstractString; kwargs...) = open(x -> read_bson(x; kwargs...), filename, "r")


text:
abstract type AbstractAPIResponder end

struct APISpec
    fn::Function
    resp_json::Bool
    resp_headers::Dict{String,String}
end

const EndPts = Dict{String,APISpec}

"""
APIResponder holds the transport and format used for data exchange and the endpoint specifications.
"""
struct APIResponder{T<:AbstractTransport,F<:AbstractMsgFormat}
    transport::T
    format::F
    id::Union{Nothing,String}  # optional responder id to be sent back
    open::Bool  #whether the responder will process all functions, or only registered ones
    endpoints::EndPts
end

APIResponder(transport::T, format::F, id::Union{Nothing,String}=nothing, open::Bool=false) where {T<:AbstractTransport,F<:AbstractMsgFormat} = APIResponder(transport, format, id, open, EndPts())

"""
APIResponder holds the transport and format used for data exchange and the endpoint specifications.
This method creates an APIResponder over ZMQ transport using JSON message format.
(provided for backward compatibility)
"""
APIResponder(addr::String, ctx::Context=Context(), bound::Bool=true, id=nothing, open=false) = APIResponder(ZMQTransport(addr, REP, bound, ctx), JSONMsgFormat(), id, open, EndPts())
APIResponder(ip::IPv4, port::Int, ctx::Context=Context()) = APIResponder("tcp://$ip:$port", ctx)

function Base.show(io::IO, x::APIResponder)
    print(io, "JuliaWebAPI.APIResponder with $(length(x.endpoints)) endpoints (", join(keys(x.endpoints), ","), ")")
end

function default_endpoint(f::Function)
    endpt = string(f)
    # separate the module (more natural URL, assumes 'using Module')
    if '.' in endpt
        endpt = rsplit(endpt, '.', limit=2)[2]
    end
    endpt
end

"""
Register a function as API call.
TODO: validate method belongs to module?
"""
function register(conn::APIResponder, f::Function;
                  resp_json::Bool=false,
                  resp_headers::Dict=Dict{String,String}(), endpt=default_endpoint(f))
    @info("registering", endpt)
    conn.endpoints[endpt] = APISpec(f, resp_json, resp_headers)
    return conn # make fluent api possible
end

"""send a response over the transport in the specified format"""
function respond(conn::APIResponder, code::Int, headers::Dict, resp)
    resp = wireformat(conn.format, code, headers, resp, conn.id)
    sendresp(conn.transport, resp)
end

respond(conn::APIResponder, api::Union{Nothing,APISpec}, status::Symbol, resp=nothing) =
    respond(conn, ERR_CODES[status][1], get_hdrs(api), get_resp(api, status, resp))

get_hdrs(api::Nothing) = Dict{String,String}()
get_hdrs(api::APISpec) = api.resp_headers

function get_resp(api::Union{Nothing,APISpec}, status::Symbol, resp=nothing)
    st = ERR_CODES[status]
    stcode = st[2]
    stresp = ((stcode != 0) && (resp === nothing)) ? string(st[3], " : ", st[2]) : resp

    if (api !== nothing) && api.resp_json
        return Dict{String, Any}("code"=>stcode, "data"=>stresp)
    else
        return stresp
    end
end

# Note: needs to be changed when https://github.com/JuliaLang/julia/pull/22646 is merged
function dynamic_invoke(conn::APIResponder, f, args...; kwargs...)
    if conn.open && isdefined(Core, :_apply_latest)
        inner() = f(args...; kwargs...)
        Core._apply_latest(inner)
    else
        f(args...; kwargs...)
    end
end

"""call the actual API method, and send the return value back as response"""
function call_api(api::APISpec, conn::APIResponder, args, data::Dict{Symbol,Any})
    try
        if !applicable(api.fn, args...) || ((api.fn === (*) || api.fn === (/) || api.fn === (\)) && all(x->isa(x,Vector), args))
            narrow_args!(args)
        end
        result = dynamic_invoke(conn, api.fn, args...; data...)
        respond(conn, api, :success, result)
    catch ex
        logerr("api_exception: ", ex)
        respond(conn, api, :api_exception, string(ex))
    end
end

#=
narrow_args! if a private function that processes an Array transferred
over JSON to a properly typed and dimensioned Julia array. Arrays
serialised from JSON are stored as an Array{Any}, even if all the elements
are members of the same concrete type, eg, Float64. This function will
transform such an array to an Array{Float64} type.

Further, since JSON does not have true multidimensional arrays, they are
transmitted as arrays containing arrays. This function will convert them to
a true multidimensional array in Julia.
=#
promote_arr(x) = Base.collect_to!(similar(x, typeof(x[1])), x, 1, 1)

function narrow_args!(x)
    for (i, v) in enumerate(x)
        if (typeof(v) <: AbstractArray)
            if (length(v) > 0 && typeof(v[1]) <: Array)
                x[i] = hcat(x[i]...)
            end
            x[i] = promote_arr(x[i])
        end
    end
end

"""start processing as a server"""
function process(conn::APIResponder; async::Bool=false)
    if async
        @debug("processing async...")
        @async process(conn)
    else
        @debug("processing...")
        while true
            msg = juliaformat(conn.format, recvreq(conn.transport))

            command = cmd(conn.format, msg)
            @info("received", command)

            if startswith(command, ':')    # is a control command
                ctrlcmd = Symbol(command[2:end])
                if ctrlcmd === :terminate
                    respond(conn, nothing, :terminate, "")
                    break
                else
                    @error("invalid control command ", command)
                    continue
                end
            end

            if !haskey(conn.endpoints, command)
                if !conn.open || !isdefined(Main, Symbol(command))
                    respond(conn, nothing, :invalid_api)
                    continue
                else
                    _add_spec(getfield(Main, Symbol(command)), conn)
                end
            end

            try
                call_api(conn.endpoints[command], conn, args(conn.format, msg), data(conn.format, msg))
            catch ex
                logerr("exception ", ex)
                respond(conn, conn.endpoints[command], :invalid_data)
            end
        end
        close(conn.transport)
        @info("stopped processing.")
    end
    conn
end

_add_spec(fn::Function, api::APIResponder) = register(api, fn, resp_json=false, resp_headers=Dict{String,String}())

function _add_spec(spec::Tuple, api::APIResponder)
    fn = spec[1]
    resp_json = (length(spec) > 1) ? spec[2] : false
    resp_headers = (length(spec) > 2) ? spec[3] : Dict{String,String}()
    api_name = (length(spec) > 3) ? spec[4] : default_endpoint(fn)
    register(api, fn, resp_json=resp_json, resp_headers=resp_headers, endpt=api_name)
end

function create_responder(apispecs::Array, addr, bind, nid, open=false)
    api = APIResponder(ZMQTransport(addr, REP, bind, Context()), JSONMsgFormat(), nid, open)
    for spec in apispecs
        _add_spec(spec, api)
    end
    api
end

function logerr(msg, ex)
    iob = IOBuffer()
    write(iob, msg)
    showerror(iob, ex)
    @error(String(take!(iob)))
end


text:
using Azure
using Test

# on public CI, just testing whether package loads at the moment


text:
export CuDim3, CuDim

"""
    CuDim3(x)

    CuDim3((x,))
    CuDim3((x, y))
    CuDim3((x, y, x))

A type used to specify dimensions, consisting of 3 integers for respectively the `x`, `y`
and `z` dimension. Unspecified dimensions default to `1`.

Often accepted as argument through the `CuDim` type alias, eg. in the case of
[`cudacall`](@ref) or [`CUDA.launch`](@ref), allowing to pass dimensions as a plain integer
or a tuple without having to construct an explicit `CuDim3` object.
"""
struct CuDim3
    x::Cuint
    y::Cuint
    z::Cuint
end

CuDim3(dims::Integer)             = CuDim3(dims,    Cuint(1), Cuint(1))
CuDim3(dims::NTuple{1,<:Integer}) = CuDim3(dims[1], Cuint(1), Cuint(1))
CuDim3(dims::NTuple{2,<:Integer}) = CuDim3(dims[1], dims[2],  Cuint(1))
CuDim3(dims::NTuple{3,<:Integer}) = CuDim3(dims[1], dims[2],  dims[3])

# Type alias for conveniently specifying the dimensions
# (e.g. `(len, 2)` instead of `CuDim3((len, 2))`)
const CuDim = Union{Integer,
                    Tuple{Integer},
                    Tuple{Integer, Integer},
                    Tuple{Integer, Integer, Integer}}


text:
# Command: Delete Session
"""
	delete!(session::Session)
Deletes the session from the Remote Driver.
"""
function delete!(session::Session)
    try
        @unpack addr, id = session
        response =
            HTTP.delete("$addr/session/$id", [("Content-Type" => "application/json")])
        @assert response.status == 200
        @assert JSON3.read(response.body).status == 0
        id
    catch err
        throw(WDError(err))
    end
end


text:
__precompile__(true)

module EFTfitter

using BAT
using Distributions
using LinearAlgebra
using StatsBase
using NamedTupleTools
using Setfield
using Parameters
using RecipesBase
using Requires
using ValueShapes

include("datatypes.jl")
include("EFTfitterModel.jl")
include("EFTfitterDensity.jl")
include("ranking/ranking.jl")
include("plotting/plotting.jl")
include("utils.jl")

const _PLOTS_MODULE = Ref{Union{Module,Nothing}}(nothing)
_plots_module() = _PLOTS_MODULE[]

function __init__()
    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" _PLOTS_MODULE[] = Plots
end

end # module


text:
module NLP
    using JuMP
    using Ipopt
    using Base

    function find_optimal_physical_parameters(GraphParam, ConfParam, output_file, type_energy, print_function = print)
        """
        Solves a nonlinear optimization of the network's physical parameters
        Hypothesis: temperature of the fluid in the supply and return networks
        is assumed to be constant

        INPUT:
        - GraphParam, ConfParam, output_dir

        OUTPUT:
        - Diameter (dict): pipe(u,v) -> diameter of the pipe in meters
        - AnnualCost (scalar): equivalent annual cost of the network
        """

        ########################################################################
        # Parameters
        ########################################################################
        LENGTH = GraphParam["LENGTH"]
        ELEVATION = GraphParam["ELEVATION"]
        OUTFLOW = GraphParam["OUTFLOW"]
        SUPPLY_MAX_INFLOW =  GraphParam["SUPPLY_MAX_INFLOW"]
        LB_DIAM = GraphParam["LB_DIAM"]
        MAX_CAPACITY = GraphParam["MAX_CAPACITY"]

        # DIAMETER RANGE
        MIN_DIAMETER = ConfParam["MIN_DIAMETER"]
        MAX_DIAMETER = ConfParam["MAX_DIAMETER"]

        # PRESSURE RANGE
        MIN_PRESSURE = ConfParam["MIN_PRESSURE"]
        MAX_PRESSURE = ConfParam["MAX_PRESSURE"]

        # THERM PARAM
        if type_energy == "Heating"
            SUPPLY_HC_TEMPERATURE = ConfParam["SUPPLY_HEAT_TEMPERATURE"]
            RETURN_HC_TEMPERATURE = ConfParam["RETURN_HEAT_TEMPERATURE"]
            HC_PERIOD = ConfParam["HEATING_PERIOD"]
            HC_COST = ConfParam["HEAT_COST"]
        elseif type_energy == "Cooling"
            SUPPLY_HC_TEMPERATURE = ConfParam["SUPPLY_COOL_TEMPERATURE"]
            RETURN_HC_TEMPERATURE = ConfParam["RETURN_COOL_TEMPERATURE"]
            HC_PERIOD = ConfParam["COOLING_PERIOD"]
            HC_COST = ConfParam["COOL_COST"]
        end
        GROUND_TEMPERATURE = ConfParam["GROUND_TEMPERATURE"]

        # ECO PARAM
        ELECTRICITY_COST = ConfParam["ELECTRICITY_COST"]
        INTEREST_RATE = ConfParam["INTEREST_RATE"]
        LIFETIME = ConfParam["LIFETIME"]
        DEPRECIATION_MAINTENANCE_RATE = ConfParam["DEPRECIATION_MAINTENANCE_RATE"]
        EFFECT_COEF = INTEREST_RATE / (1 - (1 + INTEREST_RATE)^(-LIFETIME)) + DEPRECIATION_MAINTENANCE_RATE

        # PUMPING PARAM
        PUMPING_EFFICIENCY = ConfParam["PUMPING_EFFICIENCY"]

        # PHYSICAL PARAM
        CP = ConfParam["CP"]
        RHO = ConfParam["RHO"]
        ROUGHNESS = ConfParam["ROUGHNESS"]
        FRICTION_COEF = (2 * log(ROUGHNESS / (3.71 * (MIN_DIAMETER + MAX_DIAMETER)/2 * 1e-2)) / log(10))^(-2)  #Fanning factor

        # REGRESSIONS
        A_MAX_VELOCITY = ConfParam["A_MAX_VELOCITY"]
        B_MAX_VELOCITY = ConfParam["B_MAX_VELOCITY"]

        A_LINEAR_COST = ConfParam["A_LINEAR_COST"]
        B_LINEAR_COST = ConfParam["B_LINEAR_COST"]

        if type_energy == "Heating"
            A_TRANSIT_COEF = ConfParam["A_HEAT_TRANSIT_COEF"]
            B_TRANSIT_COEF = ConfParam["B_HEAT_TRANSIT_COEF"]
        elseif type_energy == "Cooling"
            A_TRANSIT_COEF = ConfParam["A_COOL_TRANSIT_COEF"]
            B_TRANSIT_COEF = ConfParam["B_COOL_TRANSIT_COEF"]
        end

        # CONVERSION
        W_TO_MW = 1e-6
        DAY_TO_HOUR = 24
        BAR_TO_PA = 1e5
        CM_TO_M = 1e-2
        CM_TO_MM = 10
        G = 9.81 # gravitation constant (m/s2)

        # Sets of edges and nodes
        PIPE = keys(LENGTH)
        NODE = [n for n=Set(vcat([p[1] for p in PIPE], [p[2] for p in PIPE]))]
        SUPPLY_NODE = keys(SUPPLY_MAX_INFLOW)
        BUILDING_NODE = keys(OUTFLOW)

        # Set the bounds for mass flows at each node
        LOWER_BOUND_IO = Dict(n => 0.0 for n in NODE)
        UPPER_BOUND_IO = Dict(n => 0.0 for n in NODE)
        for n in BUILDING_NODE
            LOWER_BOUND_IO[n] = - OUTFLOW[n]
            UPPER_BOUND_IO[n] = - OUTFLOW[n]
        end
        for n in SUPPLY_NODE
            LOWER_BOUND_IO[n] = 0.0
            UPPER_BOUND_IO[n] = SUPPLY_MAX_INFLOW[n]
        end

        #print_function("1")
        ########################################################################
        # Optimization
        ########################################################################

        m = Model(with_optimizer(Ipopt.Optimizer, output_file = output_file,
                                        mu_strategy = "adaptive",  # more successfull
                                        tol = 1e-3, # tolerance hasn't to be very high
                                        max_cpu_time=1.0*ConfParam["MAXTIME"])) # time limit

        ############################ Variables #################################
        @variable(m, MIN_DIAMETER <= varDiameter[PIPE] <= MAX_DIAMETER) # Pipe diameters (cm)
        @variable(m, varVelocity[PIPE]) # Velocity in pipes (m/s),
        # Positive flow (oriented edges) p = (n1, n2) <=> the flow is going from n1 to n2
        @variable(m, varMassFlow[PIPE]) # Flow in pipes (kg/s)
        @variable(m, MIN_PRESSURE <= varPressure[PIPE] <= MAX_PRESSURE) # Pressure in pipe (bars)
        @variable(m, varPressureFriction[NODE]) # The pressure drop due to friction
        @variable(m, varConstBernoulli) # variable that plays the role of the constant in the Bernoulli equation

        ############################ Constraints ##############################
        for p in PIPE

            # MassFlow = rho * section * velecity where section = pi*r^2
            @NLconstraint(m, 4*varMassFlow[p] == RHO * pi * (varDiameter[p] * CM_TO_M)^2 * varVelocity[p])

            # Max Velocity
            @constraint(m, varVelocity[p] <= A_MAX_VELOCITY * varDiameter[p] * CM_TO_MM + B_MAX_VELOCITY)
            @constraint(m, -varVelocity[p] <= A_MAX_VELOCITY * varDiameter[p] * CM_TO_MM + B_MAX_VELOCITY)

            # Pressure loss in pipes (just due to friction)
            @NLconstraint(m, (varPressureFriction[p[1]]-varPressureFriction[p[2]]) * BAR_TO_PA == 0.5 * RHO *
                LENGTH[p] * varVelocity[p] * abs(varVelocity[p]) * FRICTION_COEF / (varDiameter[p]*CM_TO_M))

            # Bernoulli's equation
            @NLconstraint(m, varPressure[p]*BAR_TO_PA + 1/2*RHO*varVelocity[p]^2 + RHO*G*(ELEVATION[p[1]]+ELEVATION[p[2]])/2 ==
                varConstBernoulli + (varPressureFriction[p[1]]+varPressureFriction[p[2]])/2*BAR_TO_PA)

            # Diameter to respect
            @constraint(m, varDiameter[p] >= get(LB_DIAM, p, 0))
            @constraint(m, varDiameter[p] >= get(LB_DIAM, (p[2],p[1]), 0))
            if p in keys(MAX_CAPACITY)
                if !isinf(MAX_CAPACITY[p]) && MAX_CAPACITY[p] >= MIN_DIAMETER
                    # need to be at least equals to MIN_DIAMETER, we put 1 cm more than real value
                    @constraint(m, varDiameter[p] <= get(LB_DIAM, p, Inf)+1)
                    @constraint(m, varDiameter[p] <= get(LB_DIAM, (p[2],p[1]), Inf)+1)
                end
            end
        end

        # Flow conservation
        for n in NODE
            @constraint(m, sum(varMassFlow[p] for p=PIPE if p[1]==n) - sum(varMassFlow[p] for p=PIPE if p[2]==n) >= LOWER_BOUND_IO[n])
            @constraint(m, sum(varMassFlow[p] for p=PIPE if p[1]==n) - sum(varMassFlow[p] for p=PIPE if p[2]==n) <= UPPER_BOUND_IO[n])
        end

        ############################## Objective ###############################

         @objective(m, Min, (A_LINEAR_COST * CM_TO_MM * EFFECT_COEF + # Construction cost
            2 * A_TRANSIT_COEF * CM_TO_MM * W_TO_MW * (HC_PERIOD * DAY_TO_HOUR) * HC_COST ) * # Heat loss cost
            sum(varDiameter[p] * LENGTH[p] for p=PIPE)

            # Pumping cost (the problem become much longer with pumping cost)
            + sum(SUPPLY_MAX_INFLOW[n] * sum(varPressure[p] for p=PIPE if n in p) for n=SUPPLY_NODE) *
            BAR_TO_PA / RHO * W_TO_MW / PUMPING_EFFICIENCY * (HC_PERIOD * DAY_TO_HOUR) * ELECTRICITY_COST)

        ########################################################################
        # Solve and collect the solution
        ########################################################################
        optimize!(m)
        status = termination_status(m)

        NLP_Output = Dict()
        #print_function("2")
        if status == MOI.OPTIMAL
            NLP_Output["Diameter"] = Dict(p => value(varDiameter[p]) for p=PIPE)
            NLP_Output["Velocity"] = Dict(p => value(varVelocity[p]) for p=PIPE)
            NLP_Output["MassFlow"] = Dict(p => value(varMassFlow[p]) for p=PIPE)
            NLP_Output["Pressure"] = Dict(p => value(varPressure[p]) for p=PIPE)
            NLP_Output["PressureFriction"] = Dict(n => value(varPressureFriction[n]) for n=NODE)
            NLP_Output["TotalCost"] = getobjectivevalue(m)

            NLP_Output["ConstructionCost"] = Dict(p => (A_LINEAR_COST * value(varDiameter[p]) * CM_TO_MM + B_LINEAR_COST) * LENGTH[p] * EFFECT_COEF for p=PIPE)
            if type_energy == "Heating"
                NLP_Output["HeatLossCost"] = Dict(p => 2 * (A_TRANSIT_COEF * value(varDiameter[p]) * CM_TO_MM + B_TRANSIT_COEF) * LENGTH[p] *
                                    W_TO_MW * (HC_PERIOD * DAY_TO_HOUR) * HC_COST for p=PIPE)
            elseif type_energy == "Cooling"
                NLP_Output["CoolLossCost"] = Dict(p => 2 *(A_TRANSIT_COEF * value(varDiameter[p]) * CM_TO_MM + B_TRANSIT_COEF) * LENGTH[p] *
                                    W_TO_MW * (HC_PERIOD * DAY_TO_HOUR) * HC_COST for p=PIPE)
            end
            NLP_Output["PumpingCost"] = Dict(n => value(sum(varPressure[p] for p=PIPE if n in p) - 1) * BAR_TO_PA *
                SUPPLY_MAX_INFLOW[n] / RHO * W_TO_MW / PUMPING_EFFICIENCY * (HC_PERIOD * DAY_TO_HOUR) * ELECTRICITY_COST for n=SUPPLY_NODE)

        end
        status_str = "Unkown"
        if status == MOI.OPTIMAL
            status_str = "Optimal"
        elseif status == MOI.INFEASIBLE
            status_str = "Infeasible"
        end
        return NLP_Output, status_str
    end

    export find_optimal_physical_parameters
end

max_stars_repo_name: Planheat/Planheat-Tool

max_issues_repo_name: Planheat/Planheat-Tool

max_forks_repo_name: Planheat/Planheat-Tool


text: zerofuel(distancetopump, mpg, fuelleft) = distancetopump/mpg <= fuelleft

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
using DelimitedFiles

const eye3 = [true false false;
              false true false;
              false false true]

"""
    writetrc(filename, f; <keyword arguments>)

Write the C3DFile `f` to a .trc format at `filename`.

# Keyword arguments
- `delim::Char='\\t'`: The text delimiter to use
- `strip_prefixes::Bool=true`: Strip marker label prefixes if they exist
- `subject::String=""`: The subject (among multiple subjects) in the C3DFile to write
- `prefixes::Vector{String}=[subject]`: Marker label prefixes to strip if
  `strip_prefixes == true`
- `remove_unlabeled_markers::Bool=true`: Remove markers with empty labels, labels
  matching `r"(\\*\\d+|M\\d\\d\\d)"`
- `lab_orientation::Matrix=Matrix(I, (3,3))`: Rotation to apply to markers before writing
- `precision::Int=6`: Number of decimal places to print
- `virtual_markers::Dict{String,Matrix}`: (Optional) Virtual markers (calculated with markers
  from `f`) to write to the .trc file
"""
function writetrc(filename::String, f::C3DFile;
    delim::Char='\t',
    # TODO: Add units kwarg and predicate precision on units
    precision::Int=6,
    strip_prefixes::Bool=true,
    remove_unlabeled_markers::Bool=true,
    subject::String="",
    prefixes::Vector{String}=[subject],
    lab_orientation::AbstractMatrix{T}=eye3,
    virtual_markers::Dict{String,Matrix{U}}=Dict{String,Matrix{Float32}}()
) where {T,U}
    open(filename, "w") do io
        writetrc(io, f; delim, precision, strip_prefixes, remove_unlabeled_markers, subject,
            prefixes, lab_orientation, virtual_markers)
    end
end

function writetrc(io, f::C3DFile;
    delim::Char='\t',
    precision::Int=6,
    strip_prefixes::Bool=true,
    remove_unlabeled_markers::Bool=true,
    subject::String="",
    prefixes::Vector{String}=[subject],
    lab_orientation::AbstractMatrix{T}=eye3,
    virtual_markers::Dict{String,Matrix{U}}=Dict{String,Matrix{Float32}}()
) where {T,U}
    if subject !== ""
        if haskey(f.groups, :SUBJECTS)
            any(subject .== f.groups[:SUBJECTS][Vector{String}, :NAMES]) || throw(ArgumentError("subject $subject does not exist in $f.groups[:SUBJECTS]"))
        elseif strip_prefixes && prefixes == [subject]
            @warn "subject $subject does not exist in $f.groups[:SUBJECTS] and may not be the correct prefix"
        end
    end

    len = (f.groups[:POINT][Int, :FRAMES] == typemax(UInt16)) ?
        f.groups[:POINT][Int, :LONG_FRAMES] : f.groups[:POINT][Int, :FRAMES]
    period = inv(f.groups[:POINT][Float64, :RATE])

    mkrnames = collect(keys(f.point))
    if subject !== ""
        filter!(x -> startswith(x, subject), mkrnames)
        isempty(mkrnames) && @warn "no markers matched subject $subject"
    end

    if strip_prefixes
        # Assume that LABEL_PREFIXES are used if present despite absence of USES_PREFIXES
        # prompted by a c3d file from OptiTrack Motive 2.2.0
        if haskey(f.groups, :SUBJECTS) && ((haskey(f.groups[:SUBJECTS], :USES_PREFIXES) &&
            f.groups[:SUBJECTS][Int, :USES_PREFIXES] == 1) ||
            (!haskey(f.groups[:SUBJECTS], :USES_PREFIXES) &&
            haskey(f.groups[:SUBJECTS], :LABEL_PREFIXES)))
            if subject !== ""
                subi = findfirst(==(subject), f.groups[:SUBJECTS][Vector{String}, :NAMES])
                r = Regex("("*f.groups[:SUBJECTS][Vector{String}, :LABEL_PREFIXES][subi] *
                        ")(?<label>\\w*)")
                mkrnames_stripped = replace.(mkrnames, r => s"\g<label>")
            else
                r = Regex("("*join([f.groups[:SUBJECTS][Vector{String}, :LABEL_PREFIXES];
                        prefixes], '|')*")(?<label>\\w*)")
                mkrnames_stripped = replace.(mkrnames, r => s"\g<label>")
            end
        elseif subject !== "" || prefixes !== [""]
            if any(subject .== prefixes)
                r = Regex("("*join(prefixes, '|')*")(?<label>\\w*)")
            else
                r = Regex("("*join([subject; prefixes], '|')*")(?<label>\\w*)")
            end
            mkrnames_stripped = replace.(mkrnames, r => s"\g<label>")
        else
            mkrnames_stripped = mkrnames
        end
    else
        mkrnames_stripped = mkrnames
    end

    if remove_unlabeled_markers
        ids = findall(x -> isempty(x) ||
            match(r"(\*\d+|M\d\d\d)", x) !== nothing, mkrnames)
        ids_strip = findall(x -> isempty(x) ||
            match(r"(\*\d+|M\d\d\d)", x) !== nothing, mkrnames_stripped)

        # mkrnames and mkrnames_stripped need to maintain the same order
        @assert ids == ids_strip

        deleteat!(mkrnames, ids)
        deleteat!(mkrnames_stripped, ids)
    end

    ord = sortperm(mkrnames_stripped)
    mkrnames = mkrnames[ord]
    mkrnames_stripped = mkrnames_stripped[ord]

    # Header
    line1 = ["PathFileType", "4", "(X/Y/Z)", basename(f.name)*'\n']
    line2 = ["DataRate", "CameraRate", "NumFrames", "NumMarkers", "Units", "OrigDataRate",
             "OrigDataStartFrame", "OrigNumFrames\n"]
    join(io, line1, delim)
    join(io, line2, delim)
    print(io, f.groups[:POINT][Float32, :RATE], delim)
    print(io, f.groups[:POINT][Float32, :RATE], delim)
    print(io, len, delim)
    print(io, length(mkrnames) + length(virtual_markers), delim)
    print(io, strip(f.groups[:POINT][String, :UNITS], Char(0x00)), delim)
    print(io, f.groups[:POINT][Float32, :RATE], delim)
    print(io, 1, delim)
    print(io, len, delim, '\n')

    # Header line
    write(io, string("Frame#", delim, "Time", delim))
    join(io, mkrnames_stripped, delim^3)
    print(io, delim^3)
    if !isempty(virtual_markers)
        extra_mkrnames = sort!(collect(keys(virtual_markers)))
        join(io, extra_mkrnames, delim^3)
        print(io, delim^3)
    else
        extra_mkrnames = String[]
    end
    println(io)

    # Coordinate line
    write(io, delim^2)
    join(io, (string('X', n, delim, 'Y', n, delim, 'Z', n)
              for n in 1:(length(mkrnames)+length(virtual_markers))), delim)
    write(io, '\n'^2)

    # Core data block
    et = promote_type(Float32, T, U)
    nanet = convert(et, NaN)
    nummkr = length(mkrnames)
    numxmkr = length(extra_mkrnames)
    data = Matrix{Union{Missing,et}}(undef, len, 1+3*(nummkr + numxmkr))

    data[:,1] .= range(zero(et), step=period, length=len)
    for (i, name) in enumerate(mkrnames)
        idxs = (2:4).+(i-1)*3
        data[:, idxs] .= coalesce.(f.point[name]*lab_orientation, nanet)
    end
    if !isempty(virtual_markers)
        for (i, name) in enumerate(extra_mkrnames)
            idxs = (2:4).+(i-1+nummkr)*3
            data[:, idxs] .= coalesce.(virtual_markers[name]*lab_orientation, nanet)
        end
    end

    data .= round.(data; digits=precision)
    data = Any[ 1:len data fill("", len) ]
    replace!(data, missing=>"")
    replace!(data, NaN=>"")

    writedlm(io, data, delim)

    return nothing
end


text:
################### acquisition functions for active learning ##################
random_sampling(G::Gaussian, x::AbstractVector) = rand(x)
random_sampling(G::Gaussian, x::AbstractMatrix) = rand(eachcol(x))

# choose the x with the highest uncertainty, according to G
function uncertainty_sampling(G::Gaussian, x::AbstractVector)
    σ = zeros(length(x))
    @threads for i in eachindex(x)
        σ[i] = std(G)(x[i])
    end
    i = argmax(σ)
    xs = x[i]
    return xs
end

# returns integrated uncertainty at integration nodes x with weights w
function integrated_uncertainty(G::Gaussian, x::AbstractVector, w::AbstractVector = fill(1/length(x), length(x)))
    dot(w, var(G, x))
end

function integrated_uncertainty_sampling(σ::Real)
    (G, x) -> integrated_uncertainty_sampling(G, x, σ)
end

# could introduce weights
function integrated_uncertainty_sampling(G::Gaussian, x::AbstractVector, σ::Real, w = ones(size(x)))
    n = length(x)
    K = Kernel.gramian(cov(G), x)
    K = Matrix(K) # this is the bottleneck
    u = zeros(n) # posterior uncertainty
    for (i, k) in enumerate(eachcol(K))
        u[i] = -sum(abs2, k) / (K[i, i] + σ^2) # technically have to add tr(K), but it doesn't change the minimum
    end
    return x[argmin(u)]
end

# ternary mul
using WoodburyIdentity
function Base.:*(x::AbstractMatrix{<:Real}, W::Woodbury, y::AbstractVecOrMat)
    x*(W*y)
    # xU = x*W.U
    # Vy = (x ≡ y' && W.U ≡ W.V') ? xU' : W.V*y
    # *(x, W.A, y) + W.α * *(xU, W.C, Vy) # can avoid two temporaries
end

############### acquisition function for SARA's inner loop #####################
# picks one-step Bayes optimal point to minimize derivative uncertainty
function inner_sampling(G::Gaussian, x::AbstractVector, σ::Real)
    n = length(x)
    D = GaussianDistributions.value_derivative(G) # multi-ouput process
    K = Kernel.gramian(cov(D), x)
    K11, K22, K21 = zeros(n, n), zeros(n, n), zeros(n, n)
    @threads for i in 1:n # derivative w.r.t. first component
        for j in 1:n
            kij = K.A[i, j] # ijth block
            K11[i, j] = kij[1, 1]
            K21[i, j] = kij[2, 1]
            K22[i, j] = kij[2, 2]
        end
    end
    u = zeros(size(x)) # uncertainty in derivative
    for i in 1:n
        u[i] = - sum(abs2, @view K21[:, i]) / (K11[i, i] + σ^2) # + K22
    end
    return x[argmin(u)]
end
inner_sampling(σ::Real) = (G, x) -> inner_sampling(G, x, σ)

################ upper-confidence-bound sampling ###############################
ucb_acquisition(G::Gaussian, α::Real) = x -> mean(G)(x) + α * std(G)(x)
ucb_acquisition(α::Real) = (G::Gaussian, x) -> ucb_acquisition(G, α)(x)

# choose a point among x with highest confidence bound score
# if α > 0 we are sample at points with high uncertainty
# if α = 0 it reduces to optimization of the predictive mean
# the larger α > 0, the more we exploit known regions with high function value
function ucb(G::Gaussian, x::AbstractVector, α::Real)
    a = ucb_acquisition(G, α)
    score = zeros(length(x))
    @threads for i in eachindex(x)
        score[i] = a(x[i])
    end
    i = argmax(score)
    return x[i]
end
ucb_sampling(α) = (G, x) -> ucb(G, x, α)

# optimizes x
# function ucb(G::Gaussian, x, α::Real, optimize::Val{true})
#     score(x) = mean(M)(x) + α * std(M)(x)
#     # IDEA: maximize score with respect to x
#     return -1
# end

####################### outer loop sampling strategies #########################
# modified sampling policy, taking into account that SARA
# is measuring entire lines in T (stripes) through lateral gradient LSA.
# G is result of get_gradient_map calculated from inner loop data,
# x is a vector of potential (T_max, log10_τ)-conditions
# returns most promising condition
function stripe_sampling(G::Gaussian, x::AbstractVector{<:Tuple}, relevant_T, acquisition)
    u = zeros(length(x)) # acquisition value
    @threads for i in eachindex(x)
        T_max, log10_τ = x[i]
        x_other = x[i][3:end] # other dimensions (e.g. composition)
        T = relevant_T(T_max, log10_τ)
        ci = tuple.(T, log10_τ, x_other...)
        u[i] = sum(x->acquisition(G, x), ci) # IDEA: could further parallize this sum
    end
    return x[argmax(u)]
end

# if X is a matrix, interpret columns as conditions to be measured NOT T peak conditions
function stripe_sampling(G::Gaussian, X::AbstractMatrix{<:Tuple}, acquisition)
    n = size(X, 2)
    u = zeros(n) # acquisition value
    @threads for i in 1:n
        xi = @view X[:, i]
        u[i] = sum(x->acquisition(G, x), xi) # IDEA: could further parallize this sum
    end
    return x[argmax(u)]
end

# returns outer_sampling acquisition function
# relevant_T calculates the relevant temperatures for each stripe condition
function stripe_sampling(relevant_T, acquisition)
    function (G::Gaussian, conditions::AbstractVector{<:Tuple})
        stripe_sampling(G, conditions, relevant_T, acquisition)
    end
end
# uncertainty and upper confidence bound sampling generalized to stripes
stripe_uncertainty_sampling(relevant_T) = stripe_sampling(relevant_T, (G, x)->var(G)(x))
function stripe_ucb_sampling(relevant_T, α::Real)
    stripe_sampling(relevant_T, ucb_acquisition(α))
end

################################################################################
# relevant temperatures with constant offset from T_max
function get_relevant_T(T_max::Real, log10_τ::Real, constant::Val{true},
                        c_low::Real, c_high::Real, n::Int,)
    collect(range(T_max - c_low, T_max - c_high, length = n))
end

# relevant temperatures with proportional offset from T_max
function get_relevant_T(T_max::Real, log10_τ::Real, constant::Val{false},
                        p_min::Real, p_max::Real, n::Int)
    collect(T_max * range(p_min, p_max, length = n))
end

# lazifies evaluation on T_max and log10_τ
function get_relevant_T(constant::Union{Val{true}, Val{false}}, c_min::Real, c_max::Real, n::Int)
    get_T(T_max::Real, log10_τ::Real) = get_relevant_T(T_max, log10_τ, constant, c_min, c_max, n)
end

# makes constant offset the default
function get_relevant_T(c_min::Real, c_max::Real, n::Int)
    constant = Val(true)
    get_T(T_max::Real, log10_τ::Real) = get_relevant_T(T_max, log10_τ, constant, c_min, c_max, n)
end


text:
@testset "Sampler" begin
    Random.seed!(0)
    d, l = hcat(randn(5, 95), randn(5, 5) .+ 2), fill(:inlier, 100)
    l[end-4:end] .= :outlier
    gamma = 0.5
    c = KDECache(d, gamma)
    chunk_mask = trues(length(l))
    chunk_mask[1:5] .= false
    threshold_strat = MaxDensityThresholdStrategy(0.5)
    threshold_strat_gt = GroundTruthThresholdStrategy()

    init_strat = SVDD.SimpleCombinedStrategy(SVDD.FixedGammaStrategy(MLKernels.GaussianKernel(gamma)),
                                             SVDD.FixedCStrategy(0.5))

    @testset "RandomSampler" begin
        @test_throws ArgumentError RandomRatioSampler(-1)
        s = RandomRatioSampler(1)
        sample_mask = sample(s, d, l)
        @test length(sample_mask) == length(l)
        @test all(sample_mask)
        s = RandomRatioSampler(0.5)
        sample_mask = sample(s, d, l)
        @test length(sample_mask) == length(l)
        @test !all(sample_mask)
    end

    @testset "RandomNSampler" begin
        @test_throws ArgumentError RandomNSampler(-1)
        s = RandomNSampler(length(l))
        sample_mask = sample(s, d, l)
        @test length(sample_mask) == length(l)
        @test all(sample_mask)
        s = RandomNSampler(length(l) - 5)
        sample_mask = sample(s, d, l)
        @test length(sample_mask) == length(l)
        @test !all(sample_mask)
    end

    @testset "BPS" begin
        @test_throws ArgumentError BPS(-1, 0.1)
        @test_throws ArgumentError BPS(0, 0.1)
        @test_throws ArgumentError BPS(5, -1)
        s = BPS(nothing, 0.1)
        sample_mask = sample(s, d, l)
        @test length(sample_mask) == length(l)
        @test !all(sample_mask)
    end

    @testset "CSSVDD" begin
        @test_throws ArgumentError CSSVDD(-1, 1, 1, 1, 0.05, init_strat, TEST_SOLVER)
        @test_throws ArgumentError CSSVDD(0, 1, 1, 1, 0.05, init_strat, TEST_SOLVER)
        @test_throws ArgumentError CSSVDD(1, -1, 1, 1, 0.05, init_strat, TEST_SOLVER)
        @test_throws ArgumentError CSSVDD(1, 0, 1, 1, 0.05, init_strat, TEST_SOLVER)
        @test_throws ArgumentError CSSVDD(1, 1, -1, 1, 0.05, init_strat, TEST_SOLVER)
        @test_throws ArgumentError CSSVDD(1, 1, 0, 1, 0.05, init_strat, TEST_SOLVER)
        @test_throws ArgumentError CSSVDD(1, 1, 1, -1, 0.05, init_strat, TEST_SOLVER)
        @test_throws ArgumentError CSSVDD(1, 1, 1, 0, 0.05, init_strat, TEST_SOLVER)
        @test_throws ArgumentError CSSVDD(1, 1, 1, -1, -0.05, init_strat, TEST_SOLVER)
        @test_throws ArgumentError CSSVDD(1, 1, 1, 0, 1.05, init_strat, TEST_SOLVER)
        s = CSSVDD(0.1, 0.05, init_strat, TEST_SOLVER)
        sample_mask = sample(s, d, l)
        @test length(sample_mask) == length(l)
        @test !all(sample_mask)
    end

    @testset "DAEDS" begin
        @test_throws ArgumentError DAEDS(-1, 0.1, 0.3)
        @test_throws ArgumentError DAEDS(0, 0.1, 0.3)
        @test_throws ArgumentError DAEDS(5, -1, 0.3)
        @test_throws ArgumentError DAEDS(5, 0.1, -1)
        s = DAEDS(5, 0.1, 0.3)
        sample_mask = sample(s, d, l)
        @test length(sample_mask) == length(l)
        @test !all(sample_mask)
    end

    @testset "DBSRSVDD" begin
        @test_throws ArgumentError DBSRSVDD(-1, 0.1)
        @test_throws ArgumentError DBSRSVDD(0, 0.1)
        @test_throws ArgumentError DBSRSVDD(5, -1)
        s = DBSRSVDD(5, 0.5)
        sample_mask = sample(s, d, l)
        @test length(sample_mask) == length(l)
        @test !all(sample_mask)
    end

    @testset "FBPE" begin
        @test_throws ArgumentError FBPE(-1)
        @test_throws ArgumentError FBPE(0)
        s = FBPE(2)
        sample_mask = sample(s, d, l)
        @test length(sample_mask) == length(l)
        @test !all(sample_mask)
    end

    @testset "HSC" begin
        @test_throws ArgumentError HSC(-1, 1, threshold_strat, c, init_strat, TEST_SOLVER)
        @test_throws ArgumentError HSC(0, 1, threshold_strat, c, init_strat, TEST_SOLVER)
        @test_throws ArgumentError HSC(1, -1, threshold_strat, c, init_strat, TEST_SOLVER)
        @test_throws ArgumentError HSC(1, 0, threshold_strat, c, init_strat, TEST_SOLVER)
        s = HSC(5, 2, threshold_strat, :test, init_strat, TEST_SOLVER)
        @test_throws ArgumentError sample(s, d, l)
        s = HSC(5, 2, threshold_strat, c, init_strat, TEST_SOLVER)
        sample_mask = sample(s, d, l)
        @test length(sample_mask) == length(l)
        @test !all(sample_mask)
    end

    @testset "HSR" begin
        @test_throws ArgumentError HSR(-1, 0.1)
        @test_throws ArgumentError HSR(0, 0.1)
        @test_throws ArgumentError HSR(5, -1)
        s = HSR(5, 0.5)
        sample_mask = sample(s, d, l)
        @test length(sample_mask) == length(l)
        @test !all(sample_mask)
    end

    @testset "IESRSVDD" begin
        @test_throws ArgumentError IESRSVDD(c, -0.1)
        @test_throws ArgumentError IESRSVDD(c, 1.1)
        s = IESRSVDD(:foo, 0.1)
        @test_throws ArgumentError sample(s, d, l)
        s = IESRSVDD(c, 0.5)
        sample_mask = sample(s, d, l)
        @test length(sample_mask) == length(l)
        @test !all(sample_mask)
    end

    @testset "KFNCBD" begin
        @test_throws ArgumentError KFNCBD(-1, 0.1)
        @test_throws ArgumentError KFNCBD(0, 0.1)
        @test_throws ArgumentError KFNCBD(5, -0.1)
        @test_throws ArgumentError KFNCBD(5, 1.5)
        s = KFNCBD(20, 0.5)
        sample_mask = sample(s, d, l)
        @test length(sample_mask) == length(l)
        @test !all(sample_mask)
    end

    @testset "KMSVDD" begin
        @test_throws ArgumentError KMSVDD(-1, init_strat, TEST_SOLVER)
        @test_throws ArgumentError KMSVDD(0, init_strat, TEST_SOLVER)
        s = KMSVDD(2, init_strat, TEST_SOLVER)
        sample_mask = sample(s, d, l)
        @test length(sample_mask) == length(l)
        @test !all(sample_mask)
    end

    @testset "NDPSR" begin
        @test_throws ArgumentError NDPSR(-1, 5)
        @test_throws ArgumentError NDPSR(0, 5)
        @test_throws ArgumentError NDPSR(5, -1)
        s = NDPSR(20, 15)
        sample_mask = sample(s, d, l)
        @test length(sample_mask) == length(l)
        @test !all(sample_mask)
    end

    @testset "RSSVDD" begin
        @test_throws ArgumentError RSSVDD(-1, 1, 1, 1, init_strat, TEST_SOLVER)
        @test_throws ArgumentError RSSVDD(0, 1, 1, 1, init_strat, TEST_SOLVER)
        @test_throws ArgumentError RSSVDD(1, -1, 1, 1, init_strat, TEST_SOLVER)
        @test_throws ArgumentError RSSVDD(1, 0, 1, 1, init_strat, TEST_SOLVER)
        @test_throws ArgumentError RSSVDD(1, 1, -1, 1, init_strat, TEST_SOLVER)
        @test_throws ArgumentError RSSVDD(1, 1, 0, 1, init_strat, TEST_SOLVER)
        @test_throws ArgumentError RSSVDD(1, 1, 1, -1, init_strat, TEST_SOLVER)
        @test_throws ArgumentError RSSVDD(1, 1, 1, 0, init_strat, TEST_SOLVER)
        s = RSSVDD(10, 2, init_strat, TEST_SOLVER)
        sample_mask = sample(s, d, l)
        @test length(sample_mask) == length(l)
        @test !all(sample_mask)
    end

    @testset "RAPID" begin
        s = RAPID(threshold_strat, :test)
        @test_throws ArgumentError sample(s, d, l)
        s = RAPID(threshold_strat, c)
        sample_mask = sample(s, d, l)
        @test length(sample_mask) == length(l)
        @test !all(sample_mask)
        s = RAPID(threshold_strat_gt, c)
        sample_mask = sample(s, d, l)
        @test length(sample_mask) == length(l)
        @test !all(sample_mask)
        _, inlier_mask, _ = split_masks(threshold_strat_gt, c, l)
        @test all(inlier_mask[sample_mask])
    end

    @testset "HDS" begin
        @test_throws ArgumentError HDS(1.5, c)
        @test_throws ArgumentError HDS(-1.0, c)
        s = HDS(1.0, :foo)
        @test_throws ArgumentError sample(s, d, l)
        s = HDS(1.0, c, nothing)
        sample_mask1 = sample(s, d, l)
        @test length(sample_mask1) == length(l)
        @test all(sample_mask1)
        s = HDS(1.0, c, threshold_strat)
        sample_mask2 = sample(s, d, l)
        @test length(sample_mask2) == length(l)
        @test !all(sample_mask2)
        @test count(sample_mask1) > count(sample_mask2)
        s = HDS(0.1, c, threshold_strat)
        sample_mask3 = sample(s, d, l)
        @test length(sample_mask3) == length(l)
        @test count(sample_mask2) > count(sample_mask3)
    end

    @testset "PreFilteringWrapper" begin
        sampler = RandomRatioSampler(1.0)
        s = PreFilteringWrapper(sampler, :test, threshold_strat)
        @test_throws ArgumentError sample(s, d, l)
        s1 = PreFilteringWrapper(sampler, c, nothing)
        sample_mask1 = sample(s1, d, l)
        @test length(sample_mask1) == length(l)
        @test all(sample_mask1)
        s2 = PreFilteringWrapper(sampler, c, threshold_strat)
        sample_mask2 = sample(s2, d, l)
        @test length(sample_mask2) == length(l)
        @test count(sample_mask1) > count(sample_mask2)
        sampler = RandomRatioSampler(0.5)
        s3 = PreFilteringWrapper(sampler, c, threshold_strat)
        sample_mask3 = sample(s3, d, l)
        @test length(sample_mask3) == length(l)
        @test count(sample_mask2) > count(sample_mask3)
    end
end


text:
using Flux, Trebuchet
using Zygote: forwarddiff
using Statistics: mean
using Random

lerp(x, lo, hi) = x*(hi-lo)+lo

function shoot(wind, angle, weight)
  Trebuchet.shoot((wind, Trebuchet.deg2rad(angle), weight))[2]
end

shoot(ps) = forwarddiff(p -> shoot(p...), ps)

Random.seed!(0)

model = Chain(Dense(2, 16, σ),
              Dense(16, 64, σ),
              Dense(64, 16, σ),
              Dense(16, 2)) |> f64

θ = params(model)

function aim(wind, target)
  angle, weight = model([wind, target])
  angle = σ(angle)*90
  weight = weight + 200
  angle, weight
end

distance(wind, target) =
  shoot(Tracker.collect([wind, aim(wind, target)...]))

function loss(wind, target)
  try
    (distance(wind, target) - target)^2
  catch e
    # Roots.jl sometimes give convergence errors, ignore them
    param(0)
  end
end

DIST  = (20, 100)	# Maximum target distance
SPEED =   5 # Maximum wind speed

target() = (randn() * SPEED, lerp(rand(), DIST...))

meanloss() = mean(sqrt(loss(target()...)) for i = 1:100)

opt = ADAM()

dataset = (target() for i = 1:100_000)
cb = Flux.throttle(() -> @show(meanloss()), 10)

Flux.train!(loss, θ, dataset, opt, cb = cb)

max_stars_repo_name: c-p-murphy/model-zoo

max_issues_repo_name: c-p-murphy/model-zoo

max_forks_repo_name: c-p-murphy/model-zoo


text:
__precompile__(false)

# @osx_only begin
#     const qtlibdir = "/Users/kfischer/Projects/qt-everywhere-opensource-src-5.3.1/qtbase/~/usr/lib/"
#     const QtCore = joinpath(qtlibdir,"QtCore.framework/")
#     const QtWidgets = joinpath(qtlibdir,"QtWidgets.framework/")

#     addHeaderDir(qtlibdir; isFramework = true, kind = C_System)

#     dlopen(joinpath(QtCore,"QtCore_debug"))
#     addHeaderDir(joinpath(QtCore,"Headers"), kind = C_System)
#     addHeaderDir("/Users/kfischer/Projects/qt-everywhere-opensource-src-5.3.1/qtbase/lib/QtCore.framework/Headers/5.3.1/QtCore")

#     cxxinclude("/Users/kfischer/Projects/qt-everywhere-opensource-src-5.3.1/qtbase/lib/QtCore.framework/Headers/5.3.1/QtCore/private/qcoreapplication_p.h")

#     dlopen(joinpath(QtWidgets,"QtWidgets"))
#     addHeaderDir(joinpath(QtWidgets,"Headers"), kind = C_System)
# end

using Cxx

@linux_only begin
    const qtincdir = "/usr/include/qt5"
    const qtlibdir = "/usr/lib64"
    const QtCore = joinpath(qtincdir, "QtCore/")
    const QtWidgets = joinpath(qtincdir, "QtWidgets/")

    addHeaderDir(qtincdir, kind = C_System)
    addHeaderDir(QtCore, kind = C_System)
    addHeaderDir(QtWidgets, kind = C_System)

    Libdl.dlopen(joinpath(qtlibdir,"libQt5Core.so"), Libdl.RTLD_GLOBAL)
    Libdl.dlopen(joinpath(qtlibdir,"libQt5Gui.so"), Libdl.RTLD_GLOBAL)
    Libdl.dlopen(joinpath(qtlibdir,"libQt5Widgets.so"), Libdl.RTLD_GLOBAL)
end

cxx"""
    #include <QtCore>
    #include <QApplication>
    #include <QMessageBox>
    #include <QPushButton>
"""

# # init the QApplication object
# function initApp()
#     # notes: This is pretty stupid, but it seems QApplication is capturing the pointer
#     # to the reference, so we can't just # stack allocate it because that won't
#     # be valid for exec
#     ac = [Int32(1)]
#     global const a = "julia"
#     x = Ptr{UInt8}[pointer(a),C_NULL]
#     app = @cxxnew QApplication(*(pointer(ac)),pointer(x))
# end
# const app = initApp()

const ac = [Int32(1)]
const a = "julia"
x = Ptr{UInt8}[pointer(a),C_NULL]
const app = @cxxnew QApplication(*(pointer(ac)),pointer(x))

# this will be called from the julia main loop
function update_loop(_timer)
    icxx"""
        $app->processEvents();
    """
end

function startEventLoop()
    # kick off the gui loop callbacks
    Base.Timer( update_loop, 0.1, 0.005 )
end

# BUG: this version doesn't work - unresponsive gui
#update_loop(_timer) = @cxx app.processEvents() # default QEventFlags::AllEvents

# add silly test
say_hi() = println("Hi!")::Void
cxx"""
#include <iostream>
void handle_hi()
{
    $:(say_hi());
}
"""

function createMessageBox()
    # create messagebox
    mb = @cxxnew QMessageBox(@cxx(QMessageBox::Information),
                          pointer("Hello World"),
                          pointer("This is a QMessageBox"))

    # add buttons
    @cxx mb->addButton(@cxx(QMessageBox::Ok))

    hibtn = @cxxnew QPushButton(pointer("Say Hi!"))
    @cxx mb->addButton(hibtn, @cxx(QMessageBox::ActionRole))

    setup(hibtn)

    # display the window
    @cxx mb->setWindowModality(@cxx(Qt::NonModal))
    @cxx mb->show()

    mb
end

# BUGS?
# - type translation doesn't work right for $:(btn) if I call connect from a cxx""" block
# - I get isexprs assertion failure if I try to use a lambda. think it might be a
#   block parsing issue though.
function setup(btn)
    icxx"""
        QObject::connect($btn, &QPushButton::clicked, handle_hi );
    """
end
# setup(hibtn)

# exit if not interactive shell
!isinteractive() && stop_timer(timer)


text:
using YaoCompiler
using Test

@testset "interpreter" begin
    include("interp.jl")
end

@testset "emulate" begin
    include("emulate.jl")
end

@testset "printing" begin
    include("printing.jl")
end

@testset "zx optimization" begin
    include("zx.jl")
end


text:
module StaticThreadSchedules

import Base.Threads: threadid, threading_run, nthreads, _threadsfor
export @threads

function _threadsfor(iter, lbody, schedule, cost_model)
    lidx = iter.args[1]         # index
    range = iter.args[2]

    quote
        local threadsfor_fun
        # compute the normalization of the cost
        per_thread_cost = 0.0
        for i in $(esc(range))
            per_thread_cost += Float64($(esc(cost_model))(i))
        end
        per_thread_cost = per_thread_cost / nthreads() + 10eps()
        # @show per_thread_cost
        chunk_indices = zeros(Int, nthreads()+1)
        current_chunk = 1
        chunk_indices[1] = 1
        accumulated_cost = 0.0
        for i in $(esc(range))
            # @show i, accumulated_cost
            accumulated_cost += Float64($(esc(cost_model))(i))
            if accumulated_cost > per_thread_cost
                accumulated_cost -= per_thread_cost
                current_chunk += 1
                chunk_indices[current_chunk] = i
            end
        end
        chunk_indices[end] = lastindex($(esc(range))) + 1
        # @show chunk_indices

        let range = $(esc(range)), chunks_ = chunk_indices
            function threadsfor_fun(onethread=false)
                r = range # Load into local variable
                chunks = chunks_
                tid = threadid()

                f = chunks[tid]
                l = chunks[tid+1]-1

                # run this thread's iterations
                for i = f:l
                    local $(esc(lidx)) = @inbounds r[i]
                    $(esc(lbody))
                end
            end
        end
        if threadid() != 1 || ccall(:jl_in_threaded_region, Cint, ()) != 0
            $(if schedule === :static
              :(error("`@threads :static` can only be used from thread 1 and not nested"))
              else
              # only use threads when called from thread 1, outside @threads
              :(Base.invokelatest(threadsfor_fun, true))
              end)
        else
            threading_run(threadsfor_fun)
        end
        nothing
    end
end

"""
    StaticThreadSchedules.@threads [schedule] [optional scaling] for ... end

A macro to parallelize a `for` loop to run with multiple threads. Splits the iteration
space among multiple tasks and runs those tasks on threads according to a scheduling
policy.
A barrier is placed at the end of the loop which waits for all tasks to finish
execution.

The `schedule` argument can be used to request a particular scheduling policy.
If one specifies `:static` and a cost function, then the work will be split such that the cost
is approximately equal per thread. Specifying `:static` is an error
if used from inside another `@threads` loop or from a thread other than 1.

The default schedule (used when no `schedule` argument is present) is subject to change.

!!! compat "Julia 1.5"
    The `schedule` argument is available as of Julia 1.5.
"""
macro threads(args...)
    na = length(args)
    if na == 3
        sched, cost_model, ex = args
        if sched isa QuoteNode
            sched = sched.value
        elseif sched isa Symbol
            # for now only allow quoted symbols
            sched = nothing
        end
        if sched !== :static
            throw(ArgumentError("unsupported schedule argument in @threads"))
        end
        return _threadsfor(ex.args[1], ex.args[2], sched, cost_model)
    elseif na == 2
        sched, ex = args
        if sched isa QuoteNode
            sched = sched.value
        elseif sched isa Symbol
            # for now only allow quoted symbols
            sched = nothing
        end
        if sched !== :static
            throw(ArgumentError("unsupported schedule argument in @threads"))
        end
    elseif na == 1
        sched = :default
        ex = args[1]
    else
        throw(ArgumentError("wrong number of arguments in @threads"))
    end
    if !(isa(ex, Expr) && ex.head === :for)
        throw(ArgumentError("@threads requires a `for` loop expression"))
    end
    if !(ex.args[1] isa Expr && ex.args[1].head === :(=))
        throw(ArgumentError("nested outer loops are not currently supported by @threads"))
    end
    return _threadsfor(ex.args[1], ex.args[2], sched)
end

end


text:
# Structs for stress and velocity fields (used for multiple dispatch)

abstract type Fields end

"""
Return a list of fields which contain `s`.
Retruns all the fields defined here by default.
"""
Fields(s = "") = filter(x -> contains(string(x), string(s)), nameof.(subtypes(Fields)))

# define pressure for acoustic
struct p <: Fields end
# define stress and velocity (tauxx, tauxy, tauxz,)
for f in vcat(dim_names(3, "v"), dim_names(3, "tau"; order = 2))
    @eval struct $f <: Fields end
end
# define spatial derivatives of stress, velocity, pressure
dnames = dim_names(3, "d") # derivatives
vnames = dim_names(3, "v") # velocity
taunames = dim_names(3, "tau"; order = 2)
for (i, x) in enumerate(dim_names(3))
    tn = filter(y -> contains(string(y), string(x)), taunames)
    for t in vcat(tn, vnames, [:p])
        @eval struct $(Symbol("d", string(t), "d", string(x))) <: Fields end
    end
end

#-------------------------------------------------------
# grid indices and linear-interpolation weights
#-------------------------------------------------------
"""
Output CartesianIndices and corresponding linear interpolation weights for point P in mgrid
For example, when 2D,
mgrid=[range(1.3, stop=10.6,step=0.003), range(1.2,stop=15.3,step=0.004)]
P=[5,5]
"""
function get_indices_weights(mgrid, P)
    @assert length(mgrid) == length(P)
    N = length(mgrid)
    idx = [Interpolation.indminn(mgrid[i], P[i], 2) for i = 1:N]
    denomI = inv(prod([diff(mgrid[i][idx[i]])[1] for i = 1:N]))
    c = CartesianIndices(Tuple(broadcast(x -> x[1]:x[2], idx)))
    weights = zeros(length(c))
    for (i, cc) in enumerate(c)
        weights[i] = prod([abs(mgrid[i][cc[i]] - P[i]) for i = 1:N]) * denomI
    end
    return c, weights
end

for dimnames in [zip([:1, :2, :3], [:z, :y, :x]), zip([:1, :2], [:z, :x])]
    grids = broadcast(x -> Symbol(string("m", x)), getindex.(collect(dimnames), 2))
    for (idim, dim) in dimnames
        v = Symbol("v", string(dim))
        g = Symbol("m", string(dim))
        @eval function get_indices_weights(::$v, $(grids...), P)
            # velocity grids are stagerred in respective dimensions, so generate half grids in respective dimensions
            $g = range(
                $g[1] - step($g) * (_fd.order - 1) * 0.5,
                step = step($g),
                length = length($g) + _fd.order - 1,
            )
            return get_indices_weights([$(grids...)], P)
        end
    end
    for f in [:p, :tauxx, :tauyy, :tauzz]
        @eval function get_indices_weights(::$f, $(grids...), P)
            # pressure/ stress grids unchanged
            return get_indices_weights([$(grids...)], P)
        end
    end
end

#-------------------------------------------------------
# Acoustic
#-------------------------------------------------------
Base.zeros(::p, ::FdtdAcou, nz, nx) = Data.Array(zeros(nz, nx))
Base.zeros(::p, ::FdtdAcou, nz, ny, nx) = Data.Array(zeros(nz, ny, nx))
Base.zeros(::dpdx, ::FdtdAcou, nz, nx; pml::Bool = false) =
    Data.Array(zeros(nz - 2(_fd.order - 1), pml ? 2 * npml : nx - 1(_fd.order - 1)))
Base.zeros(::dpdx, ::FdtdAcou, nz, ny, nx; pml::Bool = false) = Data.Array(
    zeros(
        nz - 2(_fd.order - 1),
        ny - 2(_fd.order - 1),
        pml ? 2 * npml : nx - 1(_fd.order - 1),
    ),
)
Base.zeros(::dpdz, ::FdtdAcou, nz, nx; pml::Bool = false) =
    Data.Array(zeros(pml ? 2 * npml : nz - 1(_fd.order - 1), nx - 2(_fd.order - 1)))
Base.zeros(::dpdz, ::FdtdAcou, nz, ny, nx; pml::Bool = false) = Data.Array(
    zeros(
        pml ? 2 * npml : nz - 1(_fd.order - 1),
        ny - 2(_fd.order - 1),
        nx - 2(_fd.order - 1),
    ),
)
Base.zeros(::dpdy, ::FdtdAcou, nz, ny, nx; pml::Bool = false) = Data.Array(
    zeros(
        nz - 2(_fd.order - 1),
        pml ? 2 * npml : ny - 1(_fd.order - 1),
        nx - 2(_fd.order - 1),
    ),
)
Base.zeros(::vx, ::FdtdAcou, nz, nx) = Data.Array(zeros(nz, nx + _fd.order - 1))
Base.zeros(::vz, ::FdtdAcou, nz, nx) = Data.Array(zeros(nz + _fd.order - 1, nx))
Base.zeros(::vx, ::FdtdAcou, nz, ny, nx) = Data.Array(zeros(nz, ny, nx + _fd.order - 1))
Base.zeros(::vz, ::FdtdAcou, nz, ny, nx) = Data.Array(zeros(nz + _fd.order - 1, ny, nx))
Base.zeros(::vy, ::FdtdAcou, nz, ny, nx) = Data.Array(zeros(nz, ny + 1, nx))
Base.zeros(::dvxdx, ::FdtdAcou, nz, nx; pml::Bool = false) =
    Data.Array(zeros(nz, pml ? 2 * npml : nx))
Base.zeros(::dvzdz, ::FdtdAcou, nz, nx; pml::Bool = false) =
    Data.Array(zeros(pml ? 2 * npml : nz, nx))
Base.zeros(::dvxdx, ::FdtdAcou, nz, ny, nx; pml::Bool = false) =
    Data.Array(zeros(nz, ny, pml ? 2 * npml : nx))
Base.zeros(::dvzdz, ::FdtdAcou, nz, ny, nx; pml::Bool = false) =
    Data.Array(zeros(pml ? 2 * npml : nz, ny, nx))
Base.zeros(::dvydy, ::FdtdAcou, nz, ny, nx; pml::Bool = false) =
    Data.Array(zeros(nz, pml ? 2 * npml : ny, nx))
Base.zeros(::dvxdz, ::FdtdAcou, nz, nx; pml::Bool = false) = Data.Array(zeros(1, 1))
Base.zeros(::dvzdx, ::FdtdAcou, nz, nx; pml::Bool = false) = Data.Array(zeros(1, 1))
Base.zeros(::dvxdz, ::FdtdAcou, nz, ny, nx; pml::Bool = false) = Data.Array(zeros(1, 1, 1))
Base.zeros(::dvydx, ::FdtdAcou, nz, ny, nx; pml::Bool = false) = Data.Array(zeros(1, 1, 1))
Base.zeros(::dvzdx, ::FdtdAcou, nz, ny, nx; pml::Bool = false) = Data.Array(zeros(1, 1, 1))
Base.zeros(::dvydz, ::FdtdAcou, nz, ny, nx; pml::Bool = false) = Data.Array(zeros(1, 1, 1))
Base.zeros(::dvxdy, ::FdtdAcou, nz, ny, nx; pml::Bool = false) = Data.Array(zeros(1, 1, 1))
Base.zeros(::dvzdy, ::FdtdAcou, nz, ny, nx; pml::Bool = false) = Data.Array(zeros(1, 1, 1))

#-------------------------------------------------------
# Elastic
#-------------------------------------------------------
Base.zeros(::tauxx, ::FdtdElastic, nz, ny, nx) = Data.Array(zeros(nz, ny, nx))
Base.zeros(::tauyy, ::FdtdElastic, nz, ny, nx) = Data.Array(zeros(nz, ny, nx))
Base.zeros(::tauzz, ::FdtdElastic, nz, ny, nx) = Data.Array(zeros(nz, ny, nx))
Base.zeros(::tauxy, ::FdtdElastic, nz, ny, nx) =
    Data.Array(zeros(nz - 2(_fd.order - 1), ny - 1(_fd.order - 1), nx - 1(_fd.order - 1)))
Base.zeros(::tauxz, ::FdtdElastic, nz, ny, nx) =
    Data.Array(zeros(nz - 1(_fd.order - 1), ny - 2(_fd.order - 1), nx - 1(_fd.order - 1)))
Base.zeros(::tauyz, ::FdtdElastic, nz, ny, nx) =
    Data.Array(zeros(nz - 1(_fd.order - 1), ny - 1(_fd.order - 1), nx - 2(_fd.order - 1)))

Base.zeros(::tauxx, ::FdtdElastic, nz, nx) = Data.Array(zeros(nz, nx))
Base.zeros(::tauzz, ::FdtdElastic, nz, nx) = Data.Array(zeros(nz, nx))
Base.zeros(::tauxz, ::FdtdElastic, nz, nx) =
    Data.Array(zeros(nz - 1(_fd.order - 1), nx - 1(_fd.order - 1)))

Base.zeros(::vx, ::FdtdElastic, nz, ny, nx) = Data.Array(zeros(nz, ny, nx + _fd.order - 1))
Base.zeros(::vy, ::FdtdElastic, nz, ny, nx) = Data.Array(zeros(nz, ny + _fd.order - 1, nx))
Base.zeros(::vz, ::FdtdElastic, nz, ny, nx) = Data.Array(zeros(nz + _fd.order - 1, ny, nx))

Base.zeros(::vx, ::FdtdElastic, nz, nx) = Data.Array(zeros(nz, nx + _fd.order - 1))
Base.zeros(::vz, ::FdtdElastic, nz, nx) = Data.Array(zeros(nz + _fd.order - 1, nx))

Base.zeros(::dvxdx, ::FdtdElastic, nz, ny, nx; pml::Bool = false) =
    Data.Array(zeros(nz, ny, pml ? 2 * npml : nx))
Base.zeros(::dvydy, ::FdtdElastic, nz, ny, nx; pml::Bool = false) =
    Data.Array(zeros(nz, pml ? 2 * npml : ny, nx))
Base.zeros(::dvzdz, ::FdtdElastic, nz, ny, nx; pml::Bool = false) =
    Data.Array(zeros(pml ? 2 * npml : nz, ny, nx))
Base.zeros(::dvxdy, ::FdtdElastic, nz, ny, nx; pml::Bool = false) = Data.Array(
    zeros(
        nz - 2(_fd.order - 1),
        pml ? 2 * npml : ny - 1(_fd.order - 1),
        nx - 1(_fd.order - 1),
    ),
)
Base.zeros(::dvxdz, ::FdtdElastic, nz, ny, nx; pml::Bool = false) = Data.Array(
    zeros(
        pml ? 2 * npml : nz - 1(_fd.order - 1),
        ny - 2(_fd.order - 1),
        nx - 1(_fd.order - 1),
    ),
)
Base.zeros(::dvydx, ::FdtdElastic, nz, ny, nx; pml::Bool = false) = Data.Array(
    zeros(
        nz - 2(_fd.order - 1),
        ny - 1(_fd.order - 1),
        pml ? 2 * npml : nx - 1(_fd.order - 1),
    ),
)
Base.zeros(::dvydz, ::FdtdElastic, nz, ny, nx; pml::Bool = false) = Data.Array(
    zeros(
        pml ? 2 * npml : nz - 1(_fd.order - 1),
        ny - 1(_fd.order - 1),
        nx - 2(_fd.order - 1),
    ),
)
Base.zeros(::dvzdx, ::FdtdElastic, nz, ny, nx; pml::Bool = false) = Data.Array(
    zeros(
        nz - 1(_fd.order - 1),
        ny - 2(_fd.order - 1),
        pml ? 2 * npml : nx - 1(_fd.order - 1),
    ),
)
Base.zeros(::dvzdy, ::FdtdElastic, nz, ny, nx; pml::Bool = false) = Data.Array(
    zeros(
        nz - 1(_fd.order - 1),
        pml ? 2 * npml : ny - 1(_fd.order - 1),
        nx - 2(_fd.order - 1),
    ),
)

Base.zeros(::dvxdx, ::FdtdElastic, nz, nx; pml::Bool = false) =
    Data.Array(zeros(nz, pml ? 2 * npml : nx))
Base.zeros(::dvzdz, ::FdtdElastic, nz, nx; pml::Bool = false) =
    Data.Array(zeros(pml ? 2 * npml : nz, nx))
Base.zeros(::dvxdz, ::FdtdElastic, nz, nx; pml::Bool = false) =
    Data.Array(zeros(pml ? 2 * npml : nz - 1(_fd.order - 1), nx - 1(_fd.order - 1)))
Base.zeros(::dvzdx, ::FdtdElastic, nz, nx; pml::Bool = false) =
    Data.Array(zeros(nz - 1(_fd.order - 1), pml ? 2 * npml : nx - 1(_fd.order - 1)))

Base.zeros(::dtauxxdx, ::FdtdElastic, nz, ny, nx; pml::Bool = false) = Data.Array(
    zeros(
        nz - 2(_fd.order - 1),
        ny - 2(_fd.order - 1),
        pml ? 2 * npml : nx - 1(_fd.order - 1),
    ),
)
Base.zeros(::dtauyydy, ::FdtdElastic, nz, ny, nx; pml::Bool = false) = Data.Array(
    zeros(
        nz - 2(_fd.order - 1),
        pml ? 2 * npml : ny - 1(_fd.order - 1),
        nx - 2(_fd.order - 1),
    ),
)
Base.zeros(::dtauzzdz, ::FdtdElastic, nz, ny, nx; pml::Bool = false) = Data.Array(
    zeros(
        pml ? 2 * npml : nz - 1(_fd.order - 1),
        ny - 2(_fd.order - 1),
        nx - 2(_fd.order - 1),
    ),
)
Base.zeros(::dtauxydx, ::FdtdElastic, nz, ny, nx; pml::Bool = false) = Data.Array(
    zeros(
        nz - 2(_fd.order - 1),
        ny - 1(_fd.order - 1),
        pml ? 2 * npml : nx - 2(_fd.order - 1),
    ),
)
Base.zeros(::dtauxydy, ::FdtdElastic, nz, ny, nx; pml::Bool = false) = Data.Array(
    zeros(
        nz - 2(_fd.order - 1),
        pml ? 2 * npml : ny - 2(_fd.order - 1),
        nx - 1(_fd.order - 1),
    ),
)
Base.zeros(::dtauxzdx, ::FdtdElastic, nz, ny, nx; pml::Bool = false) = Data.Array(
    zeros(
        nz - 1(_fd.order - 1),
        ny - 2(_fd.order - 1),
        pml ? 2 * npml : nx - 2(_fd.order - 1),
    ),
)
Base.zeros(::dtauxzdz, ::FdtdElastic, nz, ny, nx; pml::Bool = false) = Data.Array(
    zeros(
        pml ? 2 * npml : nz - 2(_fd.order - 1),
        ny - 2(_fd.order - 1),
        nx - 1(_fd.order - 1),
    ),
)
Base.zeros(::dtauyzdy, ::FdtdElastic, nz, ny, nx; pml::Bool = false) = Data.Array(
    zeros(
        nz - 1(_fd.order - 1),
        pml ? 2 * npml : ny - 2(_fd.order - 1),
        nx - 2(_fd.order - 1),
    ),
)
Base.zeros(::dtauyzdz, ::FdtdElastic, nz, ny, nx; pml::Bool = false) = Data.Array(
    zeros(
        pml ? 2 * npml : nz - 2(_fd.order - 1),
        ny - 1(_fd.order - 1),
        nx - 2(_fd.order - 1),
    ),
)

Base.zeros(::dtauxxdx, ::FdtdElastic, nz, nx; pml::Bool = false) =
    Data.Array(zeros(nz - 2(_fd.order - 1), pml ? 2 * npml : nx - 1(_fd.order - 1)))
Base.zeros(::dtauzzdz, ::FdtdElastic, nz, nx; pml::Bool = false) =
    Data.Array(zeros(pml ? 2 * npml : nz - 1(_fd.order - 1), nx - 2(_fd.order - 1)))
Base.zeros(::dtauxzdx, ::FdtdElastic, nz, nx; pml::Bool = false) =
    Data.Array(zeros(nz - 1(_fd.order - 1), pml ? 2 * npml : nx - 2(_fd.order - 1)))
Base.zeros(::dtauxzdz, ::FdtdElastic, nz, nx; pml::Bool = false) =
    Data.Array(zeros(pml ? 2 * npml : nz - 2(_fd.order - 1), nx - 1(_fd.order - 1)))

# Dont need pressure for FdtdElastic, so no initialization
for f in Fields("p")
    @eval Base.zeros(::$f, ::FdtdElastic, args1...; args2...) =
        Data.Array(zeros(fill(1, length(args1))...))
end
# remove stress from acoustic simulations
for f in Fields("tau")
    @eval Base.zeros(::$f, ::FdtdAcou, args1...; args2...) =
        Data.Array(zeros(fill(1, length(args1))...))
end
# no "y" for 2D simulations
for f in Fields("y")
    @eval Base.zeros(::$f, attrib_mod, nz, nx; args...) = Data.Array(zeros(fill(1, 2)...))
end


text:
# * Ponderomotive potential
@doc raw"""
    ponderomotive_potential(f)

Return the [ponderomotive
potential](https://en.wikipedia.org/wiki/Ponderomotive_energy)
``U_p``, which is the cycle-average quiver energy of a free electron
in an electromagnetic field `f`. It is given by

```math
U_p =
\frac{e^2E_0^2}{4m\omega^2}=\frac{2e^2}{c\varepsilon_0m}\times\frac{I}{4\omega^2},
```

or, in atomic units,

```math
U_p = \frac{I}{4\omega^2}.
```
"""
ponderomotive_potential(f::Union{LinearField,TransverseField,WrappedField}) =
    params(f)[:Uₚ]

ponderomotive_potential(f::SumField) =
    ponderomotive_potential(f.a) + ponderomotive_potential(f.b)

# * Keldysh parameter

@doc raw"""
    keldysh(f, Iₚ)

The [Keldysh
parameter](https://en.wikipedia.org/wiki/Tunnel_ionization) relates
the strength of a dynamic electric field to that of the binding
potential of an atom. It is given by

```math
\gamma = \sqrt{\frac{I_p}{2U_p}},
```

where ``I_p`` is the ionization potential of the atom and ``U_p`` is
the ponderomotive potential of the dynamic field.
"""
keldysh(f::AbstractField, Iₚ::Unitful.Energy) =
    √(Iₚ/2ponderomotive_potential(f)) |> NoUnits

# * Free oscillation amplitude

@doc raw"""
    free_oscillation_amplitude(F)

Compute the free oscillation amplitude of an electric field `F`,
i.e. the mean excursion length during one cycle of the field, defined
as

```math
\alpha \defd \frac{F}{\omega^2}
```

where `F` is the peak amplitude, i.e. this is defined for one cycle of
a monochrome field.
"""
free_oscillation_amplitude(F::AbstractField) =
    amplitude(F)/photon_energy(F)^2

free_oscillation_amplitude(F::Union{WrappedField,NegatedField,DelayedField}) =
    free_oscillation_amplitude(parent(F))

# This is just a heuristic, i.e. if `F.a` and `F.b` are of the same
# frequency, but exactly out of phase, the free oscillation amplitude
# is actually zero.
free_oscillation_amplitude(F::SumField) =
    free_oscillation_amplitude(F.a) +
    free_oscillation_amplitude(F.b)

# * Print info

function show_strong_field_properties(io::IO, F::Union{LinearField,TransverseField})
    Uₚ = austrip(ponderomotive_potential(F))
    α = austrip(free_oscillation_amplitude(F))
    printfmt(io, "Uₚ = {1:.4f} Ha = {2:s} => α = {3:.4f} Bohr = {4:s}",
             Uₚ, au2si_round(Uₚ, u"eV"),
             α, au2si_round(α, u"nm"))
end

# * Exports

export ponderomotive_potential, keldysh, free_oscillation_amplitude


text:
module WhittleLikelihoodInference

using FFTW, LinearAlgebra, ToeplitzMatrices, StaticArrays, RecipesBase
using Distributions, SpecialFunctions, LazyArrays
using ChainRulesCore
using Optim, DSP

import Base: ndims, show, size, getindex, @propagate_inbounds
import StaticArrays: triangularnumber

export
    TimeSeriesModel,
    UnknownAcvTimeSeriesModel,
    AdditiveTimeSeriesModel,
    npars,
    ndims,
    parameternames,
    lowerbounds,
    upperbounds,
    parameter,
    sdf,
    asdf,
    acv,
    EI,
    coherancy,
    coherance,
    groupdelay,
    simulate_gp,
    fit,
    ## non-parametrics
    Periodogram,
    BartlettPeriodogram,
    SpectralEstimate,
    CoherancyEstimate,
    ## models
    OU,
    CorrelatedOU,
    Matern,
    Matern1D,
    Matern2D,
    Matern3D,
    Matern4D,
    ## Whittle
    DebiasedWhittleLikelihood,
    WhittleLikelihood

    include("typestructure.jl")
    include("memoryallocation.jl")
    include("properties/spectraldensity.jl")
    include("properties/autocovariance.jl")
    include("properties/expectedperiodogram.jl")
    include("properties/coherancy.jl")

    include("simulation.jl")
    include("nonparametric.jl")
    include("plotting.jl")

    include("whittle/whittledata.jl")
    include("whittle/generalwhittle.jl")
    include("whittle/standardwhittle.jl")
    include("whittle/debiasedwhittle.jl")
    include("whittle/chainrules.jl")
    include("fit.jl")

    include("models/OU.jl")
    include("models/CorrelatedOU.jl")
    include("models/OUUnknown.jl")
    include("models/CorrelatedOUUnknown.jl")
    include("models/Matern.jl")
    include("models/MaternAcv.jl")

end


text:
include("sim_perf.jl")
include("cumsim_perf.jl")


text:
# AUTO GENERATED FILE - DO NOT EDIT

export title

"""
    title(;kwargs...)
    title(children::Any;kwargs...)
    title(children_maker::Function;kwargs...)

A Title component.
h1-h6 headings. For more information, see: https://mantine.dev/core/title/
Keyword arguments:
- `children` (String; optional): Primary content
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `align` (a value equal to: "left", "right", "center"; optional): Sets text-align css property
- `className` (String; optional): Often used with CSS to style elements with common properties
- `order` (Real; optional): Defines component and styles which will be used
- `style` (Dict; optional): Inline style override
"""
function title(; kwargs...)
        available_props = Symbol[:children, :id, :align, :className, :order, :style]
        wild_props = Symbol[]
        return Component("title", "Title", "dash_mantine_components", available_props, wild_props; kwargs...)
end

title(children::Any; kwargs...) = title(;kwargs..., children = children)
title(children_maker::Function; kwargs...) = title(children_maker(); kwargs...)