text:
export IQNLearner, ImplicitQuantileNet

"""
    ImplicitQuantileNet(;ψ, ϕ, header)

```
        quantiles (n_action, n_quantiles, batch_size)
           ↑
         header
           ↑
feature ↱  ⨀   ↰ transformed embedding
       ψ       ϕ
       ↑       ↑
       s        τ
```
"""
Base.@kwdef struct ImplicitQuantileNet{A,B,C}
    ψ::A
    ϕ::B
    header::C
end

Flux.@functor ImplicitQuantileNet

function (net::ImplicitQuantileNet)(s, emb)
    features = net.ψ(s)  # (n_feature, batch_size)
    emb_aligned = net.ϕ(emb)  # (n_feature, N * batch_size)
    merged =
        Flux.unsqueeze(features, 2) .*
        reshape(emb_aligned, size(features, 1), :, size(features, 2))  # (n_feature, N, batch_size)
    quantiles = net.header(flatten_batch(merged))
    reshape(quantiles, :, size(merged, 2), size(merged, 3))  # (n_action, N, batch_size)
end

"""
    IQNLearner(;kwargs)

See [paper](https://arxiv.org/abs/1806.06923)

# Keyword arguments
- `approximator`, a [`ImplicitQuantileNet`](@ref)
- `target_approximator`, a [`ImplicitQuantileNet`](@ref), must have the same structure as `approximator`
- `κ = 1.0f0`,
- `N = 32`,
- `N′ = 32`,
- `Nₑₘ = 64`,
- `K = 32`,
- `γ = 0.99f0`,
- `stack_size = 4`,
- `batch_size = 32`,
- `update_horizon = 1`,
- `min_replay_history = 20000`,
- `update_freq = 4`,
- `target_update_freq = 8000`,
- `update_step = 0`,
- `default_priority = 1.0f2`,
- `β_priority = 0.5f0`,
- `rng = Random.GLOBAL_RNG`,
- `device_seed = nothing`,
"""
mutable struct IQNLearner{A,T,R,D} <: Any
    approximator::A
    target_approximator::T
    sampler::NStepBatchSampler
    κ::Float32
    N::Int
    N′::Int
    Nₑₘ::Int
    K::Int
    min_replay_history::Int
    update_freq::Int
    target_update_freq::Int
    update_step::Int
    default_priority::Float32
    β_priority::Float32
    rng::R
    device_rng::D
    loss::Float32
end

Functors.functor(x::IQNLearner) =
    (Z = x.approximator, Zₜ = x.target_approximator, device_rng = x.device_rng),
    y -> begin
        x = @set x.approximator = y.Z
        x = @set x.target_approximator = y.Zₜ
        x = @set x.device_rng = y.device_rng
        x
    end

function IQNLearner(;
    approximator,
    target_approximator,
    κ = 1.0f0,
    N = 32,
    N′ = 32,
    Nₑₘ = 64,
    K = 32,
    γ = 0.99f0,
    stack_size = 4,
    batch_size = 32,
    update_horizon = 1,
    min_replay_history = 20000,
    update_freq = 4,
    target_update_freq = 8000,
    update_step = 0,
    default_priority = 1.0f2,
    β_priority = 0.5f0,
    rng = Random.GLOBAL_RNG,
    device_rng = CUDA.CURAND.RNG(),
    traces = SARTS,
    loss = 0.0f0,
)
    copyto!(approximator, target_approximator)  # force sync
    if device(approximator) !== device(device_rng)
        throw(
            ArgumentError(
                "device of `approximator` doesn't match the device of `device_rng`: $(device(approximator)) !== $(device_rng)",
            ),
        )
    end
    sampler = NStepBatchSampler{traces}(;
        γ = γ,
        n = update_horizon,
        stack_size = stack_size,
        batch_size = batch_size,
    )
    IQNLearner(
        approximator,
        target_approximator,
        sampler,
        κ,
        N,
        N′,
        Nₑₘ,
        K,
        min_replay_history,
        update_freq,
        target_update_freq,
        update_step,
        default_priority,
        β_priority,
        rng,
        device_rng,
        loss,
    )
end

function (learner::IQNLearner)(env)
    s = send_to_device(device(learner), state(env))
    s = Flux.unsqueeze(s, ndims(s) + 1)
    τ = rand(learner.device_rng, Float32, learner.K, 1)
    τₑₘ = embed(τ, learner.Nₑₘ)
    quantiles = learner.approximator(s, τₑₘ)
    vec(mean(quantiles; dims = 2)) |> send_to_host
end

embed(x, Nₑₘ) = cos.(Float32(π) .* (1:Nₑₘ) .* reshape(x, 1, :))

function RLBase.update!(learner::IQNLearner, batch::NamedTuple)
    Z = learner.approximator
    Zₜ = learner.target_approximator
    N = learner.N
    N′ = learner.N′
    Nₑₘ = learner.Nₑₘ
    κ = learner.κ
    β = learner.β_priority
    batch_size = learner.sampler.batch_size

    D = device(Z)
    s, r, t, s′ =
        (send_to_device(D, batch[x]) for x in (:state, :reward, :terminal, :next_state))

    τ′ = rand(learner.device_rng, Float32, N′, batch_size)  # TODO: support β distribution
    τₑₘ′ = embed(τ′, Nₑₘ)
    zₜ = Zₜ(s′, τₑₘ′)
    avg_zₜ = mean(zₜ, dims = 2)

    if haskey(batch, :next_legal_actions_mask)
        masked_value = fill(typemin(Float32), size(batch.next_legal_actions_mask))
        masked_value[batch.next_legal_actions_mask] .= 0
        avg_zₜ .+= send_to_device(D, masked_value)
    end

    aₜ = argmax(avg_zₜ, dims = 1)
    aₜ = aₜ .+ typeof(aₜ)(CartesianIndices((0:0, 0:N′-1, 0:0)))
    qₜ = reshape(zₜ[aₜ], :, batch_size)
    target =
        reshape(r, 1, batch_size) .+
        learner.sampler.γ * reshape(1 .- t, 1, batch_size) .* qₜ  # reshape to allow broadcast

    τ = rand(learner.device_rng, Float32, N, batch_size)
    τₑₘ = embed(τ, Nₑₘ)
    a = CartesianIndex.(repeat(batch.action, inner = N), 1:(N*batch_size))

    is_use_PER = haskey(batch, :priority)  # is use Prioritized Experience Replay
    if is_use_PER
        updated_priorities = Vector{Float32}(undef, batch_size)
        weights = 1.0f0 ./ ((batch.priority .+ 1.0f-10) .^ β)
        weights ./= maximum(weights)
        weights = send_to_device(D, weights)
    end

    gs = Zygote.gradient(Flux.params(Z)) do
        z = flatten_batch(Z(s, τₑₘ))
        q = z[a]

        TD_error = reshape(target, N′, 1, batch_size) .- reshape(q, 1, N, batch_size)
        # can't apply huber_loss in RLCore directly here
        abs_error = abs.(TD_error)
        quadratic = min.(abs_error, κ)
        linear = abs_error .- quadratic
        huber_loss = 0.5f0 .* quadratic .* quadratic .+ κ .* linear

        # dropgrad
        raw_loss =
            abs.(reshape(τ, 1, N, batch_size) .- Zygote.dropgrad(TD_error .< 0)) .*
            huber_loss ./ κ
        loss_per_quantile = reshape(sum(raw_loss; dims = 1), N, batch_size)
        loss_per_element = mean(loss_per_quantile; dims = 1)  # use as priorities
        loss =
            is_use_PER ? dot(vec(weights), vec(loss_per_element)) * 1 // batch_size :
            mean(loss_per_element)
        ignore() do
            # @assert all(loss_per_element .>= 0)
            is_use_PER && (
                updated_priorities .=
                    send_to_host(vec((loss_per_element .+ 1.0f-10) .^ β))
            )
            learner.loss = loss
        end
        loss
    end

    update!(Z, gs)

    is_use_PER ? updated_priorities : nothing
end


text:
using Test
using Example2018

# Polynomials form a ring
# 1. addition is associative
# 2. additive identity
# 3. additive inverse
# 4. addition commues
# 5. multiplication is...

@testset "Polynomials form a ring" begin

	z=Polynomial([0])

	for n in 1:100

	      	deg=rand(0:100)

		p=Polynomial(rand(-1000000:1000000,deg+1))

		@test p+z==p
		@test z+p==p
		@test p+(-p)==z
	end
end

ext: jl


text:
"""
    LuaLaTeX

The `lualatex` engine.

Command line options `--output-format=dvi`, `--interaction=batchmode`, and `--halt-on-error` are passed to the `lualatex` executable, and both the input and output directories are also handled. If you would like to pass additional command line arguments, use the [`latexsvg`](@ref) function.
"""
struct LuaLaTeX <: LaTeXEngine end

function runlatex(
    ::LuaLaTeX,
    input_file::AbstractString,
    output_path::AbstractString;
    extra_args::Vector{String}=String[],
    runcommand::Bool=true
)
    command = `lualatex $extra_args --output-format=dvi --output-directory="$output_path" --interaction=batchmode --halt-on-error "$input_file"`
    if runcommand
        run(command)
    end
    return command
end

dviext(::LuaLaTeX) = "dvi"


text:
# Difference between SVGDraw and Luxor:
#   1. svgdraw does not depend on cairo or any other rendering library
#   2. svgdraw can do incremental drawing (in a standalone figure window)
#   3. svgdraw is more oriented toward math drawing, less toward artistic drawing
#   4. svgdraw provides a MetaPost-inspired interface
#
# Warning! we use canvas-coordinate system which flips the y-axis of
# the canonical cartesion-coordinate system
#using PrettyPrint        # for debug
import Base: convert, +, -
import Base: setproperty!, getproperty

# settings
const allow_outer_stylesheet = false

mutable struct R²
    x::Float64
    y::Float64
end

mutable struct R³
    x::Float64
    y::Float64
    z::Float64
end

mutable struct R⁴
    x::Float64
    y::Float64
    w::Float64
    h::Float64
end
BoundingBox = R⁴

convert(::Type{R²}, x::Tuple{Real, Real}) = R²(Float64(x[1]), Float64(x[2]))
convert(::Type{R³}, x::Tuple{Real, Real, Real}) = R³(Float64.(x)...)

+(a::R², b::R²) = R²(a.x + b.x, a.y + b.y)
+(a::R², b::Tuple{Real, Real}) = a + convert(R², b)
+(a::R³, b::R³) = R³(a.x + b.x, a.y + b.y, a.z + b.z)

-(a::R², b::R²) = R²(a.x - b.x, a.y - b.y)
-(a::R³, b::R³) = R³(a.x - b.x, a.y - b.y)

mutable struct Attribute
    id::String
    class::String
    label::String
    isUpdated::Bool     # if the position is up-to-date
end
Attribute(id::String, class::String, label::String) = Attribute(id, class, label, true)

mutable struct Rect <: VGElement
    x::Float64
    y::Float64
    w::Float64
    h::Float64
    a::Union{Attribute, Nothing}
end

function Rect(x::Real, y::Real, w::Real, h::Real; id="", class="", label="")
    x, y, w, h = Float64.((x, y, w, h))
    return Rect(x, y, w, h, Attribute(id, class, label))
end

mutable struct Circ <: VGElement
    x::Float64
    y::Float64
    r::Float64
    a::Union{Attribute, Nothing}
end

function Circ(x::Real, y::Real, r::Real; id="", class="", label="")
    x, y, r = Float64.((x, y, r))
    return Circ(x, y, r, Attribute(id, class, label))
end

function Circ(c::R², r::Number; id="", class="", label="")
    x, y = Float64.((c.x, c.y))
    return Circ(x, y, r, Attribute(id, class, label))
end

# TODO
mutable struct Point <: VGElement
    x::Float64
    y::Float64
    a::Union{Attribute, Nothing}
end

mutable struct Polyline <: VGElement
    points::Vector{R²}
    a::Union{Attribute, Nothing}
end

function Polyline(pt1::R², points::R²...; id="", class="", label="")
    pts = [pt1]
    append!(pts, points...)
    return Polyine(pts, Attribute(id, class, label))
end

function Polyline(pt1::Tuple{Real, Real}, points::Tuple{Real,Real}...;
                  id="", class="", label="")
    pts = [convert(R², pt1)]
    append!(pts, convert.(R², points))
    return Polyline(pts, Attribute(id, class, label))
end

# TODO
mutable struct Text <: VGElement
    x::Float64
    y::Float64
    a::Union{Attribute, Nothing}
end

function Text(x::Real, y::Real, w::Real, h::Real;
              id="", class="", label="")
    x, y = Float64.((x, y))
    return Text(x, y, Attribute(id, class, label))
end

mutable struct Brace <: VGElement
    x::Float64
    y::Float64
    w::Float64
    h::Float64
    se::Union{VGElement, Nothing}      # starting element
    ee::Union{VGElement, Nothing}      # ending element
    a::Union{Attribute, Nothing}
end

function Brace(x::Real, y::Real, w::Real, h::Real;
               id="", class="", label="")
    x, y, w, h = Float64.((x, y, w, h))
    return Brace(x, y, w, h, nothing, nothing, Attribute(id, class, label))
end

function Brace(se::VGElement; id="", class="", label="")
    return Brace(NaN, NaN, NaN, NaN, se, nothing, Attribute(id, class, label, false))
end

function Brace(se::VGElement, ee::VGElement;id="", class="", label="")
    return Brace(NaN, NaN, NaN, NaN, se, ee, Attribute(id, class, label, false))
end

function setproperty!(e::VGElement, name::Symbol, x)
    name === :id && (e.a.id = x; return)
    name === :class && (e.a.class = x; return)
    name === :label && (e.a.label = x; return)
    name === :isUpdated && (e.a.isUpdated = x; return)
    return Base.setfield!(e, name, x)
end

function getproperty(e::VGElement, name::Symbol)
    name === :id && return e.a.id
    name === :class && return e.a.class
    name === :label && return e.a.label
    name === :isUpdated && return e.a.isUpdated
    return Base.getfield(e, name)
end

__BoxShaped = Union{BoundingBox, Rect, Brace}
# Helper functions
boundingbox(e::__BoxShaped) = R⁴(e.x, e.y, e.w, e.h)
boundingbox(e::Circ) = R⁴(e.x - e.r, e.y - e.r, 2*e.r, 2*e.r)

# In the future, Rect object could also have their positions depending on other objects.
midtop(e::__BoxShaped) = R²(e.x + e.w/2, e.y)
midbottom(e::__BoxShaped) = R²(e.x + e.w/2, e.y + e.h)
midleft(e::__BoxShaped) = R²(e.x, e.y + e.h/2)
midright(e::__BoxShaped) = R²(e.x + e.w, e.y + e.h/2)
lefttop(e::__BoxShaped) = R²(e.x, e.y)
righttop(e::__BoxShaped) = R²(e.x + e.w, e.y)
leftbottom(e::__BoxShaped) = R²(e.x, e.y + e.h)
rightbottom(e::__BoxShaped) = R²(e.x + e.w, e.y + e.h)

update_box(::Nothing, pvd::PVDict) = nothing
update_box(e::Rect, pvd::PVDict) = nothing
update_box(e::Circ, pvd::PVDict) = nothing

function update_box(e::Brace, pvd::PVDict)
    e.isUpdated && return
    update_box(e.se)
    update_box(e.ee)
    bb = calc_brace_bbox(e, pvd)
    e.x, e.y, e.w, e.h = bb.x, bb.y, bb.w, bb.h
    e.isUpdated = true
end

#TODO: more helper functions for Circ, ...
#TODO: add midtop_xy, midbottom_xy, ...

# SVG primitives
# assume rule strings in style have been normalized (no leading and traling white
# spaces), and have been cleaned (no user-defined properties)
function make_rect(x::Real, y::Real, w::Real, h::Real, pvd::PVDict)
    return "<rect x=\"$x\" y=\"$y\" width=\"$w\" height=\"$h\" style=$(nstringify(pvd)) />"
end

function make_rect(b::Rect, pvd::PVDict)
    x, y, w, h = b.x, b.y, b.w, b.h
    return "<rect x=\"$x\" y=\"$y\" width=\"$w\" height=\"$h\" style=\"$(nstringify(pvd))\" />"
end

function make_circle(x::Real, y::Real, r::Real, pvd::PVDict)
    return "<circle cx=\"$x\" cy=\"$y\" r=\"$r\" style=\"$(nstringify(pvd))\" />"
end

function make_circle(c::Circ, pvd::PVDict)
    x, y, r = c.x, c.y, c.r
    return "<circle cx=\"$x\" cy=\"$y\" r=\"$r\" style=\"$(nstringify(pvd))\" />"
end

function make_label(el::VGElement, pvd::PVDict)
    anchor = get_property(pvd, "--vg-label-anchor")
    m = eachmatch(r"\n", el.label)
    lines = length(collect(m)) + 1
    pt = calc_label_pos(el, pvd, lines)
    textattr = "x=\"$(pt.x)\" y=\"$(pt.y)\" text-anchor=\"$(anchor)\""
    return "<text $(textattr) >"*el.label*"</text>"
end

const anchor_ratio = Dict([
    ("start", 0),
    ("middle", 0.5),
    ("end", 1)
   ])

# TODO: add support for vertical text
function calc_label_pos(el::VGElement, st::PVDict, lines::Int)
    side = get_property(st, "--vg-label-side")
    dist = parse(Float64, get_property(st, "--vg-label-distance"))
    anchor = get_property(st, "--vg-label-anchor")
    xheight = parse(Float64, get_property(st, "--vg-label-font-size"))
    xvshift = parse(Float64, get_property(st, "--vg-inner-center-vshift"))

    # calculate text position based on the `--vg-label-side` and
    #    `--vg-offset` properties.
    # NOTE: the `x`, `y` properties of a <text> element are the coordinates of
    # the bottom-left of the text box.
    pt = R²(0,0)
    bb = boundingbox(el)
    if side == "above"
        offset = anchor_ratio[anchor]*bb.w
        pt = lefttop(bb) + (offset, -dist - (lines-1)*xheight)
    elseif side == "below"
        offset = anchor_ratio[anchor]*bb.w
        pt = leftbottom(bb) + (offset, dist + xheight)
    elseif side == "left"
        pt = lefttop(bb) + (-dist, (bb.h - xheight)/2)
    elseif side == "right"
        pt = righttop(bb) + (dist, (bb.h - xheight)/2)
    else # side == "center"
        #FIXME: can not get the size of the text area, we have to guess...
        # the --vg-label-anchor must be `middle`
        pt = midtop(bb) + (0, (bb.h + xheight)/2 + xvshift)
    end
    return pt
end

calc_brace_bbox(se::VGElement, pvd::PVDict) = calc_brace_bbox(se, nothing, pvd)

function calc_brace_bbox(se::VGElement, ee::Union{VGElement, Nothing}, pvd::PVDict)
    side = get_property(pvd, "--vg-side")
    d = parse(Float64, get_property(pvd, "--vg-distance"))
    off = parse(Float64, get_property(pvd, "--vg-offset")) # relative val given in percentage
    bh = parse(Float64, get_property(pvd, "--vg-brace-height", "5"))

    # ew: the width of the bounding box if --vg-side == above|below
    # eh: the height of the bounding box if --vg-side == above|below
    # ew2: the width of the bounding box if --vg-side == left| right
    # eh2: the height of the bounding box if --vg-side == left|right
    bb = R⁴(0,0,0,0)
    if ee == nothing
        if off >= 0.5
            # degenerated brace (only a line segment)
            ew = 0
            eh2 = 0
        else
            ew = se.w*(1-2*off)
            eh2 = se.h*(1-2*off)
        end
    else
        ew = ee.x + ee.w*(1-off) - se.x - se.w*off
        eh2 = ee.y + ee.h*(1-off) - se.y - se.h*off
    end
    eh, ew2 = bh, bh

    if side == "above"
        bb = R⁴(se.x + off*se.w, se.y - d - eh, ew, eh)
    elseif side == "below"
        bb = R⁴(se.x + off*se.w, ee.y + ee.h + d, ew, eh)
    elseif side == "left"
        bb = R⁴(se.x - d - ew2, se.y + off*se.h, ew2, eh2)
    else # side == "right"
        bb = R⁴(ee.x + ee.w + d, se.y + off*se.h, ew2, eh2)
    end

    return bb
end

# Note: there can be a chain-reaction: the position of an overbrace label depends
# on the position of the overbrace, which in turn depends on the position of
# the boxes that the overbrace encloses. So the overbrace position and size
# should be calculated before calculating those of its label. But the
# quesiton is, how to distinguish a standalone overbrace and an pos-size updated
# overbrace?
function make_svg_root(eles...; css=CssRules())
    spvd = cssrules_to_spvdict(css)
    content = ""
    leftxs, rightxs, topys, bottomys = [], [], [], []
    for el in eles
        if el isa Vector
            for e in el
                st = get_style_PVD(e, spvd)
                content *= make_svg_ele(e, st)
                bb = boundingbox(e)
                push!(leftxs, bb.x)
                push!(rightxs, bb.x + bb.w)
                push!(topys, bb.y)
                push!(bottomys, bb.y + bb.h)
                if !isempty(e.label)
                    # FIXME we do not know yet how to get the pixel length of a text
                    # area, thus we can not update maxxs for it
                    content *= make_label(e, st)
                end
            end
        elseif el isa VGElement
            st = get_style_PVD(el, spvd)
            content *= make_svg_ele(el, st)
            bb = boundingbox(el)
            push!(leftxs, bb.x)
            push!(rightxs, bb.x + bb.w)
            push!(topys, bb.y)
            push!(bottomys, bb.y + bb.h)
            if !isempty(e.label)
                content *= make_label(el, st)
            end
        else
            println(typeof(el))
            error("Unknown type in writing SVG elements")
        end
    end

    # calculate viewBox
    lmargin, rmargin, tmargin, bmargin = 40, 10, 60, 30
    vx = string(minimum(leftxs) - lmargin)
    vy = string(minimum(topys) - tmargin)
    vw = string(maximum(rightxs) - minimum(leftxs) + lmargin + rmargin)
    vh = string(maximum(bottomys) - minimum(topys) + tmargin + bmargin)

    return "<svg width=\"$(vw)\" height=\"$(vh)\" viewBox=\"$(vx) $(vy) $(vw) $(vh)\">"*content*"</svg>"
end

function write_svg_ele(el, pvd::PVDict)
    return "NOT implemented yet"
end

function make_svg_ele(el::Rect, pvd::PVDict)
    if !allow_outer_stylesheet
        pvd = merge(supp_properties, pvd)
    end
    make_rect(el, pvd)
end

function make_svg_ele(el::Circ, pvd::PVDict)
    if !allow_outer_stylesheet
        pvd = merge(supp_properties, pvd)
    end
    make_circle(el, pvd)
end

function make_svg_ele(el::Polyline, pvd::PVDict)
    if !allow_outer_stylesheet
        pvd = merge(supp_properties, pvd)
    end
    pts= ""
    for pt in el.points
        pts *= "$(pt.x),$(pt.y) "
    end
    st = nstringify(pvd)
    return "<polyline fill=\"none\" style=\"$st\" points=\"$(pts)\" />"
end

# four cases to deal with:
#   1. lsd == nothing, rsd == nothing: a standalone overbrace
#   2. lsd != nothing, rsd == nothing: an overbrace over one box, degenerated
#   3. lsd != nothing, rsd == nothing: an overbrace over one box, non-degenerated
#   4. lsd != nothing, rsd != nothing: an overbrace over two boxes
#
# use properties: --vg-offset, --vg-distance, --vg-shape, --vg-suppress
function make_svg_ele(el::Brace, pvd::PVDict)
    if !allow_outer_stylesheet
        pvd = merge(supp_properties, pvd)
    end
    bb = BoundingBox(0, 0, 0, 0)
    if  el.se === nothing && el.ee === nothing
        # user-defined position-dimension
        bb = BoundingBox(el.x, el.y, el.w, el.h)
    elseif el.ee === nothing
        bb = calc_brace_bbox(el.se, pvd)
    else
        bb = calc_brace_bbox(el.se, el.ee, pvd)
    end
    el.x, el.y, el.w, el.h = bb.x, bb.y, bb.w, bb.h
    el.isUpdated = true

    side = get_property(pvd, "--vg-side")
    s = ""
    # how about making curly braces?
    if side == "above"
        pline = Polyline((bb.x, bb.y+bb.h), (bb.x, bb.y),
                         (bb.x+bb.w, bb.y), (bb.x+bb.w, bb.y+bb.h))
        s = make_svg_ele(pline, pvd) # get path color, path stroke, ... from pvd
    elseif side == "below"
        pline = Polyline((bb.x, bb.y), (bb.x, bb.y+bb.h),
                         (bb.x+bb.w, bb.y+bb.h), (bb.x+bb.w, bb.y))
        s = make_svg_ele(pline, pvd)
    elseif side == "left"
        pline = Polyline((bb.x+bb.w, bb.y), (bb.x, bb.y),
                         (bb.x, bb.y+bb.h), (bb.x+bb.w, bb.y+bb.h))
        s = make_svg_ele(pline, pvd)
    else # side == "right"
        s = pline = Polyline((bb.x, bb.y), (bb.x+bb.w, bb.y),
                         (bb.x+bb.w, bb.y+bb.h), (bb.x, bb.y+bb.h))
        make_svg_ele(pline, pvd)
    end
    return s
end


text:
@testset "Arc-Hybrid" begin

    @test showstr(LeftArc()) == "LeftArc()"
    @test showstr(LeftArc("mod")) == "LeftArc(mod)"
    for t in [LeftArc(), RightArc(), Shift()]
        @test isempty(DependencyTrees.args(t)) && isempty(DependencyTrees.kwargs(t))
    end
    for t in [LeftArc("label"), RightArc("label")]
        @test !isempty(DependencyTrees.args(t)) && isempty(DependencyTrees.kwargs(t))
    end

    AH = ArcHybrid()

    tb = test_treebank("hybridtests.conll")

    trees = collect(tb)
    @test length(trees) == 4
    @test length.(trees) == [6, 6, 9, 18]

    @test DT.projective_only(ArcHybrid())

    t1 = first(trees)
    c1,c2 = initconfig(AH, CoNLLU, [t.form for t in t1]), initconfig(AH, t1)
    @test [t.form for t in tokens(c1)] == [t.form for t in tokens(c2)]

    @testset "Static Oracle" begin
        oracle = StaticOracle(AH, arc=typed)
        model(x) = nothing
        errorcb(x, ŷ, y) = nothing

        function test_oracle(gold)
            gold_xys = collect(xys(oracle, gold))
            cfg, t = last(gold_xys)
            cfg = t(cfg)
            graph = DependencyTree(tokens(cfg))
            @test all(enumerate(graph)) do (i, token)
                g = gold[i]
                token.form == g.form && token.deprel == g.deprel && token.head == g.head
            end
        end

        for tree in tb
            test_oracle(tree)
        end

        s1, s2, s3, s4 = collect(tb)

        @test last.(xys(oracle, s1)) == [Shift(), LeftArc("nsubj"), Shift(),
                                         Shift(), RightArc("dobj"), Shift(),
                                         LeftArc("case"), Shift(), RightArc("inst"),
                                         Shift(), RightArc("pu"), RightArc("root")]

        @test last.(xys(oracle, s2)) == [Shift(), LeftArc("nsubj"), Shift(), Shift(),
                                         Shift(), Shift(), RightArc("case"),
                                         RightArc("nmod"), RightArc("dobj"), Shift(),
                                         RightArc("pu"), RightArc("root")]

        @test last.(xys(oracle, s3)) == [Shift(), LeftArc("att"), Shift(), LeftArc("sbj"),
                                         Shift(), Shift(), LeftArc("att"), Shift(), Shift(),
                                         Shift(), LeftArc("att"), Shift(), RightArc("pc"),
                                         RightArc("att"), RightArc("obj"), Shift(),
                                         RightArc("pu"), RightArc("pred")]

        @test last.(xys(oracle, s4)) == [Shift(), LeftArc("NNP"), Shift(), Shift(),
                                         RightArc("P"), Shift(), LeftArc("CD"), Shift(),
                                         LeftArc("NNS"), Shift(), RightArc("JJ"), Shift(),
                                         RightArc("PU"), LeftArc("NNP"), Shift(), Shift(),
                                         Shift(), LeftArc("DT"), Shift(), RightArc("NN"),
                                         Shift(), Shift(), Shift(), LeftArc("JJ"),
                                         LeftArc("DT"), Shift(), RightArc("NN"),
                                         RightArc("IN"), Shift(), Shift(), RightArc("CD"),
                                         RightArc("NNP"), RightArc("VB"), Shift(),
                                         RightArc("P"), RightArc("MD")]
    end

    @testset "Dynamic Oracle" begin
        TS = Union{LeftArc, RightArc, Shift}
        oracle = DynamicOracle(ArcHybrid(), arc=typed)
        model(x) = Shift()
        function errorcb(x, ŷ, y)
            @test typeof(x) <: DT.Parse.ArcHybridConfig && typeof(y) <: TS
            @test ŷ == Shift()
            @test typeof(y) <: Union{Shift, LeftArc, RightArc}
        end
    end

    @test showstr(Shift()) == "Shift()"
    @test showstr(Reduce()) == "Reduce()"
    @test showstr(LeftArc()) == "LeftArc()"
    @test showstr(LeftArc("nsubj")) == "LeftArc(nsubj)"
    @test showstr(RightArc()) == "RightArc()"
    @test showstr(RightArc("nsubj")) == "RightArc(nsubj)"
    @test showstr(LeftArc("dobj", upos="NN")) == "LeftArc(dobj; upos=NN)"
end


text:
function get_a_gf0(obj::SBASOrbitalElements)
    return jcall(obj, "getAGf0", jdouble, ())
end

function get_a_gf1(obj::SBASOrbitalElements)
    return jcall(obj, "getAGf1", jdouble, ())
end

function get_iodn(obj::SBASOrbitalElements)
    return jcall(obj, "getIODN", jint, ())
end

function get_prn(obj::SBASOrbitalElements)
    return jcall(obj, "getPRN", jint, ())
end

function get_time(obj::SBASOrbitalElements)
    return jcall(obj, "getTime", jdouble, ())
end

function get_toc(obj::SBASOrbitalElements)
    return jcall(obj, "getToc", jdouble, ())
end

function get_week(obj::SBASOrbitalElements)
    return jcall(obj, "getWeek", jint, ())
end

function get_x(obj::SBASOrbitalElements)
    return jcall(obj, "getX", jdouble, ())
end

function get_x_dot(obj::SBASOrbitalElements)
    return jcall(obj, "getXDot", jdouble, ())
end

function get_x_dot_dot(obj::SBASOrbitalElements)
    return jcall(obj, "getXDotDot", jdouble, ())
end

function get_y(obj::SBASOrbitalElements)
    return jcall(obj, "getY", jdouble, ())
end

function get_y_dot(obj::SBASOrbitalElements)
    return jcall(obj, "getYDot", jdouble, ())
end

function get_y_dot_dot(obj::SBASOrbitalElements)
    return jcall(obj, "getYDotDot", jdouble, ())
end

function get_z(obj::SBASOrbitalElements)
    return jcall(obj, "getZ", jdouble, ())
end

function get_z_dot(obj::SBASOrbitalElements)
    return jcall(obj, "getZDot", jdouble, ())
end

function get_z_dot_dot(obj::SBASOrbitalElements)
    return jcall(obj, "getZDotDot", jdouble, ())
end


text:
#Demo 9: Integration/differentiation
#-------------------------------------------------------------------------------

using CMDimCircuits
CMDimCircuits.@using_CData()

#Get a demo display:
include(CMDimCircuits.demoplotcfgscript)
#Normally use something like:
#CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
vvst = cons(:a, labels = set(yaxis="Amplitude (V)", xaxis="Time (s)"))
sigvscosref = cons(:a, labels = set(yaxis="Signal Amplitude (V)", xaxis="Ref cos(t) Value"))
noline = cons(:a, line = set(style=:none))

#==Input data
===============================================================================#
osr = 50 #(fundamental) oversampling ratio
ncycles = 5

#==Computations
===============================================================================#
t = DataF1(0:(2pi/osr):(ncycles*(2pi)))
y = cos(t)
dydt = deriv(y)
∫ydt = iinteg(y)
x2x=xcross1(y - ∫ydt, n=2)
x2y=ycross1(y, ∫ydt, n=2)
println("crossing#2: ($x2x, $x2y)")

#Generate parameter sweeps:
sweeplist = PSweep[
	PSweep("phi", pi*[0, .25, .5, .75])
]

𝜙 = parameter(DataHR, sweeplist, "phi")
ydel = cos(t+𝜙)
yref = y

#==Generate plot
===============================================================================#
plot1 = push!(cons(:plot, vvst, title="Integration/Differentiation"),
	cons(:wfrm, y, label="sig"),
	cons(:wfrm, dydt, label="dy/dt"),
	cons(:wfrm, ∫ydt, label="integ(y)dt"),
	cons(:wfrm, ycross(y, dydt), noline, glyph=set(shape=:o, color=:red), label="crossings"),
)
plot2 = push!(cons(:plot, sigvscosref, title="Lissajous Curves: yvsx(cos(t+ϕ), cos(t))"),
	cons(:wfrm, yvsx(ydel,yref)),
)

pcoll = push!(cons(:plotcoll, title="Integration/Differentiation/yvsx"), plot1, plot2)
	pcoll.ncolumns = 1

#==Display results in pcoll
===============================================================================#
EasyPlot.displaygui(pcoll)

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll


text:
using Test, FormulaOneMethod

#@testset "Testing in a diffusion system" begin

    using LinearAlgebra, SparseArrays, SuiteSparse, SciMLBase, FormulaOneMethod

    # 2D Laplacian
    function laplacian_2D(nx, ny, k)
        # the 4 neighbours in 2D, in Cartesian coordinates
        Ns = [CartesianIndex(ntuple(x -> x == d ? 1 : 0, 2)) for d in 1:2]
        neighbors = [Ns; -Ns]
        # Laplacian
        Δ = sparse([], [], Vector{Float64}(), nx*ny, nx*ny)
        R = CartesianIndices((nx,ny))
        # Fill the Laplacian within the borders
        for d in 1:2
            N = Ns[d]
            # R′ is the range of indices without the borders in dimension `d`
            R′ = CartesianIndices((nx,ny) .- 2 .* N.I)
            R′ = [r + N for r in R′]
            # Convert to linear indices to build the Laplacian
            u = vec(LinearIndices((nx,ny))[R′])
            in = LinearIndices((nx,ny))[first(R) + N] - LinearIndices((nx,ny))[first(R)]
            # Build the Laplacian (not the fastest way but easier-to-read code)
            Δ += sparse(u, u      , -2k[d], nx*ny, nx*ny)
            Δ += sparse(u, u .- in, +k[d], nx*ny, nx*ny)
            Δ += sparse(u, u .+ in, +k[d], nx*ny, nx*ny)
        end
        return Δ
    end

    # Create a small diffusion problem
    nx = 50
    ny = 40
    n = nx * ny
    m = 4 # parameters

    T(p) = -laplacian_2D(nx ,ny, [p[1], p[2]])

    # TODO make this cleaner, something with a comprehension list maybe
    # A local source (fixing the source at the center to be p[3]
    source_index = LinearIndices((nx,ny))[Int(round(nx/2)),Int(round(ny/2))]
    source_vec = zeros(n)
    source_vec[source_index] = 1
    source(x,p) = source_vec .* (p[3] .- x)
    ∇source(x,p) = -sparse(Diagonal(source_vec))

    # A global sink
    sink(x,p) = p[4] * x / n
    ∇sink(x,p) = p[4] * I / n

    # Define state function F(x,p) and Jacobian ∇ₓF(x,p)
    F(x,p) = ODEFunction((x,p) -> -T(p) * x + source(x,p) - sink(x,p),
                         jac = (x,p) -> -T(p) + ∇source(x,p) - ∇sink(x,p))

    # Basic Newton solver
    function newton_solve(F, x; Ftol=1e-10)
        while norm(F(x)) ≥ Ftol
            x .-= ∇ₓF(x) \ F(x)
        end
        return x
    end

    # Create a type for the solver's algorithm
    struct MyAlg <: SciMLBase.AbstractSteadyStateAlgorithm end

    # Overload SciMLBase's solve function
    function SciMLBase.solve(prob::SciMLBase.AbstractSteadyStateProblem,
                              alg::MyAlg;
                              Ftol=1e-10)
        # Define the functions according to SciMLBase.SteadyStateProblem type
        p = prob.p
        t = 0
        x0 = copy(prob.u0)
        dx, df = copy(x0), copy(x0)
        F(x) = prob.f(dx, x, p, t)
        ∇ₓF(x) = prob.f(df, dx, x, p, t)
        # Compute `u_steady` and `resid` as per SciMLBase using my algorithm
        x_steady = newton_solve(F, x0, Ftol=Ftol)
        resid = F(x_steady)
        # Return the common SciMLBase solution type
        SciMLBase.build_solution(prob, alg, x_steady, resid; retcode=:Success)
    end

    # Overload SciMLBase's SteadyStateProblem constructor
    function SciMLBase.SteadyStateProblem(F, x, p)
        f(dx, x, p, t) = F(x, p)
        f(df, dx, x, p, t) = ∇ₓF(x, p)
        return SciMLBase.SteadyStateProblem(f, x, p)
    end

    # Define objective function f(x,p) and ∇ₓf(x,p)
    function state_mismatch(x)
        δ(x) = x - ones(n)
        return 0.5δ(x)'δ(x)
    end
    function parameter_mismatch(p)
        δ(p) = log.(p)
        return 0.5δ(p)'δ(p)
    end
    f(x,p) = state_mismatch(x) + parameter_mismatch(p)
    function ∇ₓf(x,p)
        δ(x) = x - ones(n)
        return δ(x)'
    end

    # TODO test the derivatives are correct!
    F1 = FormulaOneMethod
    # Initialize the cache for storing reusable objects
    x₀ = rand(n)
    p₀ = rand(m)
    mem = F1.initialize_mem(x₀, p₀)

    # Compute the objective function, 𝑓̂(𝒑)
    objective(p) = F1.f̂(f, F, mem, p, MyAlg())
    objective(p₀)

    # Compute the gradient, ∇𝑓̂(𝒑)
    gradient(p) = F1.∇f̂(f, F, ∇ₓf, mem, p, MyAlg())
    gradient(p₀)

    # Compute the Hessian matrix, ∇²𝑓̂(𝒑)
    Hessian(p) = F1.∇²f̂(f, F, ∇ₓf, mem, p, MyAlg())
    Hessian(p₀)

#end


text:
module Life

export calls, query, completed, instances, processes
export init_life, construct_life_iteration

using Instance.Monitoring
using Master.Types
using Master.Startup_utils
using Master.Instance_utils

calls = Array{Call_info,1}() # All remote functions calls
query = Array{Requested_task,1}() # All requests to run
completed = Array{Completed_task,1}() # All completed requests
instances = Array{Instance_info,1}() # All available machines
processes = Array{Process_info,1}() # All running processes.
# The same machine will share the one process for all gpu instances.
process_table = Dict() # The table that maps hosts to process identifiers.
# The table that maps hosts to process identifiers.

function init_life()
    global query, instances, process_table
    query = read_tasks_list("tasks.csv")
    instances = read_instance_list("instances.csv")
    process_table = connect_machines!(processes, instances)
end

function connect_machines!(processes::Array{Process_info,1},
                           instances::Array{Instance_info,1},
                           remote_pwd::AbstractString="~/julia_temp",
                           remote_exe::AbstractString="~/julia/bin/julia")::Dict
    process_table = Dict()
    for instance in instances
        if (instance.userhost in keys(process_table)) == false
            proc = addprocs([instance.userhost],
                            tunnel=true,
                            dir=remote_pwd,
                            exename=remote_exe,
                            topology=:master_slave)
            push!(processes, Process_info(instance.userhost, proc[1]))
            process_table[instance.userhost] = proc[1]
        end
        instance.process_id = process_table[instance.userhost]
    end
    return process_table
end

function construct_life_iteration(instance_error_limit,
                                  logs_update_timeout_minutes,
                                  process_runtime_minutes)::Function
    iteration = 1
    logs_update_timeout = logs_update_timeout_minutes * 60 * 1000
    process_runtime = process_runtime_minutes * 60 * 1000 # from minutes to ms
    function life_iteration()
        for instance in instances
            if instance.status == "OK"
                try
                    # gpu_state = remotecall_fetch(get_gpu_state,
                    #         instance.process_id, instance.gpu_id)
                    gpu_state = 1
                    if gpu_state == 1 && length(query) > 0
                        request = pop!(query)
                        kill_python3(instance)
                        call!(calls, request, instance)
                        instance.status = "BUSY"
                    else
                        # It indicates it is occupied by the other program.
                        uh = instance.userhost
                        warn("Instance $(uh) is occupied by unknown force.")
                        instance.status = "OK"
                    end
                catch e
                    userhost = instance.userhost
                    warn("Failed to check the state of instance $(userhost).")
                    println(e)
                    instance.error_counter += 1
                    instance.status = "OK"
                    if instance.error_counter >= instance_error_limit
                        instance.status = "ERROR"
                    end
                end
            end
        end
        for call in calls
            if call.active == true
                try
                    pos_before = call.log_pos
                    get_log!(call)
                    pos_after = call.log_pos
                    if pos_after > pos_before
                        call.log_last_update = now()
                    end
                    if (now() - call.log_start).value > process_runtime
                        info("Call runtime reached the time limit.")
                        error("Call runtime reached the time limit.")
                    end
                    if (now() - call.log_last_update).value >logs_update_timeout
                        info("Call failed to update the logs in time.")
                        error("Call failed to update the logs in time.")
                    end
                    if call.process.termsignal >= 0
                        info("Call exited.")
                        error("Call exited.")
                    end
                    logfile = open(string(split(
                        split(string(call.log_file), ' ')[2], '>')[1]))
                    log_is_fine = startswith(readlines(logfile)[end],
                                             "Training finished.")
                    if log_is_fine
                        info("Logs showed that the call was finished.")
                        error("Logs showed that the call was finished.")
                    end
                catch e
                    println(e)
                    call.active = false
                    log_is_fine = false
                    try
                        logfile = open(string(split(
                                split(string(call.log_file), ' ')[2], '>')[1]))
                        log_is_fine = startswith(readlines(logfile)[end],
                                                 "Training finished.")
                    catch e
                        println(e)
                        log_is_fine = false
                    end
                    timeout_is_fine =
                            (now() - call.log_start).value > process_runtime
                    if log_is_fine || timeout_is_fine
                        info("The call $(call.task.id) succeed.")
                        push!(completed,
                              Completed_task(call.task, call.instance))
                    else
                        push!(query, call.task)
                        info("The call $(call.task.id) failed.")
                        userhost = call.instance.userhost
                        info("The error counter of instance $(userhost)." *
                             "incremented due to call fail.")
                        call.instance.error_counter += 1
                    end
                    kill_python3(call.instance)
                    kill_python3(call.instance)
                    call.instance.status = "OK"
                    if call.instance.error_counter >= instance_error_limit
                        call.instance.status = "ERROR"
                    end
                end
            end
        end
        println("Iteration $iteration:")
        iteration += 1
    end
    return life_iteration
end

end # module

max_stars_repo_name: Gollor/JuliaAutomation

max_issues_repo_name: Gollor/JuliaAutomation

max_forks_repo_name: Gollor/JuliaAutomation


text:
using Distributions
using ..PhysicalVectors

export RandomDislocation, UniformBurgersVector

@enum DislocationProperties begin
    BurgersVector = 1
    DislocationOrigin
end

"""
Interface type for the joint-distribution for the dislocations.

# Additional Information
For any given dislocation, it has a random Burgers vector and random location in the lattice.
This type should define both additional distributions for the dislocations.
"""
abstract type RandomDislocation end

"""
    system_size( rbv::RandomDislocation )

Interface function for the `RandomDislocation` type.
"""
system_size( rbv::RandomDislocation ) = error("No implementation defined for $(typeof(rbv)).")

"""
    rand_burger_vector( rbv::RandomDislocation )

Interface function for the `RandomDislocation` type.
"""
rand_burger_vector( rbv::RandomDislocation ) = error("No implementation defined for $(typeof(rbv)).")

"""
    rand_dislocation_source( rbv::RandomDislocation )

Interface function for the `RandomDislocation` type.
"""
rand_dislocation_source( rbv::RandomDislocation ) = error("No implementation defined for $(typeof(rbv)).")

"""
    rand_dislocation( rbv::UniformBurgersVector ) -> Tuple{Vector2D, Vector2D}

Generate a random dislocation by defining its Burgers vector and origin. Returns a tuple of `( Burgers_Vector, Source_Location )`.
"""
rand_dislocation( rbv::RandomDislocation ) = ( rand_burger_vector(rbv), rand_dislocation_source(rbv) )

"""
    origin_already_present( idx, origin, all_dislocations )

Determine if the `origin isa Vector2D` is already present within the `all_dislocations` array up until the `idx` element using `Vectors2d.equal`.
"""
function origin_already_present( idx, origin, all_dislocations )
    if idx == 1
        return false
    end
    already_here = false
    for dis_idx ∈ 1:(idx-1)
        already_here = equal( origin, all_dislocations[Int(DislocationOrigin), dis_idx] )
    end
    return already_here
end

"""
    collect_dislocations( rbv::RandomDislocation, num_dislocations ) -> Array{Vector2D, 2, num_dislocations}
"""
function collect_dislocations( rbv::RandomDislocation, num_dislocations )
    all_dislocations = Array{Vector2D{Float64}, 2}(undef, 2, num_dislocations)
    for idx ∈ eachindex(1:num_dislocations)
        (bob, origin) = rand_dislocation( rbv )
        # Be sure that the same origin is not used multiple times
        while origin_already_present( idx, origin, all_dislocations )
            origin = rand_dislocation( rbv )[ Int(DislocationOrigin) ]
        end
        all_dislocations[Int(BurgersVector), idx]     = Vector2D( bob )
        all_dislocations[Int(DislocationOrigin), idx] = Vector2D( origin )
    end
    return all_dislocations
end

"""
Struct containing uniform distribution for the Burger's vectors and their locations

# Struct Members
* `axes::Tuple{Int}`: a tuple containing `axes[1] = Lx` and `axes[2] = Ly`.
* `burgers_vectors::Vector`: the container with possible Burger's vectors
* `random_position::Tuple{Rand_Loc, Rand_Loc}`: a `Tuple` of `Distribution`s for generating random dislocation positions
* `random_burger_vector::Rand_BV`: the `Distribution` for generating random Burgers vectors
"""
struct UniformBurgersVector{Rand_Loc <: Distribution, Rand_BV <: Distribution} <: RandomDislocation
    axes::Tuple{Int, Int}
    burgers_vectors::Vector{Vector2D{Float64}}
    random_position::Tuple{Rand_Loc, Rand_Loc}
    random_burger_vector::Rand_BV
end

"""
    UniformBurgersVector(; Lx, Ly = Lx, burgers_vectors )

Convenient keyword constructor for the `UniformBurgersVector` struct which makes sure the distributions agree with other members.
"""
function UniformBurgersVector(; Lx, Ly = Lx, burgers_vectors )
    axes = (Lx, Ly)
    rand_pos = ( DiscreteUniform( 1, Lx ), DiscreteUniform(1, Ly) )
    rand_bv = DiscreteUniform(1, length(burgers_vectors))
    return UniformBurgersVector( axes,
                                 burgers_vectors,
                                 rand_pos,
                                 rand_bv )
end

"""
    rand_burger_vector( ubv::UniformBurgersVector ) -> Vector2D{Float64}

Returns a random `Vector2D` for a Burger's vector according to the `ubv <: UniformBurgersVector` allowed Burgers vectors.
"""
function rand_burger_vector( ubv::UniformBurgersVector )
    index = rand( ubv.random_burger_vector )
    return ubv.burgers_vectors[index]
end

"""
    rand_dislocation_source( ubv::UniformBurgersVector ) -> Vector2D{Float64}

Returns a random `Vector2D` that denotes the center of a plaquette for a dislocation origin.
"""
rand_dislocation_source( ubv::UniformBurgersVector ) = Vector2D( ( 0.5 + rand(ubv.random_position[1]), 0.5 + rand(ubv.random_position[2]) ) )

"""
    system_size( ubv::RandomDislocationDistribution )

Return the `prod`uct over the `axes` member `Tuple`.
"""
system_size( ubv::UniformBurgersVector ) = prod(ubv.axes)


text:
"""
    BootstrapFilter(model, n, [rng])

Construct a standard bootstrap particle filter.

The Bootstrap filter was first described in Gordon, N. J., Salmond, D. J., & Smith, A. F. M. "Novel approach to nonlinear / non-Gaussian Bayesian state estimation", with the added robustness of the LowVarianceResampler.

# Arguments
- `model`: a model for the prediction dynamics and likelihood reweighing, for example a `POMDP` or `ParticleFilterModel`
- `n::Integer`: number of particles
- `rng::AbstractRNG`: random number generator

For a more flexible particle filter structure see [`BasicParticleFilter`](@ref).
"""
function BootstrapFilter(model, n::Int, rng::AbstractRNG=Random.GLOBAL_RNG)
    return BasicParticleFilter(model, LowVarianceResampler(n), n, rng)
end


text:
using CairoMakie
using CSV
using Colors

color_lst = [colorant"#332288" colorant"#555555" colorant"#117733" colorant"#882255"]
# Blue, Grey, Green, Red, Yellow
# Colorblind friendly color palette

f = CSV.File("StaticSampling/StaticSampling.csv")

static_random_sampling_trace = f.random_sampling_trace
static_optimal_sampling_trace = f.optimal_sampling_trace

f = CSV.File("LinearSampling/LinearSampling.csv")

Linear_random_sampling_trace = f.random_sampling_trace
Linear_optimal_sampling_trace = f.optimal_sampling_trace
Linear_approx_sampling_trace  = f.approx_sampling_trace
samples = f.samples

f = CSV.File("HopfSampling/HopfSampling.csv")

Hopf_random_sampling_trace  = f.random_sampling_trace
Hopf_optimal_sampling_trace = f.optimal_sampling_trace
Hopf_approx_sampling_trace  = f.approx_sampling_trace
Hopf_samples = f.samples

f = CSV.File("LorenzSampling/LorenzSampling.csv")

Lorenz_random_sampling_trace = f.random_sampling_trace
Lorenz_optimal_sampling_trace = f.optimal_sampling_trace
Lorenz_approx_sampling_trace = f.approx_sampling_trace
Lorenz_samples = f.samples

max_steps = length(samples)

noto_sans = "./resources/NotoSans-Regular.ttf"
noto_sans_bold = "./resources/NotoSans-Bold.ttf"

tickfontsize    = 26
labelfontsize   = 30
basewidth = 2

f = Figure(font=noto_sans, resolution=(1200,550), figure_padding=40)

#ax0 = Axis( f[1,1],
#            xticklabelsize=tickfontsize,
#            yticklabelsize=tickfontsize,
#            yticklabelpad=2,
#            xlabel="Sample", ylabel = "Tr(Σ)",
#            xlabelsize=labelfontsize,
#            ylabelsize=labelfontsize,
#            yscale = log10,
#            yminorticksvisible = true,
#            yminorgridvisible = true,
#            yminorticks = IntervalsBetween(10),
#            xminorticksvisible = true,
#            xminorgridvisible = true,
#            xminorticks = IntervalsBetween(10),
#            title="Static Value",
#            titlesize=labelfontsize)
#
#l2 = lines!(ax0, 1:max_steps, static_random_sampling_trace, color=color_lst[1], linewidth=basewidth+2)
#l1 = lines!(ax0, 1:max_steps, static_optimal_sampling_trace, color=color_lst[3], linewidth=basewidth)
#
#ylims!(ax0,(5e-4,1))
#xlims!(ax0,(0,1000))

ax1 = Axis(  f[1,1],
            xticklabelsize=tickfontsize,
            yticklabelsize=tickfontsize,
            yticklabelpad=2,
            xlabel="Sample", ylabel = "Tr(Σ)",
            xlabelsize=labelfontsize,
            ylabelsize=labelfontsize,
            yscale = log10,
            yminorticksvisible = true,
            yminorgridvisible = true,
            yminorticks = IntervalsBetween(10),
            xminorticksvisible = true,
            xminorgridvisible = true,
            xminorticks = IntervalsBetween(10),
            title="Linear System",
            titlesize=labelfontsize)

ylims!(ax1,(5e-4,1))
xlims!(ax1,(0,1000))

l2 = lines!(ax1, 1:max_steps, Linear_random_sampling_trace, color=color_lst[1],linewidth=basewidth)
l3 = lines!(ax1, 1:max_steps, Linear_random_sampling_trace ./ 2, color=color_lst[2], linewidth=basewidth+2)
l1 = lines!(ax1, 1:max_steps, Linear_optimal_sampling_trace, color=color_lst[3], linewidth=basewidth)
l4 = lines!(ax1, 1:max_steps, Linear_approx_sampling_trace, color=color_lst[4], linewidth=basewidth)

ax2= Axis(  f[1,2],
            xticklabelsize=tickfontsize,
            yticklabelsize=tickfontsize,
            yticklabelpad=2,
            xlabel="Sample", ylabel = "Tr(Σ)",
            xlabelsize=labelfontsize,
            ylabelsize=labelfontsize,
            yscale = log10,
            yminorticksvisible = true,
            yminorgridvisible = true,
            yminorticks = IntervalsBetween(10),
            xminorticksvisible = true,
            xminorgridvisible = true,
            xminorticks = IntervalsBetween(10),
            title="Hopf Bifurcation",
            titlesize=labelfontsize)

ylims!(ax2,(5e-4,1))
xlims!(ax2,(0,1000))

l2 = lines!(ax2, 1:max_steps, Hopf_random_sampling_trace,      color=color_lst[1], linewidth=basewidth)
l3 = lines!(ax2, 1:max_steps, Hopf_random_sampling_trace ./ 2, color=color_lst[2], linewidth=basewidth+2)
l1 = lines!(ax2, 1:max_steps, Hopf_optimal_sampling_trace,     color=color_lst[3], linewidth=basewidth)
l4 = lines!(ax2, 1:max_steps, Hopf_approx_sampling_trace,      color=color_lst[4], linewidth=basewidth)

ax3= Axis(  f[1,3],
            xticklabelsize=tickfontsize,
            yticklabelsize=tickfontsize,
            yticklabelpad=2,
            xlabel="Sample", ylabel = "Tr(Σ)",
            xlabelsize=labelfontsize,
            ylabelsize=labelfontsize,
            yscale = log10,
            yminorticksvisible = true,
            yminorgridvisible = true,
            yminorticks = IntervalsBetween(10),
            xminorticksvisible = true,
            xminorgridvisible = true,
            xminorticks = IntervalsBetween(10),
            title="Lorenz System",
            titlesize=labelfontsize)

ylims!(ax3,(3e-3,0.1))
xlims!(ax3,(0,1000))

l2 = lines!(ax3, 1:max_steps, Lorenz_random_sampling_trace,           color=color_lst[1], linewidth=basewidth)
l3 = lines!(ax3, 1:max_steps, Lorenz_random_sampling_trace .* 2 ./ 3, color=color_lst[2], linewidth=basewidth+2)
l1 = lines!(ax3, 1:max_steps, Lorenz_optimal_sampling_trace,          color=color_lst[3], linewidth=basewidth)
l4 = lines!(ax3, 1:max_steps, Lorenz_approx_sampling_trace,           color=color_lst[4], linewidth=basewidth)

l1elem = LineElement(color=color_lst[3], linewidth=10)
l2elem = LineElement(color=color_lst[1], linewidth=10)
l3elem = LineElement(color=color_lst[2], linewidth=10)
l4elem = LineElement(color=color_lst[4], linewidth=10)

Legend(f[2,:],
       [l4elem, l1elem, l2elem, l3elem],
       ["1D Limit Set", "Dynamic Programming", "Random Sampling", "Dimension Reduction"],
       orientation=:horizontal,
       nbanks=2,
       labelsize=30,
       colgap=40)

colgap!(f.layout,40)

Label(f[1,1,TopLeft()], "A", font=noto_sans_bold, textsize = 32, halign=:left, valign=:bottom, padding=(25,0,10,0))
Label(f[1,2,TopLeft()], "B", font=noto_sans_bold, textsize = 32, halign=:left, valign=:bottom, padding=(25,0,10,0))
Label(f[1,3,TopLeft()], "C", font=noto_sans_bold, textsize = 32, halign=:left, valign=:bottom, padding=(50,0,10,0))
#Label(f[1,4,TopLeft()], "D", font=noto_sans_bold, textsize = 32, halign=:left, valign=:bottom, padding=(50,0,10,0))

save("out/FinalSim.pdf",f)


text:
include("check.jl")
include("move.jl")
include("hold.jl")

include("meta/index.jl")

include("angular/index.jl")
include("horizontal/index.jl")
include("vertical/index.jl")


text:
y = Dict{String,Array{String}}()
y["first"] = ["abc", "def"]
print(y["first"], "\n")
y["first"] = vcat(y["first"], ["ghi"])
print(y["first"], "\n")
push!(y["first"], "Ronn")
print(y["first"], "\n")

max_stars_repo_path: Julia Integration/testDict.jl

max_stars_repo_name: JayDweck/Ultimate-Datetime

max_issues_repo_path: Julia Integration/testDict.jl

max_issues_repo_name: JayDweck/Ultimate-Datetime

max_forks_repo_path: Julia Integration/testDict.jl

max_forks_repo_name: JayDweck/Ultimate-Datetime


text:
using TensorNetworkCodes
using TensorNetworkCodes: _find_product_indices
using LinearAlgebra: I
using Test

# BASIC FUNCTIONS

@testset "num_qubits" begin
    @test num_qubits(five_qubit_code()) == 5
    @test num_qubits(steane_code()) == 7
    @test num_qubits(random_code(4, 1)) == 4
    @test num_qubits(random_stabilizer_state(6)) == 6
end

@testset "verify" begin
    # test valid code
    valid_code = random_code(6, 2)
    @test verify(valid_code)

    # test invalid codes
    code_missing_pure_error = deepcopy(valid_code)
    pop!(code_missing_pure_error.pure_errors)
    @test !verify(code_missing_pure_error, log_warn=false)

    code_missing_logical = deepcopy(valid_code)
    pop!(code_missing_logical.logicals)
    @test !verify(code_missing_logical, log_warn=false)

    code_dependent_stabilizers = deepcopy(valid_code)
    pop!(code_dependent_stabilizers.stabilizers)
    push!(code_dependent_stabilizers.stabilizers, valid_code.stabilizers[1])
    @test !verify(code_dependent_stabilizers, log_warn=false)

    code_noncommuting_stabilizers = deepcopy(valid_code)
    code_noncommuting_stabilizers.stabilizers[2] = valid_code.pure_errors[1]
    @test !verify(code_noncommuting_stabilizers, log_warn=false)

    code_wrongorder_pure_errors = deepcopy(valid_code)
    pure_errors = code_wrongorder_pure_errors.pure_errors
    pure_errors[1], pure_errors[2] = pure_errors[2], pure_errors[1]
    @test !verify(code_wrongorder_pure_errors, log_warn=false)

    code_invalid_pure_errors = deepcopy(valid_code)
    pure_errors = code_invalid_pure_errors.pure_errors
    pure_errors[1] = pauli_product.(pure_errors[1], pure_errors[2])
    @test !verify(code_invalid_pure_errors, log_warn=false)

    code_noncommuting_logicals = deepcopy(valid_code)
    code_noncommuting_logicals.logicals[1] = valid_code.pure_errors[1]
    @test !verify(code_noncommuting_logicals, log_warn=false)

    # test logging
    log_pattern = (:warn, "number of stabilizers and pure errors don't match!")
    @test_logs log_pattern verify(code_missing_pure_error; log_warn=true)
    @test_logs verify(code_missing_pure_error; log_warn=false)
end

# EVALUATION FUNCTIONS

@testset "find_distance_logicals" begin
    d, ls = find_distance_logicals(five_qubit_code())
    @test d == 3 && all(==(3), (pauli_weight(l) for l in ls))

    d, ls = find_distance_logicals(steane_code())
    @test d == 3 && all(==(3), (pauli_weight(l) for l in ls))

    d, ls = find_distance_logicals(five_qubit_surface_code())
    @test d == 2 && all(==(2), (pauli_weight(l) for l in ls))

    @test_throws ErrorException find_distance_logicals(five_qubit_code(); max_distance=2)
end

@testset "_find_product_indices" begin
    operators = [[1, 3, 3, 1, 0], [0, 1, 3, 3, 1], [1, 0, 1, 3, 3], [3, 1, 0, 1, 3]]
    target = pauli_product_pow(operators, [1, 0, 1, 1])
    target_indices = _find_product_indices(operators, target)
    found = pauli_product(operators[target_indices])
    @test found == target
    # identity target
    @test _find_product_indices(operators, [0, 0, 0, 0, 0]) == []
end

@testset "find_pure_error/find_syndrome" begin
    code = steane_code()
    error = [0, 1, 3, 0, 0, 2, 0]
    syndrome = find_syndrome(code, error)
    pure_error = find_pure_error(code, syndrome)
    @test find_syndrome(code, pure_error) == syndrome

    #identity error
    identity = zeros(Int, 7)
    syndrome = find_syndrome(code, identity)
    @test all(==(0), syndrome)
    @test find_pure_error(code, syndrome) == identity
end

@testset "find_pure_errors" begin
    # 5-qubit code stabilizers: XZZXI, IXZZX, XIXZZ, ZXIXZ
    stabilizers = [[1, 3, 3, 1, 0], [0, 1, 3, 3, 1], [1, 0, 1, 3, 3], [3, 1, 0, 1, 3]]
    stabilizers_copy = deepcopy(stabilizers)
    pure_errors = find_pure_errors(stabilizers)
    @test stabilizers == stabilizers_copy  # stabilizers not modified
    # test commutation relations
    @test length(pure_errors) == length(stabilizers)
    @test [pauli_commutation(s, p) for s in stabilizers, p in pure_errors] == I
end

# TRANSFORMATION FUNCTIONS

@testset "gauge" begin
    # simple code
    code = five_qubit_code()
    code_copy = deepcopy(code)
    new_code = gauge(code, 1, 3)
    # test code not modified
    @test code.name == code_copy.name
    @test code.stabilizers == code_copy.stabilizers
    @test verify(code)
    # test new_code modified as expected
    @test new_code.name == "1/Z gauged $(code.name)"
    @test num_qubits(new_code) == num_qubits(code)
    @test new_code.stabilizers == vcat(code.stabilizers, [code.logicals[2]])
    @test length(new_code.logicals) == 0
    @test verify(new_code)
    # test exceptions
    @test_throws ErrorException gauge(five_qubit_code(), 2, 1)  # invalid logical qubit
    @test_throws ErrorException gauge(five_qubit_code(), 1, 4)  # invalid pauli
    @test_throws ErrorException gauge(five_qubit_code(), 1, 0)  # gauging with identity

    # tn code
    tn_code = TensorNetworkCode(five_qubit_code())
    tn_new_code = gauge(tn_code, 1, 3)
    @test verify(tn_new_code)
    @test tn_new_code.stabilizers == new_code.stabilizers
    @test tn_new_code.logicals == new_code.logicals
end

@testset "permute" begin
    code = five_qubit_code()
    code_copy = deepcopy(code)
    permutation = [5, 4, 3, 2, 1]
    new_code = permute(code, permutation)
    # test code not modified
    @test code.name == code_copy.name
    @test code.stabilizers == code_copy.stabilizers
    @test verify(code)
    # test new_code modified as expected
    @test new_code.name == "$(code.name) $(permutation)"
    @test new_code.stabilizers == getindex.(code.stabilizers, Ref(permutation))
    @test verify(new_code)
end

@testset "purify" begin
    code = five_qubit_code()
    code_copy = deepcopy(code)
    new_code = purify(code)
    # test code not modified
    @test code.name == code_copy.name
    @test code.stabilizers == code_copy.stabilizers
    @test verify(code)
    # test new_code modified as expected
    @test new_code.name == "Purified $(code.name)"
    @test num_qubits(new_code) == num_qubits(code) + length(code.logicals) / 2
    @test length(new_code.stabilizers) == num_qubits(new_code)
    @test length(new_code.logicals) == 0
    @test verify(new_code)
end


text:
export GnuplotScript
export register_data
export free_form
export plot, replot
export set_title
export add_vertical_line
export write_script

using DelimitedFiles

# Data identifier
#
struct RegisteredData_UUID
    _hash::UInt64
    RegisteredData_UUID(any) = new(hash(any))
end

# Show a well formed gnuplot id
import Base: show
Base.show(io::IO,uuid::RegisteredData_UUID) = print(io,"\$G"*string(uuid._hash))

# Gnuplot exe
# TODO: windows!
#
@static if Sys.iswindows()
    const gnuplot_exe = "gnuplot.exe"
else
    const gnuplot_exe = "gnuplot"
end

@doc raw"""
```julia
gp = GnuplotScript(;direct_plot = true)
```

Create a gnuplot script `gp`. If `direct_plot` is true, simultaneously
plot the registered operations.

# Usage example

You can perform a simple plot as follows:

```julia
gp = GnuplotScript(;direct_plot = true)

X=[-pi:0.1:pi;];
Ys = sin.(X);
Yc = cos.(X);

id = register_data(gp,hcat(X,Ys,Yc))
free_form(gp,"replot '$id' u 1:3 w l t 'cos'")
free_form(gp,"replot '$id' u 1:2 w l t 'sin'")
```

The plot will be created immediately.

# Also see
- [`register_data`](@ref)
- [`free_form`](@ref)
"""
mutable struct GnuplotScript
    _registered_data::Dict{RegisteredData_UUID,Any}
    _script::String
    _any_plot::Bool
    _direct_plot_io::Union{Nothing,IO}

    function GnuplotScript(;direct_plot = true)
        # manage direct plot
        #
        io = nothing
        if direct_plot
            io = open(Cmd([gnuplot_exe]),write=true)
            if !iswritable(io)
                @warn "Cannot start gnuplot executable $gnuplot_exe"
                close(io)
                io = nothing
            end
        end

        # other default values
        #
        registered_data=Dict{RegisteredData_UUID,Any}()
        script = String("")
        any_plot = false

        gp = new(registered_data,script,any_plot,io)
    end
end

# check if data has already been registered
#
function is_registered(gp::GnuplotScript,uuid::RegisteredData_UUID)
    haskey(gp._registered_data,uuid)
end

# write data
# ID << EOD
# data
# END
function _write_data(io::IO, uuid::RegisteredData_UUID, data::AbstractVecOrMat)
    println(io,"$uuid << EOD")
    writedlm(io,data)
    println(io,"EOD")
end

function _write_data(io::IO,gp::GnuplotScript)
    for (k,d) in gp._registered_data
        _write_data(io,k,d)
    end
end

function _write_script(io::IO,gp::GnuplotScript)
    _write_data(io,gp)
    print(io,gp._script)
end

# ================================================================

# private methods to write data:
# -> use them systematically. Reason contains the extra logic used for direct plotting
#
function _append_data(gp::GnuplotScript,data::AbstractVecOrMat)::RegisteredData_UUID
    # already registered
    uuid = RegisteredData_UUID(data)

    # no, register it
    if !is_registered(gp,uuid)
        gp._registered_data[uuid]=data
    end

    if gp._direct_plot_io != nothing
        try
            _write_data(gp._direct_plot_io,uuid,data)
        catch
            @warn "Broken pipe"
            close(gp._direct_plot_io)
            gp._direct_plot_io=nothing
        end
    end

    uuid
end

function _append_to_script(gp::GnuplotScript,line::AbstractString)
    # Append to script
    #
    line *= "\n"
    gp._script *= line

    # If direct plot is active
    # also pipe the line
    #
    if gp._direct_plot_io != nothing
        try
            write(gp._direct_plot_io,line)
        catch
            @warn "Broken pipe"
            close(gp._direct_plot_io)
            gp._direct_plot_io=nothing
        end
    end

    gp
end

@doc raw"""
```julia
register_data(gp::GnuplotScript,
              data::AbstractVecOrMat;
              copy_data::Bool=true) -> id
```

Register data and return the associated data identifier. Registered
data is embedded in the plot script file. The returned `id` is used to
reference registered data.

# Usage example

```julia
gp = GnuplotScript()

M = rand(10,3)

id = register_data(gp, M)

free_form(gp,"replot $id u 1:2")
free_form(gp,"replot $id u 1:3")
```

"""
function register_data(gp::GnuplotScript,data::AbstractVecOrMat;
                       copy_data::Bool=true)::RegisteredData_UUID
    if copy_data
        data = deepcopy(data)
    end

    _append_data(gp,data)
end

# Detect first plot when using free_form
#
function _contains_plot_p(gp_line::AbstractString)
    space = "^([ \t]*)"
    r = "$(space)splot |$(space)plot "
    r = Regex(r)

    match(r,gp_line) !== nothing
end

@doc raw"""
```julia
free_form(gp::GnuplotScript,gp_line::AbstractString)
```

Write gnuplot commands. This command line is directly forwarded to
Gnuplot. The only difference is that you can use `replot` even for the
first plot. This is convenient when you chain plots, you do not have
to worry if the current command is the first plot.

# Usage example

```julia
using GnuplotScripting

gp = GnuplotScript()

free_form(gp, "replot sin(x) lw 2 t 'a trigonometric function'")
```

"""
function free_form(gp::GnuplotScript,gp_line::AbstractString)
    # Detect replot instruction.
    # If no previous plot, do "replot" => "plot"
    #
    r = r"^([ \t]*)(replot) "
    if match(r,gp_line) !== nothing
        if gp._any_plot == false
            gp_line = replace(gp_line,r"^([ \t]*)(replot) "=>s"\g<1>plot ")
        end
    end

    # Detect plot instruction
    # If one is detected, set any_plot = true
    #
    if _contains_plot_p(gp_line)
        gp._any_plot = true
    end
    #
    _append_to_script(gp,gp_line)
end

@doc raw"""
```julia
set_title(gp::GnuplotScript,title::AbstractString;
                   enhanced::Bool = false)
```

Define plot title. If `enhanced` is true, some characters are
processed in a special way. By example `_` subscripts text.
"""
function set_title(gp::GnuplotScript,title::AbstractString;
                   enhanced::Bool = false)
    command = "set title '$title'"
    if enhanced==false
        command *= " noenhanced"
    end

    _append_to_script(gp,command)
end

function _plot(gp::GnuplotScript,uuid::RegisteredData_UUID,plot_arg::AbstractString;
               # true: plot triggered by gnuplot "plot"
               reset_plot::Bool)
    @assert is_registered(gp,uuid)

    # Detect replot
    #
    replot = false
    if reset_plot
        gp._any_plot = false
    else
        if gp._any_plot
            replot = true
        end
    end

    # create command
    #
    command = ""
    if replot
        command *= "re"
    end
    #
    command *= "plot $uuid " * plot_arg

    _append_to_script(gp,command)

    gp._any_plot = true
end

# ****************************************************************
# Dedicated function / recipes
# ****************************************************************
#

function plot(gp::GnuplotScript,uuid::RegisteredData_UUID,plot_arg::AbstractString)
    _plot(gp,uuid,plot_arg,reset_plot=true)
end

# Like gnuplot "replot", with the difference that we automatically
# switch to "plot" when there is no initial plot.
#
# Note that we can force "replot" thanks to "force_replot"
#
function replot(gp::GnuplotScript,uuid::RegisteredData_UUID,plot_arg::AbstractString;
                force_replot::Bool=false)

    if force_replot
        gp._any_plot = true
    end

    _plot(gp,uuid,plot_arg,reset_plot=false)
end

# vertical line ================
#

"""
    add_vertical_line(gp::GnuplotScript,position::Float64;name::Union{AbstractString,Nothing})

Add a vertical bar with a label
"""
function add_vertical_line(gp::GnuplotScript,position::Float64;name::Union{AbstractString,Nothing})
    command = ""
    if name != Nothing
        command  *= "set label at $position, 0.0 '$name' rotate by 90 front left offset -1,1,0 tc ls 1\n"
    end
    command *= "set arrow from $position, graph 0 to $position, graph 1 nohead front\n"

    _append_to_script(gp,command)
end

# ****************************************************************
# Write script
# ****************************************************************
#

"""
    write_script(script_file::AbstractString,gp::GnuplotScript)

Write script with embedded data for future use.

# Usage

```julia
gp = GnuplotScript()

...

write_script("gnuplot_script.gp",gp)
```

You can replay the script using Gnuplot:
```sh
    gnuplot gnuplot_script.gp
```

If you want to keep the gnuplot session opened, add a final `-`

```sh
    gnuplot gnuplot_script.gp -
```
"""
function write_script(script_file::AbstractString,gp::GnuplotScript)
    io = open(script_file, "w");

    _write_script(io,gp)

    close(io)
end


text:
# # Laub-Loomis
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/LaubLoomis.ipynb)
#
#md # !!! note "Overview"
#md #     System type: polynomial with quadratic nonlinearities\
#md #     State dimension: 7\
#md #     Application domain: Molecular Biology
#
# ## Model description
#
# The Laub-Loomis model is presented in [^LL98] for studying a class of enzymatic activities.
# The dynamics can be defined by the following ODE with 7 variables.

using ReachabilityAnalysis, Plots

@taylorize function laubloomis!(dx, x, params, t)
    dx[1] = 1.4*x[3] - 0.9*x[1]
    dx[2] = 2.5*x[5] - 1.5*x[2]
    dx[3] = 0.6*x[7] - 0.8*(x[2]*x[3])
    dx[4] = 2 - 1.3*(x[3]*x[4])
    dx[5] = 0.7*x[1] - (x[4]*x[5])
    dx[6] = 0.3*x[1] - 3.1*x[6]
    dx[7] = 1.8*x[6] - 1.6*(x[2]*x[7])
    return dx
end

# The system is asymptotically stable and the equilibrium is the origin.

# ## Specifications
#
# The specification for the analysis is given next.
# The initial sets are defined according to the ones used in [^TD13].
# They are boxes centered at
# ``x_1(0) = 1.2``, ``x_2(0) = 1.05``, ``x_3(0) = 1.5``, ``x_4(0) = 2.4``,
# ``x_5(0) = 1``, ``x_6(0) = 0.1``, ``x_7 (0) = 0.45``.
# The range of the box in the i-th dimension is defined by the interval ``[x_i(0) − W, x_i(0) + W]``.
# The width ``W`` of the initial set is vital to the difficulty of the reachability analysis job.
# The larger the initial set the harder the reachability analysis.

function laubloomis(; W=0.01)
    ## initial states
    X0c = [1.2, 1.05, 1.5, 2.4, 1.0, 0.1, 0.45]
    X0 = Hyperrectangle(X0c, fill(W, 7))

    ## initial-value problem
    prob = @ivp(x' = laubloomis!(x), dim: 7, x(0) ∈ X0)

    return prob
end

# We consider ``W = 0.01``, ``W = 0.05``, and ``W = 0.1``. The specification for each
# scenario is the following:
#
# - ``W = 0.01`` and ``W = 0.05``: the unsafe region is defined by ``x_4 ≥ 4.5
# - ``W = 0.1``: the unsafe set is defined by ``x_4 ≥ 5``.
#
# The time horizon for all cases is [0, 20].

# ## Results

# The final widths of ``x_4`` along with the computation times are provided for all three cases. A
# figure is provided in the ``(t, x_4)`` axes, with ``t ∈ [0, 20]``, ``x_4 ∈ [1.5, 5]``,
# where the three plots are overlaid.

# Since the safety condition is only related to the variable ``x_4``, we present the
# plots of projections of the overapproximations in the ``t - x_4`` plane such that
# ``t`` is the time variable.

# ### Case 1: Smaller initial states

# In this case we consider thet the width of the initial states is ``W = 0.01``.

prob = laubloomis(W=0.01)
alg = TMJets(abstol=1e-11, orderT=7, orderQ=1, adaptive=true);

sol_1 = solve(prob, T=20.0, alg=alg);

sol_1z = overapproximate(sol_1, Zonotope);

# Verifying that the specifications hold:

## canonical direction along x₄
const e4 = [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0];

ρ(e4, sol_1z)
#-
ρ(e4, sol_1z) < 4.5

# To compute the width of the final box width of the final box we use support functions:
@show ρ(e4, sol_1z[end]) + ρ(-e4, sol_1z[end])

# ### Case 2: Intermediate initial states

# In this case we consider the the width of the initial states is ``W = 0.05``.

prob = laubloomis(W=0.05)
alg = TMJets(abstol=1e-12, orderT=7, orderQ=1, adaptive=false);

sol_2 = solve(prob, T=20.0, alg=alg)
sol_2z = overapproximate(sol_2, Zonotope);

# Verifying that the specification holds:

ρ(e4, sol_2z)
#-
ρ(e4, sol_2z) < 4.5

## width of final box
@show ρ(e4, sol_2z[end]) + ρ(-e4, sol_2z[end])

# ### Case 3: Larger initial states

# In this case we consider that the width of the initial states is ``W = 0.1``.

prob = laubloomis(W=0.1)
alg = TMJets(abstol=1e-12, orderT=7, orderQ=1, adaptive=false);

sol_3 = solve(prob, T=20.0, alg=alg)
sol_3z = overapproximate(sol_3, Zonotope);

# Verifying that the specification holds:
ρ(e4, sol_3z)
#-
ρ(e4, sol_3z) < 5.0

## width of final box
@show ρ(e4, sol_3z[end]) + ρ(-e4, sol_3z[end])

using Plots, Plots.PlotMeasures, LaTeXStrings #!jl

fig = plot()
plot!(fig, sol_3z, vars=(0, 4), linecolor="green", color=:green, alpha=0.8)
plot!(fig, sol_2z, vars=(0, 4), linecolor="blue",  color=:blue, alpha=0.8)
plot!(fig, sol_1z, vars=(0, 4), linecolor="yellow", color=:yellow, alpha=0.8,
      tickfont=font(10, "Times"), guidefontsize=15,
      xlab=L"t",
      ylab=L"x_4",
      xtick=[0., 5., 10., 15., 20.], ytick=[1.5, 2., 2.5, 3., 3.5, 4., 4.5, 5.],
      xlims=(0., 20.), ylims=(1.5, 5.02),
      bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm,
      size=(600, 600))

plot!(fig, x->x, x->4.5, 0., 20., line=2, color="red", linestyle=:dash, legend=nothing)
plot!(fig, x->x, x->5., 0., 20., line=2, color="red", linestyle=:dash, legend=nothing)

# ## References

# [^LL98]: M. T. Laub and W. F. Loomis. *A molecular network that produces spontaneous oscillations in excitable cells of dictyostelium.* Molecular Biology of the Cell, 9:3521–3532, 1998.

# [^TD13]: R. Testylier and T. Dang. *NLTOOLBOX: A library for reachability computation of nonlinear dynamical systems.* In Proc. of ATVA’13, volume 8172 of LNCS, pages 469–473. Springer, 2013.


text:
import LinearAlgebra.dot
function dot(u :: AbstractLinearOperator, v :: Union{AbstractVector, AbstractLinearOperator})

 if size(u,1) == 1
     vu = Matrix(u')[:]
 elseif size(u,2) == 1
     vu = Matrix(u)[:]
 else
     throw("Wrong dimension")
 end
 if size(v,1) == 1
     vv = Matrix(v')[:]
 elseif size(v,2) == 1
     vv = typeof(v) <: AbstractVector ? v : Matrix(v)[:]
 else
     throw("Wrong dimension")
 end
 if length(vu) != length(vv)
     throw("Wrong dimension")
 end
 return dot(vu, vv)
end

import Base.-
function -(u :: AbstractVector, v :: AbstractLinearOperator)
 vv = zeros(length(u))
 if size(v,1) == 1
     vv = Matrix(v')[:]
 elseif size(v,2) == 1
     vv = Matrix(v)[:]
 else
     throw("Wrong dimension")
 end
 if length(u) != length(vv)
     throw("Wrong dimension")
 end

 return u - vv
end


text:
using ITensors
using LinearAlgebra
using Zygote
using Random

Random.seed!(1234)

function heisenberg_hamiltonian(N)
  os = OpSum()
  for j in 1:(N - 1)
    os += 0.5, "S+", j, "S-", j + 1
    os += 0.5, "S-", j, "S+", j + 1
    os += "Sz", j, "Sz", j + 1
  end
  return os
end

# Number of sites in the system
N = 6
s = siteinds("S=1/2", N)
os = heisenberg_hamiltonian(N)
Hmpo = MPO(os, s)

H = prod(Hmpo)

# Number of eigenstates to target
p = 3

D, V = eigen(H; ishermitian=true)
exact_energies = sort(ITensors.data(D))[1:p]

n = isqrt(dim(H))
sr = Index(p, "renorm")
X₀ = randomITensor(s..., sr)
U₀, _ = polar(X₀, s)
U₀ = noprime(U₀)

# Loss function for low energy eigenspace.
# Reference: https://arxiv.org/abs/2007.01287, page 9, example A
function E(H, U)
  UHU = dag(U) * H * U'
  trUHU = (UHU * dag(δ(inds(UHU))))[]
  return trUHU
end

loss(U) = E(H, U)

@show sum(exact_energies)
@show loss(U₀)

# Step size
η = 0.1

# Move to another point of the Stiefel manifold using the
# SVD retraction formula.
# Reference: https://arxiv.org/abs/2007.01287, page 8, equation (17)
RSVD(X, W, is) = noprime(polar(X + W, is)[1])

# Riemannian optimization of Stiefel Manifold (manifold of isometric matrices)
# Reference: https://arxiv.org/abs/2007.01287, page 8, algorithm steps 3-7
function update(U₀)
  ∇ = loss'(U₀)
  U₀∇ = dag(U₀) * prime(∇; inds=sr)
  ∇U₀ = swapprime(dag(U₀∇), 0 => 1)
  ∇ᴿ = ∇ - 1 / 2 * noprime(U₀ * (U₀∇ + ∇U₀))
  U = RSVD(U₀, -η * ∇ᴿ, s)
  return U
end

function update(U, N)
  for n in 1:N
    U = update(U)
    @show n, loss(U)
  end
  return U
end

U = U₀
U = update(U, 50)

@show loss(U)

# Compare against ED
D, _ = eigen(H; ishermitian=true)
energy_exact = sum(n -> D[n, n], 1:p)
@show loss(U), energy_exact


text:
@testset "SSOR" begin
	@testset "small full system" begin
		# CG: test sparse Laplacian
		A = getDivGrad(32,32,32)
		rhs = randn(size(A,1))
		tol = 1e-2

		# tests with A being matrix
		x1,flag1,relres1,iter1,resvec1       = ssor(A,rhs,tol=tol,maxIter=100)
		x2,flag2,relres2,iter2,resvec2       = ssor(A,rhs,tol=tol,maxIter=100,omega=2/3,out=1)
		x3,flag3,relres3,iter3,resvec3     = ssor(A,rhs,tol=tol,maxIter=100,omega=1/2,out=1,storeInterm=true)

		# test relative residuals
		@test norm(A*x1-rhs)/norm(rhs) <= tol
		@test norm(A*x2-rhs)/norm(rhs) <= tol
		@test norm(A*x3[:,end]-rhs)/norm(rhs) <= tol
	end

	@testset "preconditioner" begin
		omega = 1.2;
		A = getDivGrad(32,32,32)
		rhs = randn(size(A,1))
		d = omega./Array(diag(A));
		x = zeros(length(rhs)) # pre allocation for the preconditioner result.
		PC = r -> (x[:].=0.0; return ssorPrecTrans!(A,x,r,d));
		y = KrylovMethods.cg(A,rhs,tol=1e-12,maxIter=200,M=PC,out=1)[1]
		@test norm(A*y-rhs)/norm(rhs) <= 1e-12
	end

	@testset "complex " begin

		# CG: test sparse Laplacian
		A = getDivGrad(32,32,32)
		A += 1im*I
		rhs = randn(size(A,1)) + 1im * randn(size(A,1))
		tol = 1e-2

		#
		x0,flag0,relres0,iter0,resvec0       = ssor(A,rhs,tol=tol,maxIter=3,out=2)
		@test flag0==-1
		@test iter0==3
		@test length(resvec0)==3

		# tests with A being matrix
		x1,flag1,relres1,iter1,resvec1       = ssor(A,rhs,tol=tol,maxIter=100)
		x2,flag2,relres2,iter2,resvec2       = ssor(A,rhs,tol=tol,maxIter=100,omega=2/3,out=1)
		x3,flag3,relres3,iter3,resvec3       = ssor(A,rhs,tol=tol,maxIter=100,omega=1/2,out=1)

		# test relative residuals
		@test norm(A*x1-rhs)/norm(rhs) <= tol
		@test norm(A*x2-rhs)/norm(rhs) <= tol
		@test all(diff(resvec1).<0)
		@test all(diff(resvec2).<0)
		@test all(diff(resvec3).<0)

	end
end


text:
"""
```
filter(m, data, system, s_0 = [], P_0 = []; cond_type = :none,
       include_presample = true, in_sample = true,
       outputs = [:loglh, :pred, :filt])
```

Computes and returns the filtered values of states for the state-space
system corresponding to the current parameter values of model `m`.

### Inputs

- `m::AbstractDSGEModel`: model object
- `data`: `DataFrame` or `nobs` x `hist_periods` `Matrix{S}` of data for
  observables. This should include the conditional period if `cond_type in
  [:semi, :full]`
- `system::System`: `System` object specifying state-space system matrices for
  the model
- `s_0::Vector{S}`: optional `Nz` x 1 initial state vector
- `P_0::Matrix{S}`: optional `Nz` x `Nz` initial state covariance matrix

where `S<:AbstractFloat`.

### Keyword Arguments

- `cond_type::Symbol`: conditional case. See `forecast_all` for documentation of
  all `cond_type` options
- `include_presample::Bool`: indicates whether to include presample periods in
  the returned vector of `Kalman` objects
- `in_sample::Bool`: indicates whether or not to discard out of sample rows in
    `df_to_matrix` call
- `outputs::Vector{Symbol}`: which Kalman filter outputs to compute and return.
  See `?kalman_filter`

### Outputs

- `kal::Kalman`: see `?Kalman`
"""
function filter(m::AbstractDSGEModel, df::DataFrame, system::System{S},
                s_0::Vector{S} = Vector{S}(undef, 0),
                P_0::Matrix{S} = Matrix{S}(undef, 0, 0);
                cond_type::Symbol = :none, include_presample::Bool = true,
                in_sample::Bool = true,
                outputs::Vector{Symbol} = [:loglh, :pred, :filt],
                tol::Float64 = 0.0) where {S<:AbstractFloat}

    data = df_to_matrix(m, df; cond_type = cond_type, in_sample = in_sample)
    start_date = max(date_presample_start(m), df[1, :date])
    filter(m, data, system, s_0, P_0; start_date = start_date,
           include_presample = include_presample, outputs = outputs, tol = tol)
end

function filter(m::AbstractDSGEModel, data::AbstractArray, system::System{S},
                s_0::Vector{S} = Vector{S}(undef, 0),
                P_0::Matrix{S} = Matrix{S}(undef, 0, 0);
                start_date::Date = date_presample_start(m),
                include_presample::Bool = true,
                outputs::Vector{Symbol} = [:loglh, :pred, :filt],
                tol::Float64 = 0.0) where {S<:AbstractFloat}

    T = size(data, 2)

    # Partition sample into pre- and post-ZLB regimes
    # Note that the post-ZLB regime may be empty if we do not impose the ZLB
    regime_inds = zlb_regime_indices(m, data, start_date)

    # Get system matrices for each regime
    TTTs, RRRs, CCCs, QQs, ZZs, DDs, EEs = zlb_regime_matrices(m, system, start_date)

    # If the end of the first regime is assumed to be later than the data provided
    # change the regime to be up to the length of the data
    if regime_inds[1][end] > T
        regime_inds = [1:T]
        TTTs, RRRs, CCCs, QQs, ZZs, DDs, EEs = [TTTs[1]], [RRRs[1]], [CCCs[1]], [QQs[1]], [ZZs[1]], [DDs[1]], [EEs[1]]
    end

    # If s_0 and P_0 provided, check that rows and columns corresponding to
    # anticipated shocks are zero in P_0
    if !isempty(s_0) && !isempty(P_0)
        ant_state_inds = setdiff(1:n_states_augmented(m), inds_states_no_ant(m))
        @assert all(x -> x == 0, P_0[:, ant_state_inds])
        @assert all(x -> x == 0, P_0[ant_state_inds, :])
    end

    # Specify number of presample periods if we don't want to include them in
    # the final results
    Nt0 = include_presample ? 0 : n_presample_periods(m)

    # Run Kalman filter, construct Kalman object, and return
    out = kalman_filter(regime_inds, data, TTTs, RRRs, CCCs, QQs,
                        ZZs, DDs, EEs, s_0, P_0; outputs = outputs,
                        Nt0 = Nt0, tol = tol)
    return Kalman(out...)
end

function filter_likelihood(m::AbstractDSGEModel, df::DataFrame, system::System{S},
                           s_0::Vector{S} = Vector{S}(undef, 0),
                           P_0::Matrix{S} = Matrix{S}(undef, 0, 0);
                           cond_type::Symbol = :none, include_presample::Bool = true,
                           in_sample::Bool = true,
                           tol::Float64 = 0.0) where {S<:AbstractFloat}

    data = df_to_matrix(m, df; cond_type = cond_type, in_sample = in_sample)
    start_date = max(date_presample_start(m), df[1, :date])

    filter_likelihood(m, data, system, s_0, P_0; start_date = start_date,
                      include_presample = include_presample, tol = tol)
end

function filter_likelihood(m::AbstractDSGEModel, data::AbstractArray, system::System{S},
                           s_0::Vector{S} = Vector{S}(undef, 0),
                           P_0::Matrix{S} = Matrix{S}(undef, 0, 0);
                           start_date::Date = date_presample_start(m),
                           include_presample::Bool = true,
                           tol::Float64 = 0.0) where {S<:AbstractFloat}

    # Partition sample into pre- and post-ZLB regimes
    # Note that the post-ZLB regime may be empty if we do not impose the ZLB
    regime_inds = zlb_regime_indices(m, data, start_date)

    # Get system matrices for each regime
    TTTs, RRRs, CCCs, QQs, ZZs, DDs, EEs = zlb_regime_matrices(m, system, start_date)

    # If s_0 and P_0 provided, check that rows and columns corresponding to
    # anticipated shocks are zero in P_0
    if !isempty(s_0) && !isempty(P_0)
        ant_state_inds = setdiff(1:n_states_augmented(m), inds_states_no_ant(m))
        @assert all(x -> x == 0, P_0[:, ant_state_inds])
        @assert all(x -> x == 0, P_0[ant_state_inds, :])
    end

    # Specify number of presample periods if we don't want to include them in
    # the final results
    Nt0 = include_presample ? 0 : n_presample_periods(m)

    # Run Kalman filter, construct Kalman object, and return
    kalman_likelihood(regime_inds, data, TTTs, RRRs, CCCs, QQs,
                      ZZs, DDs, EEs, s_0, P_0; Nt0 = Nt0, tol = tol)
end

function filter_shocks(m::AbstractDSGEModel, df::DataFrame, system::System{S},
                       s_0::Vector{S} = Vector{S}(0), P_0::Matrix{S} = Matrix{S}(0, 0); cond_type::Symbol = :none,
                       start_date::Date = date_presample_start(m),
                       include_presample::Bool = false) where S<:AbstractFloat

    data = df_to_matrix(m, df, cond_type = cond_type)

    # Partition sample into pre- and post-ZLB regimes
    # Note that the post-ZLB regime may be empty if we do not impose the ZLB
    regime_inds = zlb_regime_indices(m, data, start_date)

    # Get system matrices for each regime
    Ts, Rs, Cs, Qs, Zs, Ds, Es = zlb_regime_matrices(m, system, start_date)

    # Initialize s_0 and P_0
    if isempty(s_0) || isempty(P_0)
        s_0, P_0 = init_stationary_states(Ts[1], Rs[1], Cs[1], Qs[1])

        # If s_0 and P_0 provided, check that rows and columns corresponding to
        # anticipated shocks are zero in P_0
        ant_state_inds = setdiff(1:n_states_augmented(m), inds_states_no_ant(m))
        # MDχ: Had to change the tolerance from exactly == 0, if below tol
        # explicitly set to 0
        @assert all(x -> abs(x) < 1e-14, P_0[:, ant_state_inds])
        @assert all(x -> abs(x) < 1e-14, P_0[ant_state_inds, :])
        P_0[:, ant_state_inds] = 0.
        P_0[ant_state_inds, :] = 0.
    end

    # Specify number of presample periods if we don't want to include them in
    # the final results
    Nt0 = include_presample ? 0 : n_presample_periods(m)

    # Dimensions
    Nt = size(data,  2) # number of periods of data
    Ns = size(Ts[1], 1) # number of states

    # Augment state space with shocks
    Ts, Rs, Cs, Zs, s_0, P_0 =
        augment_states_with_shocks(regime_inds, Ts, Rs, Cs, Qs, Zs, s_0, P_0)

    # Kalman filter stacked states and shocks stil_t
    _, _, _, stil_filt, Ptil_filt, _, _, _, _ =
        kalman_filter(regime_inds, data, Ts, Rs, Cs, Qs, Zs, Ds, Es, s_0, P_0,
                      outputs = [:filt])

    # Index out shocks
    ϵ_filt = stil_filt[Ns+1:end, :] # ϵ_{t|T}

    # Trim the presample if needed
    if Nt0 > 0
        insample = Nt0+1:Nt
        ϵ_filt = ϵ_filt[:, insample]
    end

    return ϵ_filt
end

"""
```
This section defines filter and filter_likelihood for the PoolModel type
```
"""
function filter(m::PoolModel, data::AbstractArray,
                s_0::AbstractArray{S} = Matrix{Float64}(undef,0,0);
                start_date::Date = date_presample_start(m),
                cond_type::Symbol = :none, include_presample::Bool = true,
                in_sample::Bool = true,
                tol::Float64 = 0., parallel::Bool = false,
                tuning::Dict{Symbol,Any} = Dict{Symbol,Any}()) where {S<:AbstractFloat}

    # Handle data
    Nt0 = include_presample ? 0 : n_presample_periods(m)
    if haskey(tuning, :n_presample_periods)
        tuning[:n_presample_periods] = Nt0
    end

    # Compute transition and measurement equations
    Φ, Ψ, F_ϵ, F_u, F_λ = compute_system(m)

    # Check initial states
    n_particles = haskey(tuning, :n_particles) ? tuning[:n_particles] : 1000
    if isempty(s_0)
        s_0 = reshape(rand(F_λ, n_particles), 1, n_particles)
        s_0 = [s_0; 1 .- s_0]
    elseif get_setting(m, :weight_type) == :dynamic
        if size(s_0,2) != n_particles
            error("s0 does not contain enough particles")
        end
    end

    # Check tuning
    if isempty(tuning)
        try
            tuning = get_setting(m, :tuning)
        catch
            if get_setting(m, :weight_type) == :dynamic
                @warn "no tuning parameters provided; using default tempered particle filter values"
            end
        end
    end
    if haskey(tuning, :parallel) # contains parallel? change keyword to that if so
        parallel = tuning[:parallel]
    end

    # Check if PoolModel has fixed_sched. If not, assume no tempering
    fixed_sched = [1.]
    try
        fixed_sched = get_setting(m, :fixed_sched)
    catch KeyError
    end

    weight_type = get_setting(m, :weight_type)
    if weight_type == :dynamic
        return tempered_particle_filter(data, Φ, Ψ, F_ϵ, F_u, s_0; parallel = parallel,
                                        poolmodel = true,
                                        fixed_sched = fixed_sched,
                                        tuning..., verbose = :none)
    elseif weight_type == :equal
        loglhconditional = log.(mapslices(x -> Ψ([0.], x), data, dims = 1))
        return sum(loglhconditional), loglhconditional
    elseif weight_type == :static
        loglhconditional = log.(mapslices(x -> Ψ([m[:λ].value; 1 - m[:λ].value], x), data, dims = 1))
        return sum(loglhconditional), loglhconditional
    elseif weight_type == :bma
        error("Estimation for Bayesian Model Averaging is computed directly by the estimate_bma function, so the filter function does not return anything.")
    end
end

function filter_likelihood(m::PoolModel, data::AbstractArray,
                           s_0::AbstractArray{S} = Matrix{Float64}(undef,0,0);
                           start_date::Date = date_presample_start(m),
                           cond_type::Symbol = :none, include_presample::Bool = true,
                           in_sample::Bool = true, parallel::Bool = false, tol::Float64 = 0.,
                           tuning::Dict{Symbol,Any} = Dict{Symbol,Any}()) where {S<:AbstractFloat}

    # Guarantee output settings in tuning give desired output
    tuning[:allout] = true
    tuning[:get_t_particle_dist] = false

    if get_setting(m, :weight_type) == :dynamic_weight
        ~, loglhconditional, ~ = filter(m, data, s_0; start_date = start_date,
                                        include_presample = include_presample,
                                        cond_type = cond_type, in_sample = in_sample, tol = tol,
                                        parallel = parallel, tuning = tuning)
    else
        ~, loglhconditional = filter(m, data, s_0; start_date = start_date,
                                     include_presample = include_presample,
                                     cond_type = cond_type, in_sample = in_sample, tol = tol,
                                     parallel = parallel, tuning = tuning)
    end

    return loglhconditional
end


text:
module TermStructure

using Ito.Time
using Calendar

export YieldTermStructure, VolatilityTermStructure, FlatYieldTermStructure,
		compound_factor, discount_factor, implied_rate, discount, forward_rate,
		reference_date, zero_rate, par_rate

#General Interest Rate functionality

compound_factor(r::Real, compounding::Symbol, freq::Symbol,  t::Real) = compound_factor(r, compounding, eval(freq), t)
compound_factor(r::Real, compounding::Symbol,  t::Real) = compound_factor(r, compounding, NoFrequency, t)
compound_factor(r::Real, compounding::Symbol, dc::DayCount, dates::CalendarTime...) = compound_factor(r, compounding, NoFrequency, dates...)
compound_factor(r::Real, compounding::Symbol, freq::Symbol, dc::DayCount, dates::CalendarTime...) = compound_factor(r, compounding, eval(freq), dates...)
compound_factor(r::Real, compounding::Symbol, freq::Integer, dc::DayCount, dates::CalendarTime...) = compound_factor(r,  compounding, freq, yearfraction(dc, dates...))

function compound_factor(r::Real, compounding::Symbol, freq::Integer, t::Real )
	if compounding == :Simple
		return 1 + r * t
	elseif compounding == :Compounded
		@assert(freq != NoFrequency)
		return (1 + r / freq) ^ (freq * t)
	elseif compounding == :Continuous
		return exp(r * t)
	elseif compounding == :SimpleThenCompounded
		@assert(freq != NoFrequency)
		if t < (1 / freq)
            return 1 + r * t
        else
            return (1 + r / freq) ^ (freq * t)
        end
	else
		error("Unknown compounding")
	end
end

discount_factor(x...) = 1/compound_factor(x...)

implied_rate(c::Real, compounding::Symbol, freq::Symbol,  t::Real) = implied_rate(c,compounding, eval(freq), t)
implied_rate(c::Real, compounding::Symbol, freq::Symbol, dc::DayCount, dates::CalendarTime...) = implied_rate(c,compounding, eval(freq), dates...)
implied_rate(c::Real, compounding::Symbol, freq::Integer, dc::DayCount, dates::CalendarTime...) = implied_rate(c, compounding, freq, yearfraction(dc, dates...))

function implied_rate(c::Real, compounding::Symbol, freq::Integer, t::Real)
	@assert c>0 && t>0
	if compounding == :Simple
		return (c-1) / t
	elseif compounding == :Compounded
		return (c ^ (1/(f*t)) - 1) *f
	elseif compounding == :Continuous
		return log(c)/t
	elseif compounding == :SimpleThenCompounded
		if t < (1 / freq)
            return (c-1) / t
        else
            return (c ^ (1/(f*t)) - 1) *f
        end
	else
		error("Unknown compounding")
	end
end

#abstract TermStructure

abstract VolatilityTermStructure
abstract YieldTermStructure

discount(ts::YieldTermStructure, d::CalendarTime) = discount(ts, yearfraction(reference_date(ts), d))

#Overload this method for each concrete YieldTermStructure
discount(ts::YieldTermStructure, t::Real) = error("Must be implemented by concrete Term Structure")

#Conventional implementation. Usually re-implemented by concrete term structures
reference_date(ts::YieldTermStructure) = ts.reference_date

forward_rate(ts::YieldTermStructure, compounding::Symbol, freq::Integer, d1::CalendarTime, d2::CalendarTime) = forward_rate(ts, compounding, freq, d1, d2)
function forward_rate(ts::YieldTermStructure, compounding::Symbol, freq::Integer, d1::CalendarTime, d2::CalendarTime )
	if d1==d2
		t1 = yearfraction(reference_date(ts), d1)
		t2 = t1+.0001
		c=discount(ts, t1) / discount(ts, t2)
		return implied_rate(c, compound, freq, delta)
	elseif d1<d2
		return implied_rate(discount(ts, d1)/discount(ts, d2), ts.dc, d1, d2)
	else
		error("Forward start date must be before forward end dates")
	end
end

function forward_rate(ts::YieldTermStructure, compounding::Symbol, freq::Integer, t1::Real, t2::Real )
	if (t2==t1)
		t2=t1+.0001
	end

	compound = discount(ts, t1) / discount(ts, t2)
	return implied_rate(discount(ts, t1) / discount(ts, t2), ts.dc, t1, t2)
end

zero_rate(ts::YieldTermStructure, compounding::Symbol, freq::Integer, d1::CalendarTime) = zero_rate(ts, compounding, freq, yearfraction(ts.dc, reference_date(ts), d1))
function zero_rate(ts::YieldTermStructure, compounding::Symbol, freq::Integer, t::Real)
	if (t == 0)
		c = 1/discount(ts, .0001)
		return implied_rate(c, compounding, freq, .0001)
	else
		c = 1/discount(ts, t)
		return implied_rate(c, compounding, freq, t)
	end
end

par_rate(ts::YieldTermStructure, compounding, freq, dates::AbstractVector{CalendarTime}) =
	par_rate(ts, compounding, freq, [yearfraction(ts.dc, reference_date(ts), dt) for dt in dates])
function par_rate(ts::YieldTermStructure, compounding, freq, tm::AbstractVector)
	s = sum([discount(ts, t) for t in tm])
	r=discount(ts, tm[1])
end

type ConstantVolatilityTermStructure

end

type FlatYieldTermStructure
	dc::DayCount
	rate::Real
	compounding::Symbol
	freq::Integer
	reference_date::CalendarTime

	FlatYieldTermStructure(dc::DayCount, rate::Real, compounding::Symbol, freq::Symbol, reference_date::CalendarTime) =
			new(dc, rate, compounding, eval(freq), reference_date)
end

FlatYieldTermStructure(dc::DayCount, rate::Real) = FlatYieldTermStructure(dc, rate, :Continuous, :NoFrequency, today())
discount(ts::FlatYieldTermStructure, t::Real) = discount_factor(ts.r, ts.compounding, ts.freq, t)

const NoFrequency       = -1
const Once			 	= 0
const Annual			= 1
const Semiannual		= 2
const EveryFourthMonth  = 3
const Quarterly		 	= 4
const Bimonthly		 	= 6
const Monthly			= 12
const EveryFourthWeek  	= 13
const Biweekly		 	= 26
const Weekly			= 52
const Daily			 	= 365
const OtherFrequency   	= 999

end #Module


text:
abstract type MultiFluidModel <: EmpiricHelmholtzModel end

struct MultiFluidPropertyParam <: EoSParam
    Mw::SingleParam{Float64} #SingleParam
    rhoc::SingleParam{Float64} #SingleParam
    vc::SingleParam{Float64} #SingleParam
    Tc::SingleParam{Float64} #SingleParam
    pc::SingleParam{Float64} #SingleParam
    lb_v::SingleParam{Float64} #SingleParam
end

struct MultiFluidIdealParam <: EoSParam
    n0::SingleParam{NTuple{7, Float64}}
    theta::SingleParam{NTuple{4, Float64}}
    gamma_T::PairParam{Float64}
    gamma_v::PairParam{Float64}
    beta_T::PairParam{Float64}
    beta_v::PairParam{Float64}
end

struct MultiFluidSingleParam <: EoSParam
    n::PackedVectorSingleParam{Float64}
    t::PackedVectorSingleParam{Float64}
    d::PackedVectorSingleParam{Int}
    c::PackedVectorSingleParam{Int}
end

const FIJ_TYPE = Clapeyron.PairParameter{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, SubArray{Float64, 1, Base.ReshapedArray{Float64, 1, SparseArrays.SparseMatrixCSC{Float64, Int64}, Tuple{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64}}}, Tuple{StepRange{Int64, Int64}}, false}}

struct MultiFluidPairParam <: EoSParam
    nij::PackedSparsePairParam{Float64} #SparsePairParam #done
    tij::PackedSparsePairParam{Float64} #SparsePairParam #done
    dij::PackedSparsePairParam{Int} #SparsePairParam #done
    Fij::FIJ_TYPE #SparsePairParam #done
    beta_ij::PackedSparsePairParam{Float64} #SparsePairParam
    gamma_ij::PackedSparsePairParam{Float64} #SparsePairParam
    eta_ij::PackedSparsePairParam{Float64} #SparsePairParam
    epsilon_ij::PackedSparsePairParam{Float64} #SparsePairParam
end