text:
## ExoplanetsSysSim/src/transit_detection_model.jl
## (c) 2015 Eric B. Ford

# Several functions below based on https://github.com/christopherburke/KeplerPORTs/blob/master/KeplerPORTs_utils.py
# That follows the procedure outlined in Burke et al.(2015).
# However we don't currently interpolate the CDDP or mesthreshold to the relevant duration

function real_log_choose(m::Real, n::Real)
  lgamma(m+1)-lgamma(n+1)-lgamma(m-n+1.0)
end

function real_binom(k::Real, BigM::Real, f::Real)
    F1 = real_log_choose(BigM,k)
    F2 = k*log(f)
    F3 = (BigM-k)*log(1.0-f)
    x = exp(F1+F2+F3)
    return x
end

function kepler_window_function(num_transits_no_gaps::Real, duty_cycle::Real; min_transits::Real = 3.0)
  if num_transits_no_gaps < min_transits
     return 0.0
  else
     return max(1.0 - real_binom(min_transits,num_transits_no_gaps,duty_cycle), 0.0)
  end
end

# See Christiansen et al. (2015)  This assumes a linear limbdarkening coefficient of 0.6
function frac_depth_to_tps_depth(frac_depth::Real)
    const alp = 1.0874
    const bet = 1.0187
    const REALDEPTH2TPSSQUARE = 1.0  # WARNING: Waiting for this to be confirmed
    k = sqrt(frac_depth)
    tps_depth = min( (alp-bet*k) * frac_depth* REALDEPTH2TPSSQUARE, 1.0)   # NOTE: I added the max based on common sense
    return tps_depth
end

function detection_efficiency_theory(mes::Real; min_pdet_nonzero::Float64 = 0.0)
   const muoffset =  0.0
   const sig =  1.0
   const mesthresh = 7.1
   if mes > (9.0 - mesthresh - muoffset)
      return 1.0
   else
      pdet = 0.5 + 0.5*erf((mes - mesthresh - muoffset) / sqrt(2.0*sig*sig))
      pdet = pdet >= min_pdet_nonzero ? pdet : 0.0
      return pdet
   end
end

function detection_efficiency_fressin2013(mes::Real)
   const mesmin =  6.0
   const mesmax =  16.0
   if mes <= mesmin
      return 0.0
   elseif mes >= mesmax
      return 1.0
   else
     return (mes - mesmin) / (mesmax - mesmin)
  end
end

function detection_efficiency_christiansen2015(mes::Real; mes_threshold::Real = 7.1, min_pdet_nonzero::Float64 = 0.0)
   const a =  4.65  # from code for detection_efficiency(...) at https://github.com/christopherburke/KeplerPORTs/blob/master/KeplerPORTs_utils.py
   const b =  0.98
   #const a =  4.35 # from arxiv abstract. TODO DETAIL: Figure out which is better.  Informal testing showed it didn't matter
   #const b =  1.05
   usemes = max(0.0,mes - 7.1 - (mes_threshold - 7.1))
   pdet = cdf(Gamma(a,b), usemes)
   pdet = pdet >= min_pdet_nonzero ? pdet : 0.0
   return pdet
end

function detection_efficiency_dr25_simple(mes::Real)
   const a =  30.87  # from pg 16 of https://exoplanetarchive.ipac.caltech.edu/docs/KSCI-19110-001.pdf
   const b =  0.271
   const c = 0.940
   pdet = c*cdf(Gamma(a,b), mes)
   return pdet
end

# Resume code original to SysSim

function calc_snr_if_transit(t::KeplerTarget, depth::Real, duration::Real, sim_param::SimParam; num_transit::Real = 1)
  cdpp = t.cdpp[1,1]                                                            # TODO: Make CDPP lookup based on season/quarter/month and timescale.  Is CDPP timescale IMPORTANT or DETAIL?
  depth_tps = frac_depth_to_tps_depth(depth)                                        # WARNING: Hardcoded this conversion
  snr = depth_tps*sqrt(num_transit*duration*LC_rate)/cdpp              # WARNING: Assumes measurement uncertainties are uncorrelated & CDPP based on LC
end

function calc_snr_if_transit(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
  depth = calc_transit_depth(t,s,p)
  duration = calc_transit_duration(t,s,p)
  num_transit = calc_expected_num_transits(t,s,p,sim_param)
  calc_snr_if_transit(t,depth,duration,sim_param,num_transit=num_transit)
end

function calc_prob_detect_if_transit(t::KeplerTarget, snr::Real, sim_param::SimParam; num_transit::Real = 1)
  const min_transits = 3.0                                                    # WARNING: Hard coded 3 transit minimum
  const mes_threshold = 7.1                                                   # WARNING: Assuming 7.1 for all stars, durations
  const min_pdet_nonzero = 0.0                                                # TODO OPT: Figure out how to prevent a plethora of planets that are very unlikely to be detected due to using 0.0
  wf = kepler_window_function(num_transit, t.duty_cycle, min_transits=min_transits)   # TODO SCI DETAIL: Replace statistical model with checking actual transit times for long period planets
  return wf*detection_efficiency_christiansen2015(snr, mes_threshold=mes_threshold, min_pdet_nonzero=min_pdet_nonzero)        # WARNING: I've hardcoded this choice, since this is called so frequently
  #return wf*detection_efficiency_theory(snr, min_pdet_nonzero=min_pdet_nonzero)	     # WARNING: Hardcoded version of ideal detection efficiency
end

function calc_prob_detect_if_transit(t::KeplerTarget, depth::Real, duration::Real, sim_param::SimParam; num_transit::Real = 1)
  snr = calc_snr_if_transit(t,depth,duration,sim_param, num_transit=num_transit)
  return calc_prob_detect_if_transit(t, snr, sim_param, num_transit=num_transit)
end

function calc_prob_detect_if_transit(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
  depth = calc_transit_depth(t,s,p)
  duration = calc_transit_duration(t,s,p)
  ntr = calc_expected_num_transits(t,s,p,sim_param)
  calc_prob_detect_if_transit(t,depth,duration, sim_param, num_transit=ntr)
end

# Compute probability of detection if we average over impact parameters b~U[0,1)
function calc_ave_prob_detect_if_transit(t::KeplerTarget, snr_central::Real, sim_param::SimParam; num_transit::Real = 1)
  const min_transits = 3.0                                                    # WARNING: Hard coded 3 transit minimum
  const mes_threshold = 7.1                                                   # WARNING: Assuming 7.1 for all stars, durations
  const min_pdet_nonzero = 0.0                                                # TODO OPT: Figure out how to prevent a plethora of planets that are very unlikely to be detected due to using 0.0
  wf = kepler_window_function(num_transit, t.duty_cycle, min_transits=min_transits)
  num_impact_param = 5
  ave_detection_efficiency = 0.5*detection_efficiency_christiansen2015(snr_central, mes_threshold=mes_threshold, min_pdet_nonzero=min_pdet_nonzero)    # WARNING: Hardcoded Gamma CDF detection efficiency
  #ave_detection_efficiency = 0.5*detection_efficiency_theory(snr_central, min_pdet_nonzero=min_pdet_nonzero)	  # WARNING: Hardcoded ideal detection efficiency
  for i in 1:(num_impact_param-1)
    b = i/num_impact_param
    snr = snr_central*sqrt(sqrt((1-b)*(1+b)))
    ave_detection_efficiency += detection_efficiency_christiansen2015(snr, mes_threshold=mes_threshold, min_pdet_nonzero=min_pdet_nonzero)   # WARNING: Hardcoded Gamma CDF detection efficiency
    #ave_detection_efficiency += detection_efficiency_theory(snr, min_pdet_nonzero=min_pdet_nonzero)   # WARNING: Hardcoded ideal detection efficiency
  end
  ave_detection_efficiency /= num_impact_param
  return wf*ave_detection_efficiency
end

function calc_ave_prob_detect_if_transit(t::KeplerTarget, depth::Real, duration_central::Real, sim_param::SimParam; num_transit::Real = 1)
  snr_central = calc_snr_if_transit(t,depth,duration_central,sim_param, num_transit=num_transit)
  return calc_ave_prob_detect_if_transit(t, snr_central, sim_param, num_transit=num_transit)
end

function calc_ave_prob_detect_if_transit(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
  depth = calc_transit_depth(t,s,p)
  duration_central = calc_transit_duration_central(t,s,p)
  ntr = calc_expected_num_transits(t,s,p,sim_param)
  calc_ave_prob_detect_if_transit(t,depth,duration_central, sim_param, num_transit=ntr)
end


text:
module TinyInt
    # Overloaded functions
    import Base: getindex, length
    # Ensure library is compiled
    if Sys.iswindows()
        error("Windows is unsupported at this time.")
    end
    build_dir = string(dirname(@__FILE__), "/../deps/simdcomp")
    run(`make -s -C $build_dir`)
    libsimdcomp = build_dir * "/libsimdcomp.so"
    # Global C function tuple constants
    global const simdmaxmin = (:simdmaxmin, libsimdcomp)
    global const simdbits = (:bits, libsimdcomp)
    global const simdpackFOR = (:simdpackFOR, libsimdcomp)
    global const simdunpackFOR = (:simdunpackFOR, libsimdcomp)
    global const simdselectFOR = (:simdselectFOR, libsimdcomp)
    # Type definition and overloaded functions
    include("tinyintvec.jl")
    # Functions to pack and unpack arrays
    include("compress.jl")

    export TinyIntVec, pack, unpack, unpack!
end

ext: jl


text:
###############
# API methods #
###############

const AllowedDerivativeConfig{F,H} = Union{DerivativeConfig{Tag{F,H}}, DerivativeConfig{Tag{Void,H}}}

derivative(f!, y, x, cfg::DerivativeConfig) = throw(ConfigHismatchError(f, cfg))
derivative!(out, f!, y, x, cfg::DerivativeConfig) = throw(ConfigHismatchError(f, cfg))

"""
    ForwardDiff.derivative(f, x::Real)

Return `df/dx` evaluated at `x`, assuming `f` is called as `f(x)`.

This method assumes that `isa(f(x), Union{Real,AbstractArray})`.
"""
@inline function derivative(f::F, x::R) where {F,R<:Real}
    T = typeof(Tag(F, R))
    return extract_derivative(f(Dual{T}(x, one(x))))
end

"""
    ForwardDiff.derivative(f!, y::AbstractArray, x::Real, cfg::DerivativeConfig = DerivativeConfig(f!, y, x))

Return `df!/dx` evaluated at `x`, assuming `f!` is called as `f!(y, x)` where the result is
stored in `y`.
"""
@inline function derivative(f!::F, y::AbstractArray, x::R,
                            cfg::AllowedDerivativeConfig{F,H} = DerivativeConfig(f!, y, x)) where {F,R<:Real,H}
    ydual = cfg.duals
    seed!(ydual, y)
    f!(ydual, Dual{Tag{F,H}}(x, one(x)))
    map!(value, y, ydual)
    return extract_derivative(ydual)
end

"""
    ForwardDiff.derivative!(result::Union{AbstractArray,DiffResult}, f, x::Real)

Compute `df/dx` evaluated at `x` and store the result(s) in `result`, assuming `f` is called
as `f(x)`.

This method assumes that `isa(f(x), Union{Real,AbstractArray})`.
"""
@inline function derivative!(result::Union{AbstractArray,DiffResult},
                             f::F, x::R) where {F,R<:Real}
    T = typeof(Tag(F, R))
    ydual = f(Dual{T}(x, one(x)))
    extract_value!(result, ydual)
    extract_derivative!(result, ydual)
    return result
end

"""
    ForwardDiff.derivative!(result::Union{AbstractArray,DiffResult}, f!, y::AbstractArray, x::Real, cfg::DerivativeConfig = DerivativeConfig(f!, y, x))

Compute `df!/dx` evaluated at `x` and store the result(s) in `result`, assuming `f!` is
called as `f!(y, x)` where the result is stored in `y`.
"""
@inline function derivative!(result::Union{AbstractArray,DiffResult},
                             f!::F, y::AbstractArray, x::R,
                             cfg::AllowedDerivativeConfig{F,H} = DerivativeConfig(f!, y, x)) where {F,R<:Real,H}
    ydual = cfg.duals
    seed!(ydual, y)
    f!(ydual, Dual{Tag{F,H}}(x, one(x)))
    extract_value!(result, y, ydual)
    extract_derivative!(result, ydual)
    return result
end

#####################
# result extraction #
#####################

# non-mutating #
#--------------#

@inline extract_derivative(y::Dual{T,V,1}) where {T,V} = partials(y, 1)
@inline extract_derivative(y::Real) = zero(y)
@inline extract_derivative(y::AbstractArray) = map(extract_derivative, y)

# mutating #
#----------#

extract_derivative!(result::AbstractArray, y::AbstractArray) = map!(extract_derivative, result, y)
extract_derivative!(result::DiffResult, y) = DiffBase.derivative!(extract_derivative, result, y)


text:
module Arrow

using Mmap
import Dates
using Tables, SentinelArrays

using Base: @propagate_inbounds
import Base: ==

const ALIGNMENT = 8
const FILE_FORMAT_MAGIC_BYTES = b"ARROW1"
const CONTINUATION_INDICATOR_BYTES = 0xffffffff

# vendored flatbuffers code for now
include("FlatBuffers/FlatBuffers.jl")
using .FlatBuffers

include("metadata/Flatbuf.jl")
using .Flatbuf; const Meta = Flatbuf

include("utils.jl")
include("eltypes.jl")
include("arraytypes.jl")
include("table.jl")
include("write.jl")

end  # module Arrow


text:
export  Turtle,
        Forward,
        Turn,
        Pendown,
        Penup,
        Circle,
        Orientation,
        Rectangle,
        Push,
        Pop,
        Message,
        HueShift,
        Randomize_saturation,
        Pencolor,
        Penwidth,
        Pen_opacity_random,
        Reposition

"""
With a Turtle you can command a turtle to move and draw: turtle graphics.

The functions that start with a capital letter are: Forward, Turn,
Circle, Orientation, Rectangle, Pendown, Penup, Pencolor, Penwidth, and Reposition.

There are also some other functions. To see how they might be used, see Lindenmayer.jl.
"""

type Turtle
    xpos::Float64
    ypos::Float64
    pendown::Bool
    orientation::Float64 # stored in radians
    pencolor::Tuple{Float64,Float64,Float64} # it's an RGB turtle
    function Turtle(xpos, ypos, pendown, orientation, pencolor)
        new(xpos, ypos, pendown, orientation, pencolor)
    end
end

# a stack to hold turtle positions
const queue = Array{Array{Float64,1},1}()

"""
Move the turtle forward by `d` units. The stored position is updated.

    Forward(t::Turtle, d)
"""
function Forward(t::Turtle, d)
    move(t.xpos, t.ypos)
    t.xpos = t.xpos + (d * cos(t.orientation))
    t.ypos = t.ypos + (d * sin(t.orientation))
    if t.pendown
        gsave()
        line(t.xpos, t.ypos)
        sethue(t.pencolor...)
        stroke()
        grestore()
    else
        move(t.xpos, t.ypos)
    end
end

"""
Increase the turtle's rotation by `r` radians. See also `Orientation`.

    Turn(t::Turtle, r)
"""
function Turn(t::Turtle, r)
    t.orientation = mod2pi(t.orientation + deg2rad(r))
end
"""
Set the turtle's orientation to `r` radians. See also `Turn`.

    Orientation(t::Turtle, r)
"""
function Orientation(t::Turtle, r)
    t.orientation = mod2pi(deg2rad(r))
end
"""
Put that pen down and start drawing.

    Pendown(t::Turtle)
"""
function Pendown(t::Turtle)
    t.pendown = true
end
"""
Pick that pen up and stop drawing.

    Penup(t::Turtle)
"""
function Penup(t::Turtle)
    t.pendown = false
end
"""
Draw a filled circle centred at the current position with the given radius.

    Circle(t::Turtle, radius)
"""
function Circle(t::Turtle, radius)
    gsave()
    sethue(t.pencolor...)
    circle(t.xpos, t.ypos, radius, :fill)
    grestore()
end
"""
Draw a filled rectangle centred at the current position with the given radius.

    Rectangle(t::Turtle, width, height)

"""
function Rectangle(t::Turtle, width, height)
    gsave()
    sethue(t.pencolor...)
    rect(t.xpos-width/2, t.ypos-height/2, width, height, :fill)
    grestore()
end

"""
Save the turtle's position and orientation on a stack.
"""
function Push(t::Turtle)
# save xpos, ypos, turn in a queue
    global queue
    push!(queue, [t.xpos, t.ypos, t.orientation])
end

"""
Lift the turtle's position and orientation off a stack.
"""
function Pop(t::Turtle)
# restore xpos, ypos, turn from queue
    global queue
    if isempty(queue)
        warn("turtle can't pop")
        t.xpos, t.ypos, t.orientation = 0.0, 0.0, 0.0
    else
        t.xpos, t.ypos, t.orientation = pop!(queue)
    end
end

"""
Write some text at the current position.

    Message(t::Turtle, txt)
"""
function Message(t::Turtle, txt)
    gsave()
    sethue(t.pencolor...)
    text(txt, t.xpos, t.ypos)
    grestore()
end

"""
Shift the Hue of the turtle's pen forward by inc

    HueShift(t::Turtle, inc = 1)
"""
function HueShift(t::Turtle, inc = 1)
    r, g, b = t.pencolor
    hsv = convert(Colors.HSV, Colors.RGB(r, g, b))
    newhsv = ((hsv.h + inc) % 360, hsv.s, hsv.v)
    newrgb= convert(Colors.RGB, Colors.HSV(newhsv[1], newhsv[2], newhsv[3]))
    sethue(newrgb)
    t.pencolor = (newrgb.r, newrgb.g, newrgb.b)
end

"""
Randomize saturation of the turtle's pen color.

    Randomize_saturation(t::Turtle)
"""
function Randomize_saturation(t::Turtle)
    r, g, b = t.pencolor
    hsv = convert(Colors.HSV, Colors.RGB(r, g, b))
    newhsv = (hsv.h, rand(0.5:0.05:1.0), hsv.v)
    newrgb= convert(Colors.RGB, Colors.HSV(newhsv[1], newhsv[2], newhsv[3]))
    sethue(newrgb)
    t.pencolor = (newrgb.r, newrgb.g, newrgb.b)
end

"""
Set the Red, Green, and Blue colors of the turtle:

    Pencolor(t::Turtle, r, g, b)
"""
function Pencolor(t::Turtle, r, g, b)
    sethue(r, g, b)
    t.pencolor = (r, g, b)
end

"""
Reposition: pick the turtle up and place it at another position:

    Reposition(t::Turtle, x, y)
"""
function Reposition(t::Turtle, x, y)
    t.xpos = x
    t.ypos = y
end

"""
Set the width of the line drawn.

    Penwidth(t::Turtle, w)
"""
Penwidth(t::Turtle, w) = setline(w)

"""
Change the opacity of the pen to some value at random.

"""
Pen_opacity_random(t::Turtle) = setopacity(rand())


text:
struct DynamicBudgetedBipartiteMatchingSolution{T, U} <:
       BudgetedBipartiteMatchingSolution{T, U}
    instance::BudgetedBipartiteMatchingInstance{T, U}
    solution::Vector{Edge{T}}
    value::Float64
    state::Dict{Tuple{Int, Int, Int, Int}, Float64} # Data structure built by the dynamic-programming recursion.
    solutions::Dict{Tuple{Int, Int, Int, Int}, Vector{Edge{T}}} # From the indices of state to the corresponding solution.
end

function matching_dp_budgeted(
    instance::BudgetedBipartiteMatchingInstance{T, U},
) where {T, U}
    # First, determine the size of a perfect matching.
    nr = Dict(k => 1.0 for (k, v) in instance.matching.reward)
    ni = BipartiteMatchingInstance(instance.matching, nr)
    ns = matching_dp_imperfect(ni)
    m = length(ns.solution)

    # Solve this budgeted matching to ensure the constraints are met.

    # Four indices:
    # - two vertices, like matching_dp_imperfect;
    # - the number of edges still to add to have a perfect matching, like matching_dp;
    # - the minimum budget to respect.
    state = Dict{Tuple{Int, Int, Int, Int}, Float64}(
        (i, j, k, l) => -Inf for i in instance.matching.vertex_left,
        j in instance.matching.vertex_right, k in 0:m,
        l in 0:(instance.budget)
    )
    solutions = Dict{Tuple{Int, Int, Int, Int}, Vector{Edge{T}}}(
        (i, j, m, zero(U)) => Edge{T}[] for
        i in instance.matching.vertex_left,
        j in instance.matching.vertex_right
    )

    # Initialise: k=m-1, l=0, i.e. take at most one edge, don't consider the budget yet.
    v1 = first(instance.matching.vertex_right)
    for i in instance.matching.vertex_left
        if has_edge(instance.matching.graph, i, v1)
            state[i, v1, m - 1, 0] = instance.matching.reward[Edge(i, v1)]
            solutions[i, v1, m - 1, 0] = [Edge(i, v1)]
        end
    end

    v1 = first(instance.matching.vertex_left)
    for j in instance.matching.vertex_right
        if has_edge(instance.matching.graph, v1, j)
            state[v1, j, m - 1, 0] = instance.matching.reward[Edge(v1, j)]
            solutions[v1, j, m - 1, 0] = [Edge(v1, j)]
        end
    end

    for (i, prev_i) in zip(
        instance.matching.vertex_left[2:end],
        instance.matching.vertex_left[1:(end - 1)],
    )
        for (j, prev_j) in zip(
            instance.matching.vertex_right[2:end],
            instance.matching.vertex_right[1:(end - 1)],
        )
            # Take (i, j).
            if has_edge(instance.matching.graph, i, j)
                state[i, j, m - 1, 0] = instance.matching.reward[Edge(i, j)]
                solutions[i, j, m - 1, 0] = [Edge(i, j)]
            end

            # Don't take (i, j): copy a previous solution. No mixing possible, as only one edge is considered.
            if state[prev_i, j, m - 1, 0] > state[i, j, m - 1, 0]
                state[i, j, m - 1, 0] = state[prev_i, j, m - 1, 0]
                solutions[i, j, m - 1, 0] = solutions[prev_i, j, m - 1, 0]
            end
            if state[i, prev_j, m - 1, 0] > state[i, j, m - 1, 0]
                state[i, j, m - 1, 0] = state[i, prev_j, m - 1, 0]
                solutions[i, j, m - 1, 0] = solutions[i, prev_j, m - 1, 0]
            end
        end
    end

    # Initialisation for l=0, k < m-1: create perfect matchings, regardless of budget.
    for k in (m - 2):-1:0
        # Don't handle the cases with just one vertex in one part of the graph:
        # this loop requires at least two edges to be taken, but the considered
        # subset of edges can lead to at most only one being taken.

        # Consider the other vertices.
        for (i, prev_i) in zip(
            instance.matching.vertex_left[2:end],
            instance.matching.vertex_left[1:(end - 1)],
        )
            for (j, prev_j) in zip(
                instance.matching.vertex_right[2:end],
                instance.matching.vertex_right[1:(end - 1)],
            )
                # Take (i, j).
                if has_edge(instance.matching.graph, i, j) &&
                   state[prev_i, prev_j, k + 1, 0] > -Inf
                    state[i, j, k, 0] =
                        instance.matching.reward[Edge(i, j)] +
                        state[prev_i, prev_j, k + 1, 0]
                    solutions[i, j, k, 0] = solutions[prev_i, prev_j, k + 1, 0]
                    push!(solutions[i, j, k, 0], Edge(i, j))
                end

                # Don't take (i, j): mix previous solutions (as a special case, maybe take a previous solution completely).
                for k2 in 0:m
                    for k3 in 0:m
                        if state[prev_i, j, k2, 0] + state[i, prev_j, k3, 0] >
                           state[i, j, k, 0]
                            potential_sol = union(
                                Set{Edge{T}}(solutions[prev_i, j, k2, 0]),
                                Set{Edge{T}}(solutions[i, prev_j, k3, 0]),
                            )
                            val = sum(
                                instance.matching.reward[Edge(src(e), dst(e))] for e in potential_sol
                            )

                            if length(potential_sol) == m - k &&
                               val > state[i, j, k, 0]
                                state[i, j, k, 0] = val
                                solutions[i, j, k, 0] = collect(potential_sol)
                            end
                        end
                    end
                end
            end
        end
    end

    # Main iteration, considering the first vertex of each side is already done.
    for l in 1:(instance.budget)
        # Repeat the loop as often as necessary to handle the corner case with zero-weight edges.
        while true
            changes = false

            # Deal with the first two vertices (one on each side). All constraints are imposed.
            # Just one edge to be taken, as one vertex is always fixed on one side:
            # no feasible solution for k < m-1.
            vl = first(instance.matching.vertex_left)
            vr = first(instance.matching.vertex_right)
            if has_edge(instance.matching.graph, vl, vr) &&
               instance.weight[Edge(vl, vr)] >= l
                state[vl, vr, m - 1, l] = instance.matching.reward[Edge(vl, vr)]
                solutions[vl, vr, m - 1, l] = [Edge(vl, vr)]
            end

            v1 = first(instance.matching.vertex_right)
            for i in instance.matching.vertex_left
                if has_edge(instance.matching.graph, i, v1) &&
                   instance.weight[Edge(i, v1)] >= l
                    state[i, v1, m - 1, l] =
                        instance.matching.reward[Edge(i, v1)]
                    solutions[i, v1, m - 1, l] = [Edge(i, v1)]
                end
            end

            v1 = first(instance.matching.vertex_left)
            for j in instance.matching.vertex_right
                if has_edge(instance.matching.graph, v1, j) &&
                   instance.weight[Edge(v1, j)] >= l
                    state[v1, j, m - 1, l] =
                        instance.matching.reward[Edge(v1, j)]
                    solutions[v1, j, m - 1, l] = [Edge(v1, j)]
                end
            end

            # Do the rest.
            for k in (m - 1):-1:0
                for (i, prev_i) in zip(
                    instance.matching.vertex_left[2:end],
                    instance.matching.vertex_left[1:(end - 1)],
                )
                    for (j, prev_j) in zip(
                        instance.matching.vertex_right[2:end],
                        instance.matching.vertex_right[1:(end - 1)],
                    )

                        # Take (i, j) if possible.
                        if has_edge(instance.matching.graph, i, j)
                            remaining_budget =
                                max(0, l - instance.weight[Edge(i, j)])

                            if state[prev_i, prev_j, k + 1, remaining_budget] >
                               -Inf
                                state[i, j, k, l] =
                                    instance.matching.reward[Edge(i, j)] +
                                    state[
                                        prev_i,
                                        prev_j,
                                        k + 1,
                                        remaining_budget,
                                    ]
                                solutions[i, j, k, l] = solutions[
                                    prev_i,
                                    prev_j,
                                    k + 1,
                                    remaining_budget,
                                ]
                                push!(solutions[i, j, k, l], Edge(i, j))

                                changes = true
                            end
                        end

                        # Don't take (i, j): mix previous solutions (as a special case, maybe take a previous solution completely).
                        for k2 in 0:m
                            for k3 in 0:m
                                if state[prev_i, j, k2, l] +
                                   state[i, prev_j, k3, l] > state[i, j, k, l]
                                    potential_sol = union(
                                        Set{Edge{T}}(
                                            solutions[prev_i, j, k2, l],
                                        ),
                                        Set{Edge{T}}(
                                            solutions[i, prev_j, k3, l],
                                        ),
                                    )
                                    val = sum(
                                        instance.matching.reward[Edge(
                                            src(e),
                                            dst(e),
                                        )] for e in potential_sol
                                    )
                                    wei = sum(
                                        instance.weight[Edge(src(e), dst(e))] for e in potential_sol
                                    )

                                    if length(potential_sol) == m - k &&
                                       val > state[i, j, k, l] &&
                                       wei >= l
                                        state[i, j, k, l] = val
                                        solutions[i, j, k, l] =
                                            collect(potential_sol)

                                        changes = true
                                    end
                                end
                            end
                        end
                    end
                end
            end

            # If this iteration made no change, stop and go on to the next budget.
            if !changes
                break
            end
        end
    end

    val = state[
        last(instance.matching.vertex_left),
        last(instance.matching.vertex_right),
        0,
        instance.budget,
    ]
    sol = solutions[
        last(instance.matching.vertex_left),
        last(instance.matching.vertex_right),
        0,
        instance.budget,
    ]

    return DynamicBudgetedBipartiteMatchingSolution(
        instance,
        sol,
        val,
        state,
        solutions,
    )
    # return SimpleBudgetedBipartiteMatchingSolution(instance, sol)
end


text:
using InverseLaplace
import InverseLaplace: talbot, gwr
using Test

include("weeks_test.jl")
include("interface_test.jl")

@testset "plain function interface" begin
    @test isapprox(talbot(s -> 1/s, 1.0) , 1.000000000007737)
    @test isapprox(talbot(s -> 1/s, 1), 1.0)
end

@testset "ilt wraps talbot" begin
    @test isapprox(talbot(s -> 1/s, 1), ilt(s -> 1/s, 1))
    @test isapprox(gwr(s -> 1/s, 1), 1.0)
    @test isapprox(gwr(s -> 1/s, 1.0), 1.0000000002465594)
    @test isapprox(talbot(s -> s/(1+s^2), 1), cos(1))
    @test isapprox(talbot(s -> 1/(1+s), 1), exp(-1))
    @test isapprox(gwr(s -> 1/(1+s), 1), exp(-1))
end

@testset "ilt wraps various" begin
    if Int != Int32 && VERSION >= v"0.5.0"
        @test isapprox(talbot(s -> 1/s^4, 1//10) * 6, 1e-3)
        @test isapprox(gwr(s -> 1/s^4, 1//10) * 6, 0.0010000012595365085; atol = 1e-9)

        # At one point, talbot returned rational with rational input. Now we convert to BigFloat
        @test isapprox(talbot.(s -> 1/s^3, [1//10, 2//10], 64), [5e-3, 2e-2]; atol = 1e-18)
        @test (talbot.(s -> 1/s^3, [[.1, .2] ; [.4, .5]]); true)
        @test isapprox(talbot.(s -> 1/s^3, [.1 .2; .3 .4] ),  [0.005 0.02; 0.045 0.08], atol = 1e-13)
    else
        @test isapprox(talbot(s -> 1/s^4, 1//10) * 6, 1e-3; atol = 1e-9)
        @test isapprox(gwr(s -> 1/s^4, 1//10) * 6, 0.0010000012595365085; atol = 1e-9)
        @test isapprox(talbot.(s -> 1/s^3, [.1 .2; .3 .4] ), [0.005 0.02; 0.045 0.08]; atol = 1e-9)
    end
end

@testset "Talbot array" begin
    @test isapprox(talbot.(s -> 1/s^3, [.1, .2]), [0.005, 0.02])
    @test isapprox(InverseLaplace.talbotarr(s -> 1/s^3, [.1, .2]), [0.005, 0.02])
end

# gwr is often less accurate.

@testset "broadcast gwr" begin
    if VERSION >= v"0.5.0"
        @test isapprox(gwr.(s -> 1/s^3, [.1, .2]), [0.005, 0.02]; atol = 1e-10)
    else
        @test isapprox(gwr.(s -> 1/s^3, [.1, .2]), [0.005, 0.02]; atol = 1e-7)
    end
end

@testset "check types" begin
    @test isapprox(map(x -> convert(Float64, x),
                       talbot.(s -> 1/s^3, [BigFloat(1//10), BigFloat(2//10)])), [0.005, 0.02])
    @test typeof(talbot(s -> 1/s^3, 1//10)) == BigFloat
    @test typeof(talbot(s -> 1/s^3, 1//10, 64)) == BigFloat
end

@time @testset "gaverstehfest" begin include("gaverstehfest_tests.jl") end


text:
include("base.jl")

function analyze(widget::QGroup, ctx::PlotContext)
    local dates = Array{DateTime,1}()
    local srcCtx = Dict{String,DataAttributeContext}()
    local mini = typemax(Float64)
    local maxi = typemin(Float64)
    local titles = Array{String,1}()

    for attr in widget.data
        local _dates, _srcCtx, _mini, _maxi, _titles = analyze(attr, ctx)

        if widget.ctxs !== nothing
            if widget.ctxs |> length === 0
                _srcCtx = Dict()
            end
            for c in widget.ctxs
                if "$c" in _srcCtx
                    mini = min(mini, _mini)
                    maxi = max(maxi, _maxi)
                    delete!(_srcCtx, "$c")
                end
            end
        end

        for (k, v) in _srcCtx
            srcCtx[k] = v
        end
        dates = _dates
        titles = [titles..., _titles...]
    end

    return dates, srcCtx, mini, maxi, titles
end

function renderWidget!(widget::QGroup, ctx::PlotContext)
    # Check if this group is a specialication
    local dates, _, mini, maxi, _ = analyze(widget, ctx)
    local numbers = nothing

    local oPlot = ctx.plot
    local myStyle = ctx.styles
    ctx.group += 1

    if ctx.group == 2
        plot!(oPlot, right_margin = 25mm) # increase margin in parent
        ctx.plot = twinx(ctx.plot)
        plot!(ctx.plot, y_formatter = Y_FORMATTER) # set formatter in child
    end
    plot!(ctx.plot, ylabel = widget.unit)
    if widget.log
        plot!(ctx.plot, yaxis = :log)
    end

    @assert ctx.group <= 2 "Currently only two groups are supported"
    local yMax = maxi
    local yMin = mini
    if widget.min !== nothing
        yMin = widget.min
    end
    if widget.log
        @assert widget.min === nothing || widget.min > 0
        if yMin <= 0
            yMin = 1
        end
        ctx.styles["log"] = yMin
    else
        yMin = min(yMin, 0)
    end

    if widget.max !== nothing
        yMax = widget.max
    end

    for w in widget.data
        ctx.styles = deepcopy(myStyle)
        delete!(ctx.styles, "source context")
        local _, srcCtx, _, _, _ = analyze(QGroup("", data = widget.data), ctx)

        local ctxs = []
        if widget.ctxs === nothing
            ctxs = [nothing]
        else
            ctxs = collect(values(srcCtx))
            if length(widget.ctxs) > 0
                ctxs = filter(y -> "$y" in map(x -> "$x", widget.ctxs), ctxs)
            end
        end
        # TODO: Does only support single source contexts!

        for c in ctxs
            try
                ctx.styles = deepcopy(myStyle)
                ctx.styles["ylims"] = (yMin, yMax)
                if c !== nothing
                    logger(c, "Plotting where group context is", true)
                    ctx.styles["source context"] = c
                    ctx.styles["append to label"] = ""
                    for cc in c.ctx
                        ctx.styles["append to label"] *= "$(cc[2])=$(cc[3])"
                    end
                end
                _, numbers = renderWidget!(w, ctx)
            catch exception
                @show "Exception in groups"
                bt = catch_backtrace()
                msg = sprint(showerror, exception, bt)
                logger(msg, "failed", true, failed = true)
                # TODO: only catch context-not-found exceptions
            end
        end
    end

    ctx.plot = oPlot

    return dates, numbers
end

max_stars_repo_name: Agapanthus/server-diary

max_issues_repo_name: Agapanthus/server-diary

max_forks_repo_name: Agapanthus/server-diary


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct NetworkInterfacesApi <: SwaggerApi
    client::Swagger.Client
end

"""

Gets the specified network interface ip configuration.
Param: resourceGroupName::String (required)
Param: networkInterfaceName::String (required)
Param: ipConfigurationName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NetworkInterfaceIPConfiguration
"""
function _swaggerinternal_networkInterfaceIPConfigurationsGet(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, ipConfigurationName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", NetworkInterfaceIPConfiguration, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkInterfaces/{networkInterfaceName}/ipConfigurations/{ipConfigurationName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkInterfaceName", networkInterfaceName)  # type String
    Swagger.set_param(_ctx.path, "ipConfigurationName", ipConfigurationName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkInterfaceIPConfigurationsGet(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, ipConfigurationName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfaceIPConfigurationsGet(_api, resourceGroupName, networkInterfaceName, ipConfigurationName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkInterfaceIPConfigurationsGet(_api::NetworkInterfacesApi, response_stream::Channel, resourceGroupName::String, networkInterfaceName::String, ipConfigurationName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfaceIPConfigurationsGet(_api, resourceGroupName, networkInterfaceName, ipConfigurationName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Get all ip configurations in a network interface.
Param: resourceGroupName::String (required)
Param: networkInterfaceName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NetworkInterfaceIPConfigurationListResult
"""
function _swaggerinternal_networkInterfaceIPConfigurationsList(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", NetworkInterfaceIPConfigurationListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkInterfaces/{networkInterfaceName}/ipConfigurations", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkInterfaceName", networkInterfaceName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkInterfaceIPConfigurationsList(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfaceIPConfigurationsList(_api, resourceGroupName, networkInterfaceName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkInterfaceIPConfigurationsList(_api::NetworkInterfacesApi, response_stream::Channel, resourceGroupName::String, networkInterfaceName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfaceIPConfigurationsList(_api, resourceGroupName, networkInterfaceName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

List all load balancers in a network interface.
Param: resourceGroupName::String (required)
Param: networkInterfaceName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NetworkInterfaceLoadBalancerListResult
"""
function _swaggerinternal_networkInterfaceLoadBalancersList(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", NetworkInterfaceLoadBalancerListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkInterfaces/{networkInterfaceName}/loadBalancers", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkInterfaceName", networkInterfaceName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkInterfaceLoadBalancersList(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfaceLoadBalancersList(_api, resourceGroupName, networkInterfaceName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkInterfaceLoadBalancersList(_api::NetworkInterfacesApi, response_stream::Channel, resourceGroupName::String, networkInterfaceName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfaceLoadBalancersList(_api, resourceGroupName, networkInterfaceName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Creates or updates a Tap configuration in the specified NetworkInterface.
Param: resourceGroupName::String (required)
Param: networkInterfaceName::String (required)
Param: tapConfigurationName::String (required)
Param: tapConfigurationParameters::NetworkInterfaceTapConfiguration (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NetworkInterfaceTapConfiguration
"""
function _swaggerinternal_networkInterfaceTapConfigurationsCreateOrUpdate(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, tapConfigurationName::String, tapConfigurationParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", NetworkInterfaceTapConfiguration, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkInterfaces/{networkInterfaceName}/tapConfigurations/{tapConfigurationName}", ["azure_auth"], tapConfigurationParameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkInterfaceName", networkInterfaceName)  # type String
    Swagger.set_param(_ctx.path, "tapConfigurationName", tapConfigurationName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkInterfaceTapConfigurationsCreateOrUpdate(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, tapConfigurationName::String, tapConfigurationParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfaceTapConfigurationsCreateOrUpdate(_api, resourceGroupName, networkInterfaceName, tapConfigurationName, tapConfigurationParameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkInterfaceTapConfigurationsCreateOrUpdate(_api::NetworkInterfacesApi, response_stream::Channel, resourceGroupName::String, networkInterfaceName::String, tapConfigurationName::String, tapConfigurationParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfaceTapConfigurationsCreateOrUpdate(_api, resourceGroupName, networkInterfaceName, tapConfigurationName, tapConfigurationParameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Deletes the specified tap configuration from the NetworkInterface.
Param: resourceGroupName::String (required)
Param: networkInterfaceName::String (required)
Param: tapConfigurationName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_networkInterfaceTapConfigurationsDelete(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, tapConfigurationName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkInterfaces/{networkInterfaceName}/tapConfigurations/{tapConfigurationName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkInterfaceName", networkInterfaceName)  # type String
    Swagger.set_param(_ctx.path, "tapConfigurationName", tapConfigurationName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkInterfaceTapConfigurationsDelete(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, tapConfigurationName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfaceTapConfigurationsDelete(_api, resourceGroupName, networkInterfaceName, tapConfigurationName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkInterfaceTapConfigurationsDelete(_api::NetworkInterfacesApi, response_stream::Channel, resourceGroupName::String, networkInterfaceName::String, tapConfigurationName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfaceTapConfigurationsDelete(_api, resourceGroupName, networkInterfaceName, tapConfigurationName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Get the specified tap configuration on a network interface.
Param: resourceGroupName::String (required)
Param: networkInterfaceName::String (required)
Param: tapConfigurationName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NetworkInterfaceTapConfiguration
"""
function _swaggerinternal_networkInterfaceTapConfigurationsGet(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, tapConfigurationName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", NetworkInterfaceTapConfiguration, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkInterfaces/{networkInterfaceName}/tapConfigurations/{tapConfigurationName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkInterfaceName", networkInterfaceName)  # type String
    Swagger.set_param(_ctx.path, "tapConfigurationName", tapConfigurationName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkInterfaceTapConfigurationsGet(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, tapConfigurationName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfaceTapConfigurationsGet(_api, resourceGroupName, networkInterfaceName, tapConfigurationName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkInterfaceTapConfigurationsGet(_api::NetworkInterfacesApi, response_stream::Channel, resourceGroupName::String, networkInterfaceName::String, tapConfigurationName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfaceTapConfigurationsGet(_api, resourceGroupName, networkInterfaceName, tapConfigurationName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Get all Tap configurations in a network interface.
Param: resourceGroupName::String (required)
Param: networkInterfaceName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NetworkInterfaceTapConfigurationListResult
"""
function _swaggerinternal_networkInterfaceTapConfigurationsList(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", NetworkInterfaceTapConfigurationListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkInterfaces/{networkInterfaceName}/tapConfigurations", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkInterfaceName", networkInterfaceName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkInterfaceTapConfigurationsList(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfaceTapConfigurationsList(_api, resourceGroupName, networkInterfaceName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkInterfaceTapConfigurationsList(_api::NetworkInterfacesApi, response_stream::Channel, resourceGroupName::String, networkInterfaceName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfaceTapConfigurationsList(_api, resourceGroupName, networkInterfaceName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Creates or updates a network interface.
Param: resourceGroupName::String (required)
Param: networkInterfaceName::String (required)
Param: parameters::NetworkInterface (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NetworkInterface
"""
function _swaggerinternal_networkInterfacesCreateOrUpdate(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", NetworkInterface, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkInterfaces/{networkInterfaceName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkInterfaceName", networkInterfaceName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkInterfacesCreateOrUpdate(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesCreateOrUpdate(_api, resourceGroupName, networkInterfaceName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkInterfacesCreateOrUpdate(_api::NetworkInterfacesApi, response_stream::Channel, resourceGroupName::String, networkInterfaceName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesCreateOrUpdate(_api, resourceGroupName, networkInterfaceName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Deletes the specified network interface.
Param: resourceGroupName::String (required)
Param: networkInterfaceName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_networkInterfacesDelete(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkInterfaces/{networkInterfaceName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkInterfaceName", networkInterfaceName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkInterfacesDelete(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesDelete(_api, resourceGroupName, networkInterfaceName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkInterfacesDelete(_api::NetworkInterfacesApi, response_stream::Channel, resourceGroupName::String, networkInterfaceName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesDelete(_api, resourceGroupName, networkInterfaceName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets information about the specified network interface.
Param: resourceGroupName::String (required)
Param: networkInterfaceName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: expand::String
Return: NetworkInterface
"""
function _swaggerinternal_networkInterfacesGet(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", NetworkInterface, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkInterfaces/{networkInterfaceName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkInterfaceName", networkInterfaceName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_param(_ctx.query, "\$expand", expand)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkInterfacesGet(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesGet(_api, resourceGroupName, networkInterfaceName, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkInterfacesGet(_api::NetworkInterfacesApi, response_stream::Channel, resourceGroupName::String, networkInterfaceName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesGet(_api, resourceGroupName, networkInterfaceName, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets all route tables applied to a network interface.
Param: resourceGroupName::String (required)
Param: networkInterfaceName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: EffectiveRouteListResult
"""
function _swaggerinternal_networkInterfacesGetEffectiveRouteTable(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", EffectiveRouteListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkInterfaces/{networkInterfaceName}/effectiveRouteTable", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkInterfaceName", networkInterfaceName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkInterfacesGetEffectiveRouteTable(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesGetEffectiveRouteTable(_api, resourceGroupName, networkInterfaceName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkInterfacesGetEffectiveRouteTable(_api::NetworkInterfacesApi, response_stream::Channel, resourceGroupName::String, networkInterfaceName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesGetEffectiveRouteTable(_api, resourceGroupName, networkInterfaceName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets all network interfaces in a resource group.
Param: resourceGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NetworkInterfaceListResult
"""
function _swaggerinternal_networkInterfacesList(_api::NetworkInterfacesApi, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", NetworkInterfaceListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkInterfaces", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkInterfacesList(_api::NetworkInterfacesApi, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesList(_api, resourceGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkInterfacesList(_api::NetworkInterfacesApi, response_stream::Channel, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesList(_api, resourceGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets all network interfaces in a subscription.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NetworkInterfaceListResult
"""
function _swaggerinternal_networkInterfacesListAll(_api::NetworkInterfacesApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", NetworkInterfaceListResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Network/networkInterfaces", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkInterfacesListAll(_api::NetworkInterfacesApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesListAll(_api, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkInterfacesListAll(_api::NetworkInterfacesApi, response_stream::Channel, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesListAll(_api, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets all network security groups applied to a network interface.
Param: resourceGroupName::String (required)
Param: networkInterfaceName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: EffectiveNetworkSecurityGroupListResult
"""
function _swaggerinternal_networkInterfacesListEffectiveNetworkSecurityGroups(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", EffectiveNetworkSecurityGroupListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkInterfaces/{networkInterfaceName}/effectiveNetworkSecurityGroups", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkInterfaceName", networkInterfaceName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkInterfacesListEffectiveNetworkSecurityGroups(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesListEffectiveNetworkSecurityGroups(_api, resourceGroupName, networkInterfaceName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkInterfacesListEffectiveNetworkSecurityGroups(_api::NetworkInterfacesApi, response_stream::Channel, resourceGroupName::String, networkInterfaceName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesListEffectiveNetworkSecurityGroups(_api, resourceGroupName, networkInterfaceName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Updates a network interface tags.
Param: resourceGroupName::String (required)
Param: networkInterfaceName::String (required)
Param: parameters::TagsObject (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NetworkInterface
"""
function _swaggerinternal_networkInterfacesUpdateTags(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PATCH", NetworkInterface, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkInterfaces/{networkInterfaceName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkInterfaceName", networkInterfaceName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkInterfacesUpdateTags(_api::NetworkInterfacesApi, resourceGroupName::String, networkInterfaceName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesUpdateTags(_api, resourceGroupName, networkInterfaceName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkInterfacesUpdateTags(_api::NetworkInterfacesApi, response_stream::Channel, resourceGroupName::String, networkInterfaceName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkInterfacesUpdateTags(_api, resourceGroupName, networkInterfaceName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export networkInterfaceIPConfigurationsGet, networkInterfaceIPConfigurationsList, networkInterfaceLoadBalancersList, networkInterfaceTapConfigurationsCreateOrUpdate, networkInterfaceTapConfigurationsDelete, networkInterfaceTapConfigurationsGet, networkInterfaceTapConfigurationsList, networkInterfacesCreateOrUpdate, networkInterfacesDelete, networkInterfacesGet, networkInterfacesGetEffectiveRouteTable, networkInterfacesList, networkInterfacesListAll, networkInterfacesListEffectiveNetworkSecurityGroups, networkInterfacesUpdateTags


text:
#=

  From
  "CAV 2020 Tutorial: Probabilistic Programming: A Guide for Verificationists"
  https://www.youtube.com/watch?v=yz5uUf_03Ik&t=2657s
  Around @23:30

  Summary Statistics
  parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

          pr    0.3290    0.2330     0.0023    0.0036   4184.9815    0.9999      900.7709
           p    0.0000    0.0000     0.0000    0.0000         NaN       NaN           NaN
           x   14.0153    0.4955     0.0050    0.0062   5962.5178    0.9999     1283.3659

  cf ~/webppl/gaussian_mixture_model2.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function gaussian_mixture_model2(xval=14.5)
    pr ~ Uniform(0,1)
    p ~ flip(pr)

    # x ~ p ? Normal(10.0,1.0) : Normal(14.0,0.5)
    # x ~ p ? Normal(10.0,1.0) : Normal(14.0,1)

    # Note: Both these approaches give incorrect solutions
    # x == 14.5 || begin Turing.@addlogprob! -Inf; return end
    # true ~ Dirac(x == 14.5)
    # true ~ Dirac(x == 8.5)
    # true ~ Dirac(x == 12.5)

    # This is the way, "observe externally"
    xval ~ p ? Normal(10.0,1.0) : Normal(14.0,0.5)

    # Post
    x_post ~ p ? Normal(10.0,1.0) : Normal(14.0,0.5)

end

model = gaussian_mixture_model2()

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 100_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, NUTS(), 10_000)
# chns = sample(model, HMC(0.1,5), 10_000)
# chns = sample(model, Gibbs(HMC(0.1,5,:pr,:x),PG(5,:p)), 10_000)
# chns = sample(model, Gibbs(HMC(0.1,5,:pr,:x),MH(:p)), 10_000)
# chns = sample(model, Gibbs(NUTS(1000,0.65,:pr,:x),PG(15,:p)), 10_000)

display(chns)
# display(plot(chns))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function background_text(str_array)
    colorband = diverging_palette(0, 10, 100)
    gsave()
    setopacity(0.5)
    x = -currentwidth/2
    y = -currentheight/2
    fontsize(12)
    while y < currentheight/2
        sethue(colorband[rand(1:end)])
        s = str_array[rand(1:end)]
        text(s, x, y)
        se = textextents(s)
        x += se[5] # move to the right
        if x > currentwidth/2
           x = -currentwidth/2 # next row
           y += 10
        end
    end
    grestore()
end

function heart()
    move(127, 1) # damn, it's offset from 0/0, who drew this :/
    curve(134.2, -6.5, 134.2, -6.5, 156.1, -29.6)
    curve(185.8, -60.5, 198.1, -74.3, 213.7, -95.3)
    curve(240.4, -131, 253.3, -163.7, 253.3, -194.9)
    curve(253.3, -218, 246.4, -237.8, 232.6, -253.7)
    curve(219.1, -268.7, 204.1, -275.3, 181.9, -275.3)
    curve(154, -275.3, 136.3, -265.1, 127, -243.8)
    curve(124, -252.5, 120.4, -257.6, 112.9, -263.6)
    curve(103.6, -270.8, 88.3, -275.3, 73.3, -275.3)
    curve(59.2, -275.3, 46, -271.4, 35.2, -264.5)
    curve(14.5, -250.7, 1, -223.4, 1, -194.6)
    curve(1, -167.3, 13, -136.4, 37.3, -101)
    curve(53.8, -77, 86.5, -39.8, 127, 1)
    closepath()
end

function heart_with_julias(x=0, y=0)
    gsave()
    translate(x, y)
    setcolor("lavenderblush")
    heart()
    fillpreserve()
    clip()
    translate(-50, -300)
    for y in 0:30:500
        gsave()
        for x in 0:30:250
            translate(30, 0)
            gsave()
                scale(0.1, 0.1)
                julialogo()
            grestore()
        end
        grestore()
        translate(0, 20)
    end
    clipreset()
    grestore()
end

function outlined_heart()
    gsave()
    scale(1.2, 1.2)
    translate(-127, -30) # must fix that x-offset one day
    heart_with_julias()
    heart()
    setline(4)
    setcolor(1, 0, 0)
    strokepath()
    grestore()
end

function julia_heart(fname)
    global currentwidth = 1928 # pts
    global currentheight = 1064 # pts
    Drawing(currentwidth, currentheight, fname)

    origin()
    background("black")
    namelist = map(x->string(x), names(Base)) # list of names in Base.

    background_text(namelist)
    for n in 1:5
        rotate(2pi/5)
        outlined_heart()
    end
    @test finish() == true
    println("...finished test: output in $(fname)")
end

julia_heart("heart-julia.pdf")


text:
using MAT
using NeuralNetworkAnalysis
using NeuralNetworkAnalysis: FunctionPreprocessing
using Plots

## Model
const u = 0.0001  # friction parameter
const a_lead = -2.0  # acceleration control input applied to the lead vehicle

@taylorize function ACC!(dx, x, p, t)
    v_lead = x[2]  # lead car velocity
    γ_lead = x[3]  # lead car acceleration
    v_ego = x[5]  # ego car velocity
    γ_ego = x[6]  # ego car acceleration
    a_ego = x[7]  # ego car acceleration control input

    # lead car dynamics
    dx[1] = v_lead
    dx[2] = γ_lead
    dx[3] = 2 * (a_lead - γ_lead) - u * v_lead^2

    # ego car dynamics
    dx[4] = v_ego
    dx[5] = γ_ego
    dx[6] = 2 * (a_ego - γ_ego) - u * v_ego^2
    dx[7] = zero(a_ego)
    return dx
end

## Specification
# We choose the controller with 5 hidden layers.
controller = read_nnet_mat(@modelpath("ACC", "controller_5_20.mat");
                           act_key="act_fcns");

# The initial states according to the specification are:
X₀ = Hyperrectangle(low=[90, 32, 0, 10, 30, 0],
                    high=[110, 32.2, 0, 11, 30.2, 0])
U₀ = ZeroSet(1);

# The system has 6 state variables and 1 control variable:
vars_idx = Dict(:state_vars=>1:6, :control_vars=>7)
ivp = @ivp(x' = ACC!(x), dim: 7, x(0) ∈ X₀ × U₀)
period = 0.1;  # control period

# Preprocessing function for the network input:
v_set = 30.0  # ego car's set speed
T_gap = 1.4
M = zeros(3, 6)
M[1, 5] = 1.0
M[2, 1] = 1.0
M[2, 4] = -1.0
M[3, 2] = 1.0
M[3, 5] = -1.0
function preprocess(X::LazySet)  # version for set computations
    Y1 = Singleton([v_set, T_gap])
    Y2 = linear_map(M, X)
    return cartesian_product(Y1, Y2)
end
function preprocess(X::AbstractVector)  # version for simulations
    Y1 = [v_set, T_gap]
    Y2 = M * X
    return vcat(Y1, Y2)
end
control_preprocessing = FunctionPreprocessing(preprocess)

prob = ControlledPlant(ivp, controller, vars_idx, period;
                       preprocessing=control_preprocessing);

# Safety specification
T = 5.0  # time horizon

D_default = 10.0
d_rel = [1.0, 0, 0, -1, 0, 0, 0]
d_safe = [0, 0, 0, 0, T_gap, 0, 0]

d_prop = d_rel - d_safe
safe_states = HalfSpace(-d_prop, -D_default)
predicate = X -> X ⊆ safe_states

## Results
alg = TMJets(abstol=1e-6, orderT=6, orderQ=1)
alg_nn = Ai2()

function benchmark(; silent::Bool=false)
    # We solve the controlled system:
    silent || println("flowpipe construction")
    res_sol = @timed solve(prob, T=T, alg_nn=alg_nn, alg=alg)
    sol = res_sol.value
    silent || print_timed(res_sol)

    # Next we check the property for an overapproximated flowpipe:
    silent || println("property checking")
    solz = overapproximate(sol, Zonotope)
    res_pred = @timed predicate(solz)
    silent || print_timed(res_pred)
    if res_pred.value
        silent || println("The property is satisfied.")
    else
        silent || println("The property may be violated.")
    end
    return solz
end

benchmark(silent=true)  # warm-up
res = @timed benchmark()  # benchmark
sol = res.value
println("total analysis time")
print_timed(res)

## Plot the Results
fig = plot(leg=(0.4, 0.3))
xlabel!(fig, "time")
F = flowpipe(sol)

fp_rel = linear_map(Matrix(d_rel'), F)
output_map_rel = d_rel

fp_safe = affine_map(Matrix(d_safe'), [D_default], F)
output_map_safe = vcat([D_default], d_safe)

plot!(fig, fp_rel, vars=(0, 1), c=:red, alpha=.4)
plot!(fig, fp_safe, vars=(0, 1), c=:blue, alpha=.4)


text:
module MangalTestInteractions

    using Mangal
    using Test

    @test typeof(interactions()) <: Vector{MangalInteraction}
    @test typeof(interactions(Pair("type", "mutualism"))) <: Vector{MangalInteraction}

    @test typeof(interactions(:, node(31863))) <: Vector{MangalInteraction}
    @test typeof(interactions(node(31904), :)) <: Vector{MangalInteraction}

    @test typeof(interactions(:, node(31863), Pair("type", "predation"))) <: Vector{MangalInteraction}
    @test typeof(interactions(node(31904), :, Pair("type", "predation"))) <: Vector{MangalInteraction}

    @test typeof(interactions(node(31904), node(31863))) <: Vector{MangalInteraction}
    @test typeof(interactions(node(31904), node(31863), Pair("type", "predation"))) <: Vector{MangalInteraction}

end


text:
#=
This file contains functions for the infinite alleles model in support of the Nearly Neutral paper.
Author:  Alden H. Wright, Department of Computer Science, Univeristy of Montana, Missoula, MT 59812 USA
    alden.wright@umontana.edu
Example run:
[experiments]$ julia run.jl examples/ia_example1
=#
using QuadGK
export fitness, inf_alleles, pop_counts64,
    dfe_deleterious,dfe_advantageous, dfe_mixed, dfe_mod, dfe_neutral, dfe_fixed

#=
function print_trial_result( tr::InfAlleles.trial_result )
  if tr.nn_simtype == 1
    println("\ninfinite alleleles model")
  else
    println("nn_simtype: ", tr.nn_simtype)
  end
  println("n: ", tr.n)
  println("N: ", tr.N)
  println("N_mu: ", tr.N_mu)
  println("mu: ", tr.N_mu/tr.N)
  println("ngens: ", tr.ngens)
  println("burn_in: ", tr.burn_in)
  println("dfe: ", tr.dfe)
  println("dfe_str: ", tr.dfe_str)
  println("mean fitness: ",tr.mean_fitness)
  println("expected richness: ",tr.expected_richness)
  println("average richness: ",tr.average_richness)
  println("stderr richness: ",tr.stderr_richness)
  println("expected_w_homoz: ",tr.expected_w_homoz)
  println("w_homoz: ",tr.w_homoz)
  println("stderr_w_homoz: ",tr.stderr_w_homoz)
  println("IQV: ",tr.IQV)
  println("stderr_IQV: ",tr.stderr_IQV)
  println("use_poplist: ",tr.use_poplist)
end
=#

#=
function add_expected_richness( tr::trial_result )
  sum = 0.0
  theta = 2*tr.N_mu
  for i = 0:(tr.N-1)
    sum += 1.0/(theta+i)
  end
  tr.expected_richness = theta*sum
end
=#

# Integrand of the equation 3.92 of Ewens (2004) which gives the population expectation of k (richness)
#   for the infinite alleles model
integrand(x,theta) = theta/x*(1-x)^(theta-1)
# Adds expected richness computed by eq. 3.92 of Ewens (2004) to tr.
function add_expected_richness( tr::trial_result )
  theta = 2*tr.N_mu
  println("N: ",tr.N,"  N_mu: ",tr.N_mu,"  integrand(1/tr.N,theta): ",integrand(1/tr.N,theta))
  result, err = quadgk( x->integrand(x,theta), 1.0/tr.N, 1.0 )
  try
    #println("N: ",tr.N,"  N_mu: ",tr.N_mu,"  integrand(1/tr.N,theta): ",integrand(1/tr.N,theta))
    result, err = quadgk( x->integrand(x,theta), 1.0/tr.N, 1.0 )
  catch
    println("computation of integrand failed. expected richness result will be incorrect.")
    result = 1.0
    err = 1.0
  end
  result += theta
  #println("add_expected richness result: ",result,"  err: ",err)
  tr.expected_richness = result
end

@doc """ function dfe_fitness( p::Int64, dfe::Function, fitness_table::Dict{Int64,Float64} )
Fitness function as saveed in dictionary fitness_table.
"""
function dfe_fitness( p::Int64, dfe::Function, fitness_table::Dict{Int64,Float64} )
  val = get(fitness_table,p,-1.0)
  if val == -1.0   # p is not in fitness table
    val = dfe( p )
    fitness_table[p] = val
  end
  val
end

@doc """ function inf_alleles( )
This function does the iteration over generations for the infinite alleles model.
Note:  dfe is "distribution of fitness effects" function.
In simulations for the nearly neutral paper, uniform_start and combine are both false.
"""
function inf_alleles( tr::trial_result; uniform_start::Bool=false, combine::Bool=false)
  global fitness_table = Dict{Int64,Float64}()
  g_limit = 1000000  # upper limit of generations to wait for extinctions and fixations
  int_burn_in = Int(round(tr.burn_in*tr.N/tr.N_mu+50.0))
  mu = tr.N_mu/tr.N
  if uniform_start  # All allele values start with the same value.  Start with a selective sweep.
    new_pop = Int64[1 for i = 1:tr.N]
    fit = dfe_fitness(1, tr.dfe, fitness_table )  # set fitness of 1 to be tr.dfe(1).
    new_id = 2
  else
    new_pop = collect(1:tr.N)
    new_id = 1
    for i = 1:tr.N
      # Note that dfe_advantageous(i), dfe_deleterious(i), dfe_mixed(i) do not depend on i, but dfe_mod(i) does depend on i
      fit = dfe_fitness(i, tr.dfe, fitness_table )  # set fitness of i to be tr.dfe(i).
      new_id += 1
    end
  end
  if tr.use_poplist
    poplist= Population[ new_pop ]
  end
  if combine
    pop_result = Population()
  end
  done = false
  # Initialize statistics accumulators
  sum_mean_fitness = 0.0
  sum_homoz = 0.0
  sum_sq_homoz = 0.0
  sum_IQV = 0.0
  sum_sq_IQV = 0.0
  sum_richness = 0.0
  sum_sq_richness = 0.0
  gcount = 0

  g = 2
  while !done && g < g_limit+tr.ngens+int_burn_in
    new_pop = deepcopy(new_pop)
    for i in 1:tr.N
      if rand() < mu
        new_pop[i] = new_id
        fit = dfe_fitness( new_id, tr.dfe, fitness_table )  # Set fitness of new_id
        new_id += 1
      end
    end
    new_pop = propsel( new_pop, tr.dfe, fitness_table )
    if g > int_burn_in

      # Accumulate statistics
      #println("mean fitness: ", mean( [ fitness_table[p] for p in new_pop ] ) )
      sum_mean_fitness += mean( [ fitness_table[p] for p in new_pop ] )
      pcounts = pop_counts64(new_pop)
      IQVvalue = IQV( pcounts )
      sum_IQV += IQVvalue
      sum_sq_IQV += IQVvalue^2
      homoz = watterson_homozygosity( pcounts )
      sum_homoz += homoz
      sum_sq_homoz += homoz^2
      richness = length(pcounts)
      sum_richness += richness
      sum_sq_richness += richness^2
      gcount += 1
    end

    if tr.use_poplist
      Base.push!( poplist, new_pop )
    end
    if combine && g >= int_burn_in+1 && g <= tr.ngens+int_burn_in
      pop_result = vcat( pop_result, new_pop )
    end
    g += 1
    done = (g > tr.ngens+int_burn_in)
  end
  #println("gcount: ",gcount)
  tr.mean_fitness = sum_mean_fitness/gcount
  tr.expected_w_homoz = 1.0/(1.0+2.0*tr.N_mu)
  tr.w_homoz = sum_homoz/gcount
  tr.stderr_w_homoz = sqrt((sum_sq_homoz/(gcount-1) - sum_homoz^2/gcount/(gcount-1))/gcount)
  tr.average_richness = sum_richness/gcount
  tr.stderr_richness = sqrt((sum_sq_richness/(gcount-1) - sum_richness^2/gcount/(gcount-1))/gcount)
  tr.IQV = sum_IQV/gcount
  tr.stderr_IQV = sqrt((sum_sq_IQV/(gcount-1) - sum_IQV^2/gcount/(gcount-1))/gcount)
  if tr.use_poplist
    if combine
      return [pop_result]
    else
      return poplist[int_burn_in+1:int_burn_in+tr.ngens]
    end
  end
end

@doc """ function pop_counts64( pop::Population )
Returns the sorted frequencies of the alleles of Population pop.
Example:  If pop = [5, 7, 9, 5, 4, 5, 7], then the returned list is [3, 2, 1, 1]
   because there are 3 5's, 2 7's, 1 9, and 1 4.  So the sum of the returned
   list is the length of the population.
"""
function pop_counts64( pop::Population )
  c = Dict{Int64,Int64}()
  for x in pop
    c[x] = get( c, x, 0 ) + 1
  end
  map( x->c[x], sort( unique(pop), by=x->c[x], rev=true ) )
end

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
module Metropolis

export metropolis, metropolis_uniform

using Distributions

"""
Generate a series of random samples from a given probability distribution
using the metropolis algorithm

# Arguments
- `distribution`: the target probability distribution
- `step`: the probability distribution of the steps
- `samples`: the number of samples
- `init`: the initial value of the series (starting point of the "walker")

# Returns

the series of randome samples from the distribution and the acceptance rate
"""
function metropolis(distribution::Function, step::Function, samples::Integer;
        init::Real=0.0)
    series = Vector{Float64}(undef, samples)

    position = init
    probability = distribution(position)
    series[1] = position

    acceptcount = 0
    for i in 2:samples
        newposition = position + step(position)
        newprobablity = distribution(newposition)

        if newprobablity / probability > 1 || rand() <= newprobablity / probability
            position, probability = newposition, newprobablity
            acceptcount += 1
        end

        series[i] = position
    end

    return series, acceptcount / (samples - 1)
end

"""
Metropolis algorithm with uniform step probability between `-stepsize` and `+stepsize`
"""
function metropolis_uniform(distribution::Function, stepsize::Real, samples::Integer;
        init::Real=0.0)
    return metropolis(distribution, x -> rand(Uniform(-stepsize, stepsize)), samples,
        init=init)
end

end

max_stars_repo_name: slhshamloo/comp-phys

max_issues_repo_name: slhshamloo/comp-phys

max_forks_repo_name: slhshamloo/comp-phys


text:
function currency_exchange_rate(from_currency::String, to_currency::String)
    uri = _form_uri_head("CURRENCY_EXCHANGE_RATE") * "&from_currency=$from_currency&to_currency=$to_currency&apikey=" * ENV["ALPHA_VANTAGE_API_KEY"]
    data = _get_request(uri)
    return _parse_response(data, "json")
end
export currency_exchange_rate

function fx_intraday(from_currency::String, to_currency::String, interval::String="1min"; outputsize::String="compact", datatype::String="json")
    @argcheck in(interval, ["1min", "5min", "15min", "30min", "60min"])
    @argcheck in(outputsize, ["compact", "full"])
    @argcheck in(datatype, ["json", "csv"])

    uri = _form_uri_head("FX_INTRADAY") * "&from_symbol=$from_currency&to_symbol=$to_currency&interval=$interval" * _form_uri_tail(outputsize, datatype)
    data = _get_request(uri)
    return _parse_response(data, datatype)
end
export fx_intraday

for func in (:daily, :weekly, :monthly)
    x = "fx_$(func)"
    fname = Symbol(x)
    @eval begin
        function ($fname)(from_currency::String, to_currency::String; datatype::String="json")
            @argcheck in(datatype, ["json", "csv"])
            uri = _form_uri_head(uppercase($x)) * "&from_symbol=$from_currency&to_symbol=$to_currency&datatype=$datatype&apikey=" * ENV["ALPHA_VANTAGE_API_KEY"]
            data = _get_request(uri)
            return _parse_response(data, datatype)
        end
        export $fname
    end
end


text:
abstract type AbstractScheduler end

function workerindices_from_threadids(threadids::AbstractVector{<:Integer})
    workerindices = zeros(Int, Threads.nthreads())
    for (wi, tid) in pairs(threadids)
        workerindices[tid] = wi
    end
    return workerindices
end

Base.isopen(scheduler::AbstractScheduler) = !@atomic scheduler.closed
function Base.close(scheduler::AbstractScheduler)
    @record(:close, scheduler)
    @atomic scheduler.closed = true
    wakeall!(scheduler)
    return
end

function scheduler_task()
    thunk = thunkof(current_task())
    worker = @atomic :monotonic thunk.worker
    worker::Worker
    return taskof(worker)
end

pushat!(scheduler::AbstractScheduler, _workerid::Integer, task::Task) =
    push!(scheduler, task)

"""
    taskof(thing) -> task::Core.Task

Return a `task` associated with the `thing` (scheduler, worker, `GenericTask`
etc.).
"""
taskof
taskof(task::Task) = task

baremodule WorkerStates
using Base: @enum, UInt8
@enum Kind::UInt8 begin
    ACTIVE
    PREPARE_WAITING
    WAITING
    NOTIFYING
    CREATED
    DONE
end
end  # baremodule TaskStates

mutable struct Worker
    @const task::Task
    @atomic state::WorkerStates.Kind
end

Worker(task::Task) = Worker(task, WorkerStates.CREATED)
taskof(worker::Worker) = worker.task

Threads.threadid(worker::Worker) = Threads.threadid(taskof(worker))

baremodule WaiterStates
using Base: @enum, UInt8
@enum Kind::UInt8 begin
    # Set by waiter:
    INITIAL
    WAITING
    # Set by notifier:
    DONTWAIT
    NOTIFYING
end
end

mutable struct Waiter
    @const task::Task
    next::Union{Waiter,Nothing}
    @atomic state::WaiterStates.Kind
end

Waiter() = Waiter(current_task(), nothing, WaiterStates.INITIAL)

taskof(waiter::Waiter) = waiter.task

baremodule TaskStates
using Base: @enum, UInt8
@enum Kind::UInt8 begin
    DONE
    ERROR
    STARTED
    CREATED
end
end  # baremodule TaskStates

const PriorityInt = Int

as_priority(::Nothing) = zero(PriorityInt)
function as_priority(priority::Integer)
    priority = convert(PriorityInt, priority)
    return max(typemin(PriorityInt) + 1, priority)
end

struct Closed end

mutable struct Thunk{Scheduler<:AbstractScheduler}
    @const f::OpaqueClosure{Tuple{}}
    @const scheduler::Scheduler
    priority::PriorityInt
    @atomic worker::Union{Nothing,Worker}
    result::Any
    @atomic waiter::Union{Nothing,Closed,Waiter}
    state::TaskStates.Kind

    Thunk{Scheduler}(
        @nospecialize(f),
        scheduler::Scheduler,
        priority::Union{Integer,Nothing},
    ) where {Scheduler<:AbstractScheduler} = new{Scheduler}(
        f,
        scheduler,
        as_priority(priority),
        nothing,
        nothing,
        nothing,
        TaskStates.CREATED,
    )
end

thunkof(tasklike) = taskof(tasklike).code::ConcreteThunk
schedulerof(tasklike) = thunkof(tasklike).scheduler

struct GenericTask
    task::Task
end

function generic_task(
    @nospecialize(f),
    scheduler::AbstractScheduler;
    priority::Union{Integer,Nothing} = nothing,
)
    thunk = ConcreteThunk(as_opaqueclosure(f)::OpaqueClosure{Tuple{}}, scheduler, priority)
    task = Task(thunk)
    task.sticky = false  # will be flipped in `(::Thunk)()`
    return GenericTask(task)
end

as_opaqueclosure(@nospecialize(f)) = @opaque () -> f()
as_opaqueclosure(f::OpaqueClosure) = f
# TODO: transform the thunk to `OpaqueClosure` in a frontend macro

GenericTask(waiter::Waiter) = GenericTask(taskof(waiter))

taskof(task::GenericTask) = task.task

stateof(task) = thunkof(task).state
isfinished(state::TaskStates.Kind) = state == TaskStates.DONE || state == TaskStates.ERROR
isfinished(task) = isfinished(stateof(task))

Base.istaskdone(task::GenericTask) = isfinished(task)
Base.istaskfailed(task::GenericTask) = stateof(task) == TaskStates.ERROR

Historic.taskid(task::GenericTask) = Historic.taskid(task.task)
Historic.taskid(waiter::Waiter) = Historic.taskid(waiter.task)
Historic.taskid(worker::Worker) = Historic.taskid(worker.task)

struct ThunkID
    uint::UInt
end

ThunkID(thunk::Thunk) = ThunkID(UInt(pointer_from_objref(thunk)))

function Base.print(io::IO, id::ThunkID)
    print(io, "thunk-", string(id.uint; base = 16, pad = Sys.WORD_SIZE >> 2))
end


text:
function PSI.construct_device!(
    optimization_container::PSI.OptimizationContainer,
    sys::PSY.System,
    model::PSI.DeviceModel{St, BookKeepingwReservationOutage},
    ::Type{S},
) where {St <: PSY.Storage, S <: PM.AbstractPowerModel}
    devices = PSI.get_available_components(St, sys)

    if !PSI.validate_available_devices(St, devices)
        return
    end

    # Variables
    PSI.add_variables!(
        optimization_container,
        PSI.ActivePowerInVariable,
        devices,
        BookKeepingwReservationOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.ActivePowerOutVariable,
        devices,
        BookKeepingwReservationOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.ReactivePowerVariable,
        devices,
        BookKeepingwReservationOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.EnergyVariable,
        devices,
        BookKeepingwReservationOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.ReserveVariable,
        devices,
        BookKeepingwReservationOutage(),
    )

    # Initial Conditions
    PSI.initial_conditions!(
        optimization_container,
        devices,
        BookKeepingwReservationOutage(),
    )

    # Constraints
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ActivePowerOutVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ActivePowerInVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ReactivePowerVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.energy_capacity_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    outage_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.feedforward!(optimization_container, devices, model, PSI.get_feedforward(model))

    # Energy Balanace limits
    PSI.energy_balance_constraint!(
        optimization_container,
        devices,
        model.formulation,
        S,
        PSI.get_feedforward(model),
    )

    return
end

function PSI.construct_device!(
    optimization_container::PSI.OptimizationContainer,
    sys::PSY.System,
    model::PSI.DeviceModel{St, BookKeepingwReservationOutage},
    ::Type{S},
) where {St <: PSY.Storage, S <: PM.AbstractActivePowerModel}
    devices = PSI.get_available_components(St, sys)

    if !PSI.validate_available_devices(St, devices)
        return
    end

    # Variables
    PSI.add_variables!(
        optimization_container,
        PSI.ActivePowerInVariable,
        devices,
        BookKeepingwReservationOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.ActivePowerOutVariable,
        devices,
        BookKeepingwReservationOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.EnergyVariable,
        devices,
        BookKeepingwReservationOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.ReserveVariable,
        devices,
        BookKeepingwReservationOutage(),
    )

    # Initial Conditions
    PSI.initial_conditions!(
        optimization_container,
        devices,
        BookKeepingwReservationOutage(),
    )

    # Constraints
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ActivePowerOutVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ActivePowerInVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.energy_capacity_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    outage_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.feedforward!(optimization_container, devices, model, PSI.get_feedforward(model))

    # Energy Balanace limits
    PSI.energy_balance_constraint!(
        optimization_container,
        devices,
        model.formulation,
        S,
        PSI.get_feedforward(model),
    )

    return
end

function PSI.construct_device!(
    optimization_container::PSI.OptimizationContainer,
    sys::PSY.System,
    model::PSI.DeviceModel{St, EndOfPeriodEnergyTargetOutage},
    ::Type{S},
) where {St <: PSY.Storage, S <: PM.AbstractPowerModel}
    devices = PSI.get_available_components(St, sys)

    if !PSI.validate_available_devices(St, devices)
        return
    end

    # Variables
    PSI.add_variables!(
        optimization_container,
        PSI.ActivePowerInVariable,
        devices,
        EndOfPeriodEnergyTargetOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.ActivePowerOutVariable,
        devices,
        EndOfPeriodEnergyTargetOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.ReactivePowerVariable,
        devices,
        EndOfPeriodEnergyTargetOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.EnergyVariable,
        devices,
        EndOfPeriodEnergyTargetOutage(),
    )

    # Initial Conditions
    PSI.initial_conditions!(
        optimization_container,
        devices,
        EndOfPeriodEnergyTargetOutage(),
    )

    # Constraints
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ActivePowerOutVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ActivePowerInVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ReactivePowerVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.energy_capacity_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    outage_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.feedforward!(optimization_container, devices, model, PSI.get_feedforward(model))

    # Energy Balanace limits
    PSI.energy_balance_constraint!(
        optimization_container,
        devices,
        model.formulation,
        S,
        PSI.get_feedforward(model),
    )

    # Cost Function
    PSI.cost_function!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )

    return
end

function PSI.construct_device!(
    optimization_container::PSI.OptimizationContainer,
    sys::PSY.System,
    model::PSI.DeviceModel{St, EndOfPeriodEnergyTargetOutage},
    ::Type{S},
) where {St <: PSY.Storage, S <: PM.AbstractActivePowerModel}
    devices = PSI.get_available_components(St, sys)

    if !PSI.validate_available_devices(St, devices)
        return
    end

    # Variables
    PSI.add_variables!(
        optimization_container,
        PSI.ActivePowerInVariable,
        devices,
        EndOfPeriodEnergyTargetOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.ActivePowerOutVariable,
        devices,
        EndOfPeriodEnergyTargetOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.EnergyVariable,
        devices,
        EndOfPeriodEnergyTargetOutage(),
    )

    # Initial Conditions
    PSI.initial_conditions!(
        optimization_container,
        devices,
        EndOfPeriodEnergyTargetOutage(),
    )

    # Constraints
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ActivePowerOutVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ActivePowerInVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.energy_capacity_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    outage_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.feedforward!(optimization_container, devices, model, PSI.get_feedforward(model))

    # Energy Balanace limits
    PSI.energy_balance_constraint!(
        optimization_container,
        devices,
        model.formulation,
        S,
        PSI.get_feedforward(model),
    )

    # Cost Function
    PSI.cost_function!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )

    return
end