text:
# Linear Interpolation
function (A::LinearInterpolation{<:AbstractVector{<:Number}})(t::Number)
  idx = findfirst(x->x>=t,A.t)-1
  idx == 0 ? idx += 1 : nothing
  θ = (t - A.t[idx])/ (A.t[idx+1] - A.t[idx])
  (1-θ)*A.u[idx] + θ*A.u[idx+1]
end

function (A::LinearInterpolation{<:AbstractMatrix{<:Number}})(t::Number)
  idx = findfirst(x->x>=t,A.t)-1
  idx == 0 ? idx += 1 : nothing
  θ = (t - A.t[idx])/ (A.t[idx+1] - A.t[idx])
  (1-θ)*A.u[:,idx] + θ*A.u[:,idx+1]
end

# Quadratic Interpolation
function (A::QuadraticInterpolation{<:AbstractVector{<:Number}})(t::Number)
  idx = findfirst(x->x>=t,A.t)-1
  idx == 0 ? idx += 1 : nothing
  if idx == length(A.t) - 1
    i₀ = idx - 1; i₁ = idx; i₂ = i₁ + 1;
  else
    i₀ = idx; i₁ = i₀ + 1; i₂ = i₁ + 1;
  end
  l₀ = ((t-A.t[i₁])*(t-A.t[i₂]))/((A.t[i₀]-A.t[i₁])*(A.t[i₀]-A.t[i₂]))
  l₁ = ((t-A.t[i₀])*(t-A.t[i₂]))/((A.t[i₁]-A.t[i₀])*(A.t[i₁]-A.t[i₂]))
  l₂ = ((t-A.t[i₀])*(t-A.t[i₁]))/((A.t[i₂]-A.t[i₀])*(A.t[i₂]-A.t[i₁]))
  A.u[i₀]*l₀ + A.u[i₁]*l₁ + A.u[i₂]*l₂
end

function (A::QuadraticInterpolation{<:AbstractMatrix{<:Number}})(t::Number)
  idx = findfirst(x->x>=t,A.t)-1
  idx == 0 ? idx += 1 : nothing
  if idx == length(A.t) - 1
    i₀ = idx - 1; i₁ = idx; i₂ = i₁ + 1;
  else
    i₀ = idx; i₁ = i₀ + 1; i₂ = i₁ + 1;
  end
  l₀ = ((t-A.t[i₁])*(t-A.t[i₂]))/((A.t[i₀]-A.t[i₁])*(A.t[i₀]-A.t[i₂]))
  l₁ = ((t-A.t[i₀])*(t-A.t[i₂]))/((A.t[i₁]-A.t[i₀])*(A.t[i₁]-A.t[i₂]))
  l₂ = ((t-A.t[i₀])*(t-A.t[i₁]))/((A.t[i₂]-A.t[i₀])*(A.t[i₂]-A.t[i₁]))
  A.u[:,i₀]*l₀ + A.u[:,i₁]*l₁ + A.u[:,i₂]*l₂
end

# Lagrange Interpolation
function (A::LagrangeInterpolation{<:AbstractVector{<:Number}})(t::Number)
  idxs = findRequiredIdxs(A,t)
  if A.t[idxs[1]] == t
    return A.u[idxs[1]]
  end
  N = zero(A.u[1]); D = zero(A.t[1]); tmp = N
  for i = 1:length(idxs)
    mult = one(A.t[1])
    for j = 1:(i-1)
      mult *= (A.t[idxs[i]] - A.t[idxs[j]])
    end
    for j = (i+1):length(idxs)
      mult *= (A.t[idxs[i]] - A.t[idxs[j]])
    end
    tmp = inv((t - A.t[idxs[i]]) * mult)
    D += tmp
    N += (tmp * A.u[idxs[i]])
  end
  N/D
end

function (A::LagrangeInterpolation{<:AbstractMatrix{<:Number}})(t::Number)
  idxs = findRequiredIdxs(A,t)
  if A.t[idxs[1]] == t
    return A.u[:,idxs[1]]
  end
  N = zero(A.u[:,1]); D = zero(A.t[1]); tmp = D
  for i = 1:length(idxs)
    mult = one(A.t[1])
    for j = 1:(i-1)
      mult *= (A.t[idxs[i]] - A.t[idxs[j]])
    end
    for j = (i+1):length(idxs)
      mult *= (A.t[idxs[i]] - A.t[idxs[j]])
    end
    tmp = inv((t - A.t[idxs[i]]) * mult)
    D += tmp
    @. N += (tmp * A.u[:,idxs[i]])
  end
  N/D
end

function (A::AkimaInterpolation{<:AbstractVector{<:Number}})(t::Number)
  i = searchsortedlast(A.t, t)
  i == 0 && return A.u[1]
  i == length(A.t) && return A.u[end]
  wj = t - A.t[i]
  @evalpoly wj A.u[i] A.b[i] A.c[i] A.d[i]
end

# ZeroSpline Interpolation
function (A::ZeroSpline{<:AbstractVector{<:Number}})(t::Number)
  if A.dir === :left
    # :left means that value to the left is used for interpolation
    i = searchsortedlast(A.t, t)
    return A.u[max(1, i)]
  else
    # :right means that value to the right is used for interpolation
    i = searchsortedfirst(A.t, t)
    return A.u[min(length(A.t), i)]
  end
end
 function (A::ZeroSpline{<:AbstractMatrix{<:Number}})(t::Number)
  if A.dir === :left
    # :left means that value to the left is used for interpolation
    i = searchsortedlast(A.t, t)
    return A.u[:, max(1, i)]
  else
    # :right means that value to the right is used for interpolation
    i = searchsortedfirst(A.t, t)
    return A.u[:, min(length(A.t), i)]
  end
end

# QuadraticSpline Interpolation
function (A::QuadraticSpline{<:AbstractVector{<:Number}})(t::Number)
  i = findfirst(x->x>=t,A.t)
  i == 1 ? i += 1 : nothing
  Cᵢ = A.u[i-1]
  σ = 1//2 * (A.z[i] - A.z[i-1])/(A.t[i] - A.t[i-1])
  A.z[i-1] * (t - A.t[i-1]) + σ * (t - A.t[i-1])^2 + Cᵢ
end

# CubicSpline Interpolation
function (A::CubicSpline{<:AbstractVector{<:Number}})(t::Number)
  i = findfirst(x->x>=t,A.t)
  i == nothing ? i = length(A.t) - 1 : i -= 1
  i == 0 ? i += 1 : nothing
  I = A.z[i] * (A.t[i+1] - t)^3 / (6A.h[i+1]) + A.z[i+1] * (t - A.t[i])^3 / (6A.h[i+1])
  C = (A.u[i+1]/A.h[i+1] - A.z[i+1]*A.h[i+1]/6)*(t - A.t[i])
  D = (A.u[i]/A.h[i+1] - A.z[i]*A.h[i+1]/6)*(A.t[i+1] - t)
  I + C + D
end

# BSpline Curve Interpolation
function (A::BSplineInterpolation{<:AbstractVector{<:Number}})(t::Number)
  # change t into param [0 1]
  idx = findfirst(x->x>=t,A.t) - 1
  idx == 0 ? idx += 1 : nothing
  t = A.p[idx] + (t - A.t[idx])/(A.t[idx+1] - A.t[idx]) * (A.p[idx+1] - A.p[idx])
  n = length(A.t)
  N = spline_coefficients(n,A.d,A.k,t)
  ucum = zero(eltype(A.u))
  for i = 1:n
    ucum += N[i] * A.c[i]
  end
  ucum
end

# BSpline Curve Approx
function (A::BSplineApprox{<:AbstractVector{<:Number}})(t::Number)
  # change t into param [0 1]
  idx = findfirst(x->x>=t,A.t) - 1
  idx == 0 ? idx += 1 : nothing
  t = A.p[idx] + (t - A.t[idx])/(A.t[idx+1] - A.t[idx]) * (A.p[idx+1] - A.p[idx])
  n = length(A.t)
  N = spline_coefficients(A.h,A.d,A.k,t)
  ucum = zero(eltype(A.u))
  for i = 1:A.h
    ucum += N[i] * A.c[i]
  end
  ucum
end

# Loess
function (A::Loess{<:AbstractVector{<:Number}})(t::Number)
  tmp = sort(abs.(A.t .- t))
  w = abs.(A.t .- t) ./ tmp[A.q]
  for i = 1:length(A.t)
    if w[i] <= one(A.t[1])
      w[i] = (1 - (w[i] ^ 3)) ^ 3
    else
      w[i] = zero(A.t[1])
    end
  end
  w = Diagonal(w)
  b = inv(transpose(A.x) * w * A.x) * transpose(A.x) * w * A.u
  u = zero(t[1])
  for (idx,v) in enumerate(b)
    u += v*(t^(idx-1))
  end
  u
end

# GaussianProcess
function (A::GPInterpolation{<:AbstractVector{<:Number}})(t::AbstractVector{<:Number})
  rand(A.gp,t,1)
end

function (A::GPInterpolation{<:AbstractVector{<:Number}})(t::Number)
  rand(A.gp,[t],1)[1]
end

# Curvefit
function (A::CurvefitCache{<:AbstractVector{<:Number}})(t::Union{AbstractVector{<:Number},Number})
  A.m(t,A.pmin)
end


text:
import AbstractTrees
using FlameGraphs

using FlameGraphs: Node, NodeData
using Base.StackTraces: StackFrame
using Profile: StackFrameTree

function pprof(fg::Node{NodeData},
    period::Union{Nothing, UInt64} = nothing;
    web::Bool = true,
    webhost::AbstractString = "localhost",
    webport::Integer = 57599,
    out::AbstractString = "profile.pb.gz",
    from_c::Bool = true,
    full_signatures::Bool = true,
    drop_frames::Union{Nothing, AbstractString} = nothing,
    keep_frames::Union{Nothing, AbstractString} = nothing,
    ui_relative_percentages::Bool = true,
)
    if period === nothing
        period = ccall(:jl_profile_delay_nsec, UInt64, ())
    end
    @assert !isempty(basename(out)) "`out=` must specify a file path to write to. Got unexpected: '$out'"
    if !endswith(out, ".pb.gz")
        out = "$out.pb.gz"
        @info "Writing output to $out"
    end

    string_table = OrderedDict{AbstractString, Int64}()
    enter!(string) = _enter!(string_table, string)
    enter!(::Nothing) = _enter!(string_table, "nothing")
    ValueType!(_type, unit) = ValueType(_type = enter!(_type), unit = enter!(unit))

    function _register_function(funcs, id, linfo, frame)
        # Known function
        haskey(funcs, id) && return

        # Store the function in our functions dict
        funcProto = Function()
        funcProto.id = id
        file = nothing
        simple_name = _escape_name_for_pprof(frame.func)
        local full_name_with_args
        if frame.linfo !== nothing && frame.linfo isa Core.MethodInstance
            linfo = frame.linfo::Core.MethodInstance
            meth = linfo.def
            file = string(meth.file)
            io = IOBuffer()
            Base.show_tuple_as_call(io, meth.name, linfo.specTypes)
            name = String(take!(io))
            full_name_with_args = _escape_name_for_pprof(name)
            funcProto.start_line = convert(Int64, meth.line)
        else
            # frame.linfo either nothing or CodeInfo, either way fallback
            file = string(frame.file)
            full_name_with_args = _escape_name_for_pprof(string(frame.func))
            funcProto.start_line = convert(Int64, frame.line) # TODO: Get start_line properly
        end
        # WEIRD TRICK: By entering a *different value* for the name and system_name, pprof
        # will use the supplied `name` *verbatim*, without pruning off the arguments. So
        # even when full_signatures == false, we want to generate two `enter!` ids. So we
        # achieve that by entering an _empty_ name for the system_name, so that pprof will
        # use the `name` as provided.
        funcProto.system_name = enter!("")
        if full_signatures
            funcProto.name = enter!(full_name_with_args)
        else
            funcProto.name = enter!(simple_name)
        end
        file = Base.find_source_file(file)
        funcProto.filename   = enter!(file)
        # Only keep C functions if from_c=true
        if (from_c || !frame.from_c)
            funcs[id] = funcProto
        end
    end

    # Setup:
    enter!("")  # NOTE: pprof requires first entry to be ""
    # Functions need a uid, we'll use the pointer for the method instance
    seen_funcs = Set{UInt64}()
    funcs = Dict{UInt64, Function}()

    seen_locs = Set{UInt64}()
    locs  = Dict{UInt64, Location}()
    locs_from_c  = Dict{UInt64, Bool}()

    sample_type = [
        ValueType!("events", "count"), # Mandatory
    ]

    prof = PProfile(
        sample = [], location = [], _function = [],
        mapping = [], string_table = [],
        sample_type = sample_type, default_sample_type = 1, # events
        period = period, period_type = ValueType!("cpu", "nanoseconds")
    )

    if drop_frames !== nothing
        prof.drop_frames = enter!(drop_frames)
    end
    if keep_frames !== nothing
        prof.keep_frames = enter!(keep_frames)
    end

    # start decoding backtraces
    lastwaszero = true

    # Do a top-down walk of the flame graph from root to leaves.
    # At each node, we enter a "sample" to flame graph if there is any time not covered by
    # its children. That is, if the exclusive time of any node is non-zero.
    # AND, we want to enter each "section" of exclusive time as a separate sample, in case
    # that allows PProf to do any ordering in its display.
    # So, the algorithm is that we enter a sample for every gap between the spans of the
    # children.
    function emit_tree(node)
        span = node.data.span
        start = span.start
        stop = span.stop
        child = node.child
        if child === node
            # We're a leaf node, so just emit this node directly, and then we're done.
            emit_stack_sample(child, start:stop)
            return
        end
        while true  # break when we've exhausted the children
            cspan = child.data.span
            if cspan.start > start
                # We know there's a gap before the next child starts
                # So we emit a frame for `node`
                emit_stack_sample(node, start:cspan.start-1)
            end
            # Now emit the child node
            emit_tree(child)

            start = cspan.stop + 1

            if child.sibling === child break end
            child = child.sibling
        end
        # If there was still time after the last child, that's another gap that we can emit.
        cspan = child.data.span
        if stop > cspan.stop
            # We know there's a gap before the next child starts
            # So we emit a frame for `node`
            emit_stack_sample(node, start:stop)
        end
    end

    function emit_stack_sample(leaf, span)
        location_id = Vector{Any}()

        # Walk back up the callstack, collecting all stack traces along the way, building
        # up a stack trace for this span.
        node = leaf
        while true  # do-while node.parent === node    -- this makes sure we get the top node.
            data = node.data

            if from_c || !data.sf.from_c
                frame = data.sf
                # Don't use the ip pointer, because this data is provided by the user and isn't
                # guaranteed to be unique.
                id = hash(frame)

                location = Location(;id = id, address = frame.pointer, line=[])
                push!(location_id, id)
                locs[id] = location
                linfo = data.sf.linfo

                # If the function comes with a `linfo` (meaning the profile was generated
                # in this process), we include it in the hash. If not, it's will be
                # `nothing`, which is okay to hash.
                # We include all of these in the hash though because in some contexts (such
                # as `@snoopi_deep` profiles), two frames with the same `linfo` might have
                # different func name.
                func_id = hash((frame.func, frame.file, frame.line, frame.linfo))

                _register_function(funcs, func_id, linfo, frame)
                push!(location.line, Line(function_id = func_id, line = frame.line))
            end

            if node.parent === node
                break
            end

            node = node.parent
        end

        value = [
            length(span), # Number of samples in this frame.
        ]
        push!(prof.sample, Sample(;location_id = location_id, value = value))

    end

    emit_tree(fg)

    # Build Profile
    prof.string_table = collect(keys(string_table))
    # If from_c=false funcs and locs should NOT contain C functions
    prof._function = collect(values(funcs))
    prof.location  = collect(values(locs))

    # Write to disk
    open(out, "w") do io
        writeproto(io, prof)
    end

    if web
        refresh(webhost = webhost, webport = webport, file = out,
            ui_relative_percentages = ui_relative_percentages)
    end

    out
end

ext: jl


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie
let
    xs = LinRange(-3, 3, 20)
    ys = LinRange(-3, 3, 20)
    us = [x + y for x in xs, y in ys]
    vs = [y - x for x in xs, y in ys]
    strength = vec(sqrt.(us .^2 .+ vs .^2))

    fig = Figure(resolution = (550, 450), fontsize = 18, font = "sans")
    ax = Axis(fig, xlabel = "x", ylabel = "y", aspect = 1)
    arrows!(ax, xs, ys, us, vs, arrowsize = 10, lengthscale = 0.1,
        arrowcolor = strength, linecolor = strength, colormap = :plasma)
    cbar = Colorbar(fig, limits =(minimum(strength),maximum(strength)),
    nsteps =100, colormap = :plasma, ticksize=15, width = 15, tickalign=1)
    limits!(ax, -3,3,-3,3)
    fig[1,1] = ax
    fig[1,2] = cbar
    fig
    #save("arrows.png", fig, px_per_unit = 2)
    save(joinpath(@__DIR__, "output", "arrows.png"), fig, px_per_unit = 2) # HIDE
end
using Pkg # HIDE
Pkg.status("CairoMakie") # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
function f1(x1, y1=1, z1=2; t1=3, u1=5)
    @show x1, y1, z1, t1, u1
    nothing
end

function f2(x2, y2=1, z2=2; t2=3, u2=5)
    @show x2, y2, z2, t2, u2
    nothing
end

function f3(x3, f1args = tuple(), f2args = tuple(); f1kwargs = (;), f2kwargs=(;))
    f1(x3, f1args...; f1kwargs...)

    f2(x3, f2args...; f2kwargs...)
end

f3(5)
f3(5, (2, 3))
f3(5, (4, 10), (6, 0))
f3(5, (4, 10), (6, 0), f1kwargs=(t1 = 5,))
f3(5, (4, 10), (6, 0), f1kwargs=(t1 = 5, u1 = 10))
f3(5, (4, 10), (6, 0), f1kwargs=(t1 = 5, u1 = 10), f2kwargs=(t2 = 5,))
f3(5, (4, 10), (6, 0), f1kwargs=(t1 = 5, u1 = 10), f2kwargs=(t2 = 5, u2 = 10))


text:
#
# Handles CR3BP problem states.
#

"""
    `const CR3BPFrames = Union{Synodic, Inertial}`

Coordinate frames that are valid within the Circular Restricted Three-body Problem.
"""
const CR3BPFrames = Union{Synodic, Inertial}

"""
    `const NormalizedCartesianState{F, FR<:CR3BPFrames} = CartesianState{F, NormalizedLengthUnit, NormalizedTimeUnit, FR}`

CR3BP systems often use normalized units for calculations.
This `NormalizedCartesianState` represents a Cartesian
state with no dimensioned units.
"""
const NormalizedCartesianState{F, FR<:CR3BPFrames} = CartesianState{F, NormalizedLengthUnit, NormalizedTimeUnit, FR}

"""
    `NormalizedCartesianState(r, v, t, frame)`

Convenience constructor for a normalized `CartesianState`.
"""
NormalizedCartesianState(r, v, t=0, frame=Inertial) = CartesianState(r, v, t, frame; lengthunit = NormalizedLengthUnit(), timeunit = NormalizedTimeUnit())

"""
Print `NormalizedCartesianState` instances to `io`.
"""
function Base.show(io::IO, state::NormalizedCartesianState{F,FR}) where {F,FR}
    println(io, "  ", "Normalized ", string(FR), " Cartesian State:")
    println("")
    println(io, "    t:  ", state.t)
    println(io, "    r = ", [state.r[1] state.r[2] state.r[3]])
    println(io, "    v = ", [state.v[1] state.v[2] state.v[3]])
end

"""
Print `NormalizedCartesianState` instances to `io`.
"""
Base.show(io::IO, ::MIME"text/plain", state::NormalizedCartesianState{F,FR}) where {F,FR} = show(io, state)

"""
A structure which contains the linearization of CR3BP dynamics at timepoint `t`.
"""
mutable struct SynodicCartesianSTMState{F, LU, TU} <: AbstractState{F, LU, TU, Synodic}
    cart::CartesianState{F, LU, TU, Synodic}
    stm::MMatrix{6,6,F}

    function SynodicCartesianSTMState(r::AbstractVector{<:Real}, v::AbstractVector{<:Real}, t::Real = 0, stm::AbstractMatrix{<:Real} = I(6);
                                      lengthunit = NormalizedLengthUnit(), timeunit = NormalizedTimeUnit())
        state = CartesianState(r, v, t, Synodic; lengthunit = lengthunit, timeunit = timeunit)
        F = eltype(state)
        LU = typeof(AstrodynamicsCore.lengthunit(state))
        TU = typeof(AstrodynamicsCore.timeunit(state))

        return new{F, LU, TU}(state, SMatrix{6,6,eltype(state)}(stm))
    end

    function SynodicCartesianSTMState(r::AbstractVector{Unitful.Length}, v::AbstractVector{Unitful.Time}, t::Unitful.Time = 0u"s", stm::AbstractMatrix{<:Real} = I(6))
        state = CartesianState(r, v, t, Synodic)
        return new{eltype(state), typeof(lengthunit(state)), typeof(timeunit(state))}(state, SMatrix{6,6,eltype(state)}(stm))
    end

    function SynodicCartesianSTMState(state::SynodicCartesianState{F,LU,TU}, stm = SMatrix{6,6,F}(I(6))) where {F,LU,TU}
        return new{F, LU, TU}(state, stm)
    end

end

# Overrides `Unitful` unit conversions for `AbstractUnitfulStructure` instances.
(u::Unitful.LengthFreeUnits)(state::SynodicCartesianSTMState{F,LU,TU}) where {F,LU,TU} = convert(SynodicCartesianSTMState{F, typeof(u), TU}, state)

# Overrides `Unitful` unit conversions for `AbstractUnitfulStructure` instances.
(u::Unitful.TimeFreeUnits)(state::SynodicCartesianSTMState{F,LU,TU}) where {F,LU,TU} = convert(SynodicCartesianSTMState{F, LU, typeof(u)}, state)

"""
Print `SynodicCartesianSTMState` instances to `io`.
"""
function Base.show(io::IO, state::SynodicCartesianSTMState{F,LU, TU}) where {F,LU,TU}
    println(io, state.cart,"\n")
    println(io, "  Local Linearization (Φ):\n")
    for row ∈ eachrow(state.stm)
        print(io, "    ")
        for el ∈ row
            print(io, el, " ")
        end
        println(io,"")
    end
end

"""
Print `SynodicCartesianSTMState` instances to `io`.
"""
Base.show(io::IO, ::MIME"text/plain", state::SynodicCartesianSTMState{F,LU,TU}) where {F,LU,TU} = show(io, state)

"""
Returns the timepoint associated with a `CartesianState`.
"""
epoch(state::SynodicCartesianSTMState) = epoch(state.cart)

"""
Returns the `Unitful` position vector of the `Cartesianstate`.
"""
position_vector(state::SynodicCartesianSTMState) = position_vector(state.cart)

"""
Returns the `Unitful` velocity vector of the `CartesianState`.
"""
velocity_vector(state::SynodicCartesianSTMState) = velocity_vector(state.cart)

"""
Returns the `Unitful` scalar position of the `CartesianState`.
"""
scalar_position(state::SynodicCartesianSTMState) = scalar_position(state.cart)

"""
Returns the `Unitful` scalar velocity of the `CartesianState`.
"""
scalar_velocity(state::SynodicCartesianSTMState) = scalar_velocity(state.cart)

"""
    `NormalizedCartesianState(r::AbstractVector{<:Real}, v::AbstractVector{<:Real}, frame::FR = Synodic, epoch::Real = 0) where FR <: CR3BPFrames = CartesianState(r, v, frame, epoch)`

Constructor for a `NormalizedCartesianState`. No `Unitful` arguments allowed!
"""
NormalizedCartesianState(r::AbstractVector{<:Real}, v::AbstractVector{<:Real}, frame::FR = Synodic, epoch::Real = 0) where FR <: CR3BPFrames = CartesianState(r, v, frame, epoch)

"""
    `const NormalizedSynodicCartesianSTMState{F, FR<:CR3BPFrames} = CartesianState{F, NormalizedLengthUnit, NormalizedTimeUnit, FR}`

CR3BP systems often use normalized units for calculations.
This `NormalizedCartesianState` represents a Cartesian
state with no dimensioned units.
"""
const NormalizedSynodicCartesianSTMState{F} = SynodicCartesianSTMState{F, NormalizedLengthUnit, NormalizedTimeUnit}

"""
    `MinimalCircularRestrictedThreeBodySystem(μ::Real) <: AbstractSystem{F, NormalizedLengthUnit, NormalizedTimeUnit}`

Sometimes you want to do CR3BP calculations with only one system parameter,
the nondimensional mass parameter `μ`. That's what this structure is for!
"""
mutable struct MinimalCircularRestrictedThreeBodySystem{F} <: AbstractSystem{F, NormalizedLengthUnit, NormalizedTimeUnit}
    μ::F
    name::String

    function MinimalCircularRestrictedThreeBodySystem(μ::Real, name = "")
        @assert μ ≤ 1//2 "Nondimensional mass parameter must be less than 0.5 by definition."
        F = typeof(μ)
        return new{F}(F(μ), name)
    end
end

"""
Returns the normalized (nondimensional) mass parameter for the CR3BP system.
"""
function normalized_mass_parameter(state::MinimalCircularRestrictedThreeBodySystem)
    return state.μ
end

function mass_parameter(::MinimalCircularRestrictedThreeBodySystem)
    throw(ArgumentError("Only a normalized mass parameter can be computed from a MinimalCircularRestrictedThreeBodySystem. Try `normalized_mass_parameter(sys)` instead."))
end

"""
Print `MinimalCircularRestrictedThreeBodySystem` instances to `io`.
"""
function Base.show(io::IO, sys::MinimalCircularRestrictedThreeBodySystem{F}) where {F}
    println(io, "  Normalized Circular Restricted Three-body System", isempty(sys.name) ? "\n" : " ($(sys.name))\n")
    println(io, "    μ:    ", string(normalized_mass_parameter(sys)))
end

"""
Print `MinimalCircularRestrictedThreeBodySystem` instances to `io`.
"""
Base.show(io::IO, ::MIME"text/plain", sys::MinimalCircularRestrictedThreeBodySystem{F}) where {F} = show(io, sys)

"""
    `mutable struct CircularRestrictedThreeBodySystem{F, LU, TU} <: AbstractSystem{F, LU, TU}`

This structure represents _all_ required parameters that are used to describe
a Circular Restricted Three-body System: normalized length unit `DU`,
normalized time unit `DT`, and a `Tuple` of two mass parameters, `μ`.
"""
mutable struct CircularRestrictedThreeBodySystem{F, LU, TU} <: AbstractSystem{F, LU, TU}
    DU::F
    DT::F
    μ::Tuple{F,F}
    name::String

    function CircularRestrictedThreeBodySystem(μ::Tuple{Real, Real}, DU::Real, DT::Real = time_scale_factor(DU, μ...), name = ""; lengthunit = u"km", timeunit = u"s")
        F = promote_type(typeof.(μ)..., typeof(DU), typeof(DT))
        if !(F <: AbstractFloat)
            @warn "Promoted type $(string(F)) is not a float. Defaulting to Float64."
            F = Float64
        end
        return new{F, typeof(lengthunit), typeof(timeunit)}(F(DU), F(DT), F.(μ), name)
    end

    function CircularRestrictedThreeBodySystem(μ::Tuple{MassParameter,MassParameter}, DU::Unitful.Length, DT::Unitful.Time = time_scale_factor(DU, μ...), name = "")
        lengthunit = unit(DU)
        timeunit = unit(DT)
        return CircularRestrictedThreeBodySystem(ustrip.(lengthunit^3 / timeunit^2, μ),
                                                 ustrip(lengthunit, DU), ustrip(timeunit, DT), name;
                                                 lengthunit = lengthunit, timeunit = timeunit)
    end
end

"""
Convert between `eltype`, `lengthunit`, and `timeunit` types for CR3BP systems.
"""
function Base.convert(::Type{CircularRestrictedThreeBodySystem{F, LU, TU}}, sys::CircularRestrictedThreeBodySystem) where {F, LU, TU}
    return CircularRestrictedThreeBodySystem(
        F.(ustrip.(LU()^3 / TU()^2, mass_parameters(sys))),
        F.(ustrip.(LU(), normalized_length_unit(sys))),
        F.(ustrip.(TU(), normalized_time_unit(sys))),
        sys.name; lengthunit = LU(), timeunit = TU()
    )
end

"""
Returns the normalized length unit for a CR3BP system.
"""
normalized_length_unit(sys::CircularRestrictedThreeBodySystem) = sys.DU * lengthunit(sys)

"""
Returns the normalized time unit for a CR3BP system.
"""
normalized_time_unit(sys::CircularRestrictedThreeBodySystem) = sys.DT * timeunit(sys)

"""
Print `CircularRestrictedThreeBodySystem` instances to `io`.
"""
function Base.show(io::IO, sys::CircularRestrictedThreeBodySystem{F,LU,TU}) where {F,LU,TU}
    println(io, "  Circular Restricted Three-body System", isempty(sys.name) ? "\n" : " ($(sys.name))\n")
    println(io, "        Length Unit:    ", sys.DU, " ", string(lengthunit(sys)))
    println(io, "          Time Unit:    ", sys.DT, " ", string(timeunit(sys)))
    println(io, "    Mass Parameters:    ", "(", max(sys.μ...), ", ", min(sys.μ...), ") ", string(massparameterunit(sys)))
end

"""
Print `CircularRestrictedThreeBodySystem` instances to `io`.
"""
Base.show(io::IO, ::MIME"text/plain", sys::CircularRestrictedThreeBodySystem{F,LU,TU}) where {F,LU,TU} = show(io, sys)

"""
Returns the normalized (nondimensional) mass parameter for a CR3BP system.
"""
normalized_mass_parameter(sys::CircularRestrictedThreeBodySystem) = min(sys.μ...) / reduce(+, sys.μ)

"""
Returns the dimensioned (not normalized) mass parameters for a CR3BP system.
"""
mass_parameters(sys::CircularRestrictedThreeBodySystem) = sys.μ .* massparameterunit(sys)

"""
Returns the primary body's mass parameter within a CR3BP system.
"""
primary_mass_parameter(sys::CircularRestrictedThreeBodySystem) = max(sys.μ...) * massparameterunit(sys)

"""
Returns the secondary body's mass parameter within a CR3BP system.
"""
secondary_mass_parameter(sys::CircularRestrictedThreeBodySystem) = min(sys.μ...) * massparameterunit(sys)

"""
A `Union` of all valid `CR3BP` state types.
"""
const CR3BPState{F, LU, TU, FR<:CR3BPFrames} = Union{SynodicCartesianSTMState{F, LU, TU}, CartesianState{F, LU, TU, FR}}

"""
Contains __all__ relevant state and system information within the CR3BP.
"""
mutable struct CircularRestrictedThreeBodyOrbit{
        F, LU, TU,
        C<:CR3BPState{F, <:Unitful.LengthFreeUnits, <:Unitful.TimeFreeUnits, <:CR3BPFrames},
        S<:Union{MinimalCircularRestrictedThreeBodySystem{F}, CircularRestrictedThreeBodySystem{F,LU,TU}}
    } <: AbstractOrbit{F, LU, TU}

    state::C
    system::S

    function CircularRestrictedThreeBodyOrbit(
            r::AbstractVector{<:Real}, v::AbstractVector{<:Real},
            μ::Tuple{<:Real, <:Real}, DU::Real, DT::Real, t::Real = 0;
            frame = Synodic, lengthunit = u"km", timeunit = u"s")

        F = promote_type(eltype(r), eltype(v), typeof(μ)..., typeof(DU), typeof(DT), typeof(t))
        if !(F <: AbstractFloat)
            @warn "Non-float promoted type $(string(F)) provided. Defaulting to Float64."
            F = Float64
        end

        state = CartesianState(F.(r), F.(v), F(t), frame; lengthunit = lengthunit, timeunit = timeunit)
        system = CircularRestrictedThreeBodySystem(
            F.(ustrip.(lengthunit(state)^3 / timeunit(state)^2, μ)),
            F(ustrip(lengthunit(state), DU)),
            F(ustrip(timeunit(state), DT)))

        return new{F, typeof(lengthunit(state)), typeof(timeunit(state)), typeof(state), typeof(system)}(state, system)

    end

    function CircularRestrictedThreeBodyOrbit(
            r::AbstractVector{Unitful.Length}, v::AbstractVector{Unitful.Length},
            μ::Tuple{MassParameter, MassParameter},
            DU::Unitful.Length, t::Unitful.Time = 0u"s";
            frame = Synodic, name = "")

        F = promote_type(eltype(ustrip.(r)), eltype(ustrip.(v)), typeof(ustrip.(μ))..., typeof(ustrip(DU)), typeof(ustrip(t)))
        if !(F <: AbstractFloat)
            @warn "Non-float promoted type $(string(F)) provided. Defaulting to Float64."
            F = Float64
        end
        state = CartesianState(F.(r), F.(v), F(t), frame)
        DT = time_scale_factor(DU, μ...)
        system = CircularRestrictedThreeBodySystem(
            F.(ustrip.(lengthunit(state)^3 / timeunit(state)^2, μ)),
            F(ustrip(lengthunit(state), DU)),
            F(ustrip(timeunit(state), DT)),
            name)

        return new{F, typeof(lengthunit(state)), typeof(timeunit(state)), typeof(state), typeof(system)}(state, system)
    end

    function CircularRestrictedThreeBodyOrbit(r::AbstractVector{<:Real}, v::AbstractVector{<:Real}, μ::Real, t::Real = 0; frame = Synodic, name = "")
        F = promote_type(eltype(r), eltype(v), typeof(μ), typeof(t))
        if !(F <: AbstractFloat)
            @warn "Non-float promoted type $(string(F)) provided. Defaulting to Float64."
            F = Float64
        end
        state = CartesianState(F.(r), F.(v), F(t), frame; lengthunit = NormalizedLengthUnit(), timeunit = NormalizedTimeUnit())
        system = MinimalCircularRestrictedThreeBodySystem(F(μ), name)
        return new{F, NormalizedLengthUnit, NormalizedTimeUnit, typeof(state), typeof(system)}(state, system)
    end

    function CircularRestrictedThreeBodyOrbit(r::AbstractVector{<:Real}, v::AbstractVector{<:Real}, system::CircularRestrictedThreeBodySystem, t::Real = 0; frame = Synodic)
        state = CartesianState(r, v, t, frame; lengthunit = AstrodynamicsCore.lengthunit(system), timeunit = AstrodynamicsCore.timeunit(system))
        F = promote_type(eltype(state), eltype(system))
        L = lengthunit(system) |> typeof
        T = timeunit(system)   |> typeof

        state  = convert(CartesianState{F}, state)
        system = convert(CircularRestrictedThreeBodySystem{F, L, T}, system)

        return new{F, L, T, typeof(state), typeof(system)}(state, system)
    end

    function CircularRestrictedThreeBodyOrbit(r::AbstractVector{<:Unitful.Length}, v::AbstractVector{<:Unitful.Velocity}, system::CircularRestrictedThreeBodySystem, t::S = 0u"s"; frame = Synodic) where S <: Unitful.Time
        state = CartesianState(r, v, t, frame)
        F = promote_type(eltype(state), eltype(system))
        L = lengthunit(system) |> typeof
        T = timeunit(system)   |> typeof

        state  = convert(CartesianState{F, L, T}, state)
        system = convert(CircularRestrictedThreeBodySystem{F, L, T}, system)

        return new{F, L, T, typeof(state), typeof(system)}(state, system)
    end

    function CircularRestrictedThreeBodyOrbit(r::AbstractVector{<:Real}, v::AbstractVector{<:Real}, system::MinimalCircularRestrictedThreeBodySystem, t::Real = 0; frame = Synodic)
        F = promote_type(eltype(r), eltype(v), eltype(system), typeof(t))
        state = NormalizedCartesianState(F.(r), F.(v), F(t), frame)
        system = MinimalCircularRestrictedThreeBodySystem(F(normalized_mass_parameter(system)))
        return new{F, NormalizedLengthUnit, NormalizedTimeUnit, typeof(state), typeof(system)}(state, system)
    end

    function CircularRestrictedThreeBodyOrbit(state::NormalizedCartesianState, system::MinimalCircularRestrictedThreeBodySystem, t::Real = 0; frame = Synodic)
        F = promote_type(typeof(state).parameters[1], typeof(system).parameters[1])
        return new{F, NormalizedLengthUnit, NormalizedTimeUnit, typeof(state), typeof(system)}(state, system)
    end

    function CircularRestrictedThreeBodyOrbit(state::CartesianState, sys::CircularRestrictedThreeBodySystem)
        F = promote_type(eltype(state), eltype(sys))
        state = convert(CartesianState{F, typeof(lengthunit(state)), typeof(timeunit(state))}, state)
        return new{F, typeof(lengthunit(sys)), typeof(timeunit(sys)), typeof(state), typeof(sys)}(state, sys)
    end

    function CircularRestrictedThreeBodyOrbit(state::SynodicCartesianSTMState, sys::CircularRestrictedThreeBodySystem)
        F = promote_type(eltype(state), eltype(sys))
        state = convert(SynodicCartesianSTMState{F, typeof(lengthunit(state)), typeof(timeunit(state))}, state)
        return new{F, typeof(lengthunit(sys)), typeof(timeunit(sys)), typeof(state), typeof(sys)}(state, sys)
    end
end

"""
A simpler constructor for a `Synodic` CR3BP state _with_ or
_without_ the local linearization (state transition matrix)
included.

Please note that this is a _convenience_ constructor for quick calculations.
This method is __not__ type stable, since the `kwarg` `include_stm` directly
affects the output type of the function. This sub-optimal code is great for
quickly constructing an initial state, but this method __will__ be slower
than an explicit `CircularRestrictedThreeBodyOrbit` constructor.
"""
Orbit(r, v, sys::CircularRestrictedThreeBodySystem; include_stm = false) = include_stm ? CircularRestrictedThreeBodyOrbit(SynodicCartesianSTMState(r, v), sys) : CircularRestrictedThreeBodyOrbit(r, v, sys)

"""
Returns the position vector of the CR3BP orbit.
"""
position_vector(orb::CircularRestrictedThreeBodyOrbit) = position_vector(orb.state)

"""
Returns the velocity vector of the CR3BP orbit.
"""
velocity_vector(orb::CircularRestrictedThreeBodyOrbit) = velocity_vector(orb.state)

"""
Returns the epoch associated with the `state` field.
"""
epoch(orb::CircularRestrictedThreeBodyOrbit) = epoch(orb.state)

"""
Convert between `eltype`, `lengthunit`, and `timeunit` types for CR3BP AstrodynamicsCore.
"""
function Base.convert(::Type{CircularRestrictedThreeBodyOrbit{F, LU, TU}}, orb::CircularRestrictedThreeBodyOrbit) where {F, LU, TU}
    return CircularRestrictedThreeBodyOrbit(
        F.(position_vector(orb.state)),
        F.(velocity_vector(orb.state)),
        convert(CircularRestrictedThreeBodySystem{F, LU, TU}, orb.system),
        F(epoch(orb.state));
        frame = coordinateframe(orb.state)
    )
end

"""
An alias for `CircularRestrictedThreeBodyOrbit` instances with `NormalizedCartesianState` states.
"""
const NormalizedCR3BPOrbit{F,LU,TU,ST<:NormalizedCartesianState{F,<:AbstractFrame},SR} = CircularRestrictedThreeBodyOrbit{F, LU, TU, ST, SR}

"""
An alias for `CircularRestrictedThreeBodyOrbit` instances with `NormalizedSynodicCartesianSTMState` states.
"""
const NormalizedSynodicSTMCR3BPOrbit{F,LU,TU,ST<:NormalizedSynodicCartesianSTMState{F},SR} = CircularRestrictedThreeBodyOrbit{F, LU, TU, ST, SR}

"""
A constructor for a Normalized Synodic CR3BP orbit with a local linearization.
"""
function NormalizedSynodicSTMCR3BPOrbit(orbit::CircularRestrictedThreeBodyOrbit)
    return CircularRestrictedThreeBodyOrbit(SynodicCartesianSTMState((normalize ∘ synodic)(orbit).state), orbit.system)
end

"""
An alias for `CircularRestrictedThreeBodyOrbit` instances with `NormalizedCartesianState` states in the `Synodic` frame.
"""
const NormalizedSynodicCR3BPOrbit{F,LU,TU,ST<:NormalizedCartesianState{F,Synodic}, SR} = NormalizedCR3BPOrbit{F, LU, TU, ST,SR}

"""
Returns the position of one of the two CR3BP bodies in the `state`'s reference frame.
"""
function body_position(orbit::CircularRestrictedThreeBodyOrbit, N)
    (N == 1 || N == 2) || throw(ArgumentError("The second argument must specify the primary or secondary body! Set N equal to 1 or 2."))

    if N == 1
        pos = MVector{3}(-normalized_mass_parameter(orbit.system), 0.0, 0.0)
        if coordinateframe(orbit.state) != Synodic
            pos = inertial(pos, epoch(normalize(orbit)))
        end
        pos *= lengthunit(orbit.state)
    else
        pos = MVector{3}(1-normalized_mass_parameter(orbit.system), 0.0, 0.0)
        if coordinateframe(orbit.state) != Synodic
            pos = inertial(pos, epoch(normalize(orbit)))
        end
        pos *= lengthunit(orbit.state)
        return pos
    end
end

"""
Returns the position of one of the two CR3BP bodies in the normalized `Synodic` reference frame.
"""
function body_position(sys::CircularRestrictedThreeBodySystem, N)
    (N == 1 || N == 2) || throw(ArgumentError("The second argument must specify the primary or secondary body! Set N equal to 1 or 2."))

    if N == 1
        pos = MVector{3}(-normalized_mass_parameter(sys), 0.0, 0.0)
    else
        pos = MVector{3}(1-normalized_mass_parameter(sys), 0.0, 0.0)
        return pos
    end
end

"""
Returns the position of the orbiter with respect to one of the two CR3BP bodies,
in the `state`'s reference frame.
"""
function position_wrt_body(orbit::CircularRestrictedThreeBodyOrbit, N)
    (N == 1 || N == 2) || throw(ArgumentError("The second argument must specify the primary or secondary body! Set N equal to 1 or 2."))
    body_pos = body_position(orbit, N)
    return position_vector(orbit) .- body_pos
end

"""
An alias for `CircularRestrictedThreeBodyOrbit`.
"""
const CR3BPOrbit = CircularRestrictedThreeBodyOrbit

"""
an alias for `MinimalCircularRestrictedThreeBodySystem`.
"""
const CR3BPSystem = CircularRestrictedThreeBodySystem

"""
Returns the normalized (nondimensional) length unit associated with a CR3BP orbit.
"""
normalized_length_unit(orb::CircularRestrictedThreeBodyOrbit) = normalized_length_unit(orb.system)

"""
Returns the normalized (nondimensional) time unit associated with a CR3BP orbit.
"""
normalized_time_unit(orb::CircularRestrictedThreeBodyOrbit) = normalized_time_unit(orb.system)

"""
Returns the normalized (nondimensional) mass parameter associated with a CR3BP orbit.
"""
normalized_mass_parameter(orb::CircularRestrictedThreeBodyOrbit) = normalized_mass_parameter(orb.system)

"""
Returns the dimensioned (not normalized) mass parameters associated with a CR3BP orbit.
"""
mass_parameters(orb::CircularRestrictedThreeBodyOrbit) = mass_parameters(orb.system)

"""
Returns the mass parameter of the primary body within a CR3BP orbit.
"""
primary_mass_parameter(orb::CircularRestrictedThreeBodyOrbit) = primary_mass_parameter(orb.system)

"""
Returns the mass parameter of the primary body within a CR3BP orbit.
"""
secondary_mass_parameter(orb::CircularRestrictedThreeBodyOrbit) = secondary_mass_parameter(orb.system)

"""
Prints a `CircularRestrictedThreeBodyOrbit` to `io`.
"""
function Base.show(io::IO, orb::CircularRestrictedThreeBodyOrbit{F,LU,TU}) where {F,LU,TU}
    println(io, "Circular Restricted Three-body Orbit\n")
    println(io, orb.state)
    println(io, orb.system)
end

"""
Prints a `CircularRestrictedThreeBodyOrbit` to `io`.
"""
Base.show(io::IO, ::MIME"text/plain", orb::CircularRestrictedThreeBodyOrbit{F,LU,TU}) where {F,LU,TU} = show(io, orb)

"""
Prints a `Trajectory` instance to `io`.
"""
Base.show(io::IO, traj::Trajectory{<:CircularRestrictedThreeBodyOrbit}) = println(io, "Circular Restricted Three-body trajectory with ", length(traj), " steps")

"""
Prints a `Manifold` instance to `io`.
"""
Base.show(io::IO, man::Manifold{<:Trajectory{<:CircularRestrictedThreeBodyOrbit}}) = println(io, "Circular Restricted Three-body manifold with ", length(man), " trajectories")

"""
Returns true if both systems are equal.
"""
function Base.isequal(c1::CircularRestrictedThreeBodySystem, c2::CircularRestrictedThreeBodySystem)
    return normalized_length_unit(c1)    == normalized_length_unit(c2) &&
           normalized_time_unit(c1)      == normalized_time_unit(c2)   &&
           normalized_mass_parameter(c1) == normalized_mass_parameter(c2)
end


text:
export stoppropagation,
       subscribe,
       samplesignals

# Don't allow a signal to propagate outward
immutable StopPropagation <: Tile
    tile::Tile
    name::Symbol
end

@doc """
Stop a UI signal from propagating further.
""" ->
stoppropagation(tile::Tile, name::Symbol) =
    StopPropagation(tile, name)

render(tile::StopPropagation) =
    render(tile.tile) <<
        Elem("stop-propagation",
            attributes=[:name=>tile.name])

# Send a signal update to the Julia side
immutable SignalTransport <: Tile
    tile::Tile
    name::Symbol
    signal::Input
end

subscribe(t::Tile, name, s::Input; absorb=true) =
    SignalTransport(t, name, s) |>
       (x -> absorb ? stoppropagation(x, name) : x)

render(sig::SignalTransport) =
    render(sig.tile) <<
        Elem("signal-transport",
            attributes=[:name=>sig.name, :signalId => setup_transport(sig.signal)])

setup_transport(x) =
    error("Looks like there is no trasport set up")

# Utility functions for transports
decodeJSON(sig::Input, val) = val
decodeJSON{T <: String}(sig::Input{T}, ::Nothing) = ""
decodeJSON{T <: String}(sig::Input{T}, val) = string(val)
decodeJSON{T <: Integer}(sig::Input{T}, val) = convert(T, int(val))
decodeJSON{T <: FloatingPoint}(sig::Input{T}, val) = convert(T, float(val))

istruthy(::Nothing) = false
istruthy(b::Bool) = b
istruthy(::None) = false
istruthy(x) = !isempty(x)

decodeJSON(sig::Input{Bool}, val) = istruthy(val)

import Base.Random: UUID, uuid4

const signal_to_id = Dict()
const id_to_signal = Dict()

function makeid(sig::Signal)
    if haskey(signal_to_id, sig)
        # todo ensure connection
        return signal_to_id[sig]
    else
        id = get!(() -> string(uuid4()), signal_to_id, sig)
        id_to_signal[id] = sig
        return id
    end
end

function fromid(id)
    id_to_signal[id]
end


text:
using KomePackage
using Test

@testset "KomePackage.jl group 1" begin
    @test KomePackage.greet() == "Hello Test"
    @test KomePackage.greet("Package") == "Hello Package"
end


text:
# Load ParticleDA
using ParticleDA
#
# # Save some variables for later use
test_dir = joinpath(@__DIR__)
module_src = joinpath(test_dir, "model", "model.jl")
input_file = joinpath(test_dir,"speedy_1.yaml")
# # Instantiate the test environment
using Pkg
Pkg.activate(test_dir)
Pkg.instantiate()

# Include the sample model source code and load it
include(module_src)
using .Model
run_particle_filter(Model.init, input_file, OptimalFilter())

max_stars_repo_name: Team-RADDISH/ParticleDA-UseCases

max_issues_repo_name: Team-RADDISH/ParticleDA-UseCases

max_forks_repo_name: Team-RADDISH/ParticleDA-UseCases


text:
"""
    univariate_kalman_filter(model::StateSpaceModel{Typ}; tol::Typ = Typ(1e-5)) where Typ

Kalman filter with big Kappa initialization, i.e., initializing state variances as 1e6.
"""
function univariate_kalman_filter(model::StateSpaceModel{Typ}; tol::Typ = Typ(1e-5)) where Typ

    time_invariant = model.mode == "time-invariant"

    # Predictive state and its covariance
    a = Matrix{Typ}(undef, model.dim.n+1, model.dim.m)
    P = Array{Typ, 3}(undef, model.dim.m, model.dim.m, model.dim.n+1)
    att = Matrix{Typ}(undef, model.dim.n, model.dim.m)
    Ptt = Array{Typ, 3}(undef, model.dim.m, model.dim.m, model.dim.n)

    # Innovation and its covariance
    v = Vector{Typ}(undef, model.dim.n)
    F = Vector{Typ}(undef, model.dim.n)

    # Kalman gain
    K = Matrix{Typ}(undef, model.dim.m, model.dim.n)

    # Auxiliary structures
    ZP = Vector{Typ}(undef, model.dim.m)
    P_Ztransp_invF = Vector{Typ}(undef, model.dim.m)
    RQR = Matrix{Typ}(undef, model.dim.m, model.dim.m)

    # Steady state initialization
    steadystate = false
    tsteady     = model.dim.n+1

    # Initial state: big Kappa initialization
    fill_a1!(a)
    fill_P1!(P; bigkappa = Typ(1e6))

    mul!(RQR, model.R, LinearAlgebra.BLAS.gemm('N', 'T', Typ(1.0), model.Q, model.R)) # RQR = R Q R'
    # Kalman filter
    for t = 1:model.dim.n
        if t in model.missing_observations
            steadystate  = false
            v[t] = NaN
            F[t] = NaN
            update_att!(att, a, t) # att_t = a_t
            update_Ptt!(Ptt, P, t) # Ptt_t = P_t
            update_a!(a, att, model.T, model.c, t) # a_t+1 = T * att_t + c_t
            update_P!(P, model.T, Ptt, RQR, t) # P_t+1 = T * Ptt_t * T' + RQR'
        elseif steadystate
            update_v!(v, model.y, model.Z, model.d, a, t) # v_t = y_t - Z_t * a_t - d_t
            F[t] = F[t-1]
            repeat_vector_t_plus_1!(K, t-1) # K[:, t]   = K[:, t-1]
            update_att!(att, a, P_Ztransp_invF, v, t) # att_t = a_t + P_t * Z_t * F^-1_t * v_t
            repeat_matrix_t_plus_1!(Ptt, t-1) # Ptt_t = Ptt_t-1
            update_a!(a, att, model.T, model.c, t) # a_t+1 = T * att_t + c_t
            repeat_matrix_t_plus_1!(P, t) # P_t+1 = P_t
        else
            update_v!(v, model.y, model.Z, model.d, a, t) # v_t = y_t - Z_t * a_t - d_t
            update_ZP!(ZP, model.Z, P, t) # ZP = Z[:, :, t] * P[:, :, t]
            update_F!(F, ZP, model.Z, model.H, t) # F_t = Z_t * P_t * Z_t' + H
            update_P_Ztransp_Finv!(P_Ztransp_invF, ZP, F, t) # P_Ztransp_invF   = ZP' * invF(F, t)
            update_K!(K, P_Ztransp_invF, model.T, t) # K_t = T * P_t * Z_t * F^-1_t
            update_att!(att, a, P_Ztransp_invF, v, t) # att_t = a_t + P_t * Z_t * F^-1_t * v_t
            update_Ptt!(Ptt, P, P_Ztransp_invF, ZP, t) # Ptt_t = P_t - P_t * Z_t' * F^-1_t * Z_t * P_t
            update_a!(a, att, model.T, model.c, t) # a_t+1 = T * att_t + c_t
            update_P!(P, model.T, Ptt, RQR, t) # P_t+1 = T * Ptt_t * T' + RQR'
            if time_invariant && check_steady_state(P, t, tol)
                steadystate = true
                tsteady     = t
            end
        end
    end
    # Return the auxiliary filter structre
    return UnivariateKalmanFilter(a, att, v, P, Ptt, F, steadystate, tsteady, K)
end

"""
    univariate_smoother(model::StateSpaceModel{Typ}, kfilter::UnivariateKalmanFilter{Typ}) where Typ

Smoother for state-space model.
"""
function univariate_smoother(model::StateSpaceModel{Typ}, kfilter::UnivariateKalmanFilter{Typ}) where Typ

    # Load dimensions data
    n, p, m, r = size(model)

    # Load system data
    Z, T, R = ztr(model)

    # Load filter data
    a       = kfilter.a
    v       = kfilter.v
    F       = kfilter.F
    P       = kfilter.P
    K       = kfilter.K

    # Smoothed state and its covariance
    alpha = Matrix{Typ}(undef, n, m)
    V     = Array{Typ, 3}(undef, m, m, n)
    L     = Array{Typ, 3}(undef, m, m, n)
    r     = Matrix{Typ}(undef, n, m)
    N     = Array{Typ, 3}(undef, m, m, n)

    # Initialization
    N[:, :, end] = zeros(m, m)
    r[end, :]    = zeros(m, 1)

    @views @inbounds for t = n:-1:2
        if t in model.missing_observations
            r[t-1, :]     = T' * r[t, :]
            N[:, :, t-1]  = T' * N[:, :, t] * T
        else
            Z_transp_invF = Z[1, :, t]/F[t]
            L[:, :, t]    = T - K[:, t] * Z[:, :, t]
            r[t-1, :]     = Z_transp_invF * v[t] + L[:, :, t]' * r[t, :]
            N[:, :, t-1]  = Z_transp_invF * Z[:, :, t] + L[:, :, t]' * N[:, :, t] * L[:, :, t]
        end
        alpha[t, :]   = a[t, :] + P[:, :, t] * r[t-1, :]
        V[:, :, t]    = P[:, :, t] - P[:, :, t] * N[:, :, t-1] * P[:, :, t]
    end

    # Last iteration
    if 1 in model.missing_observations
        r0  = T' * r[1, :]
        N0  = T' * N[:, :, 1] * T
    else
        Z_transp_invF = Z[:, :, 1]' * invertF(F, 1)
        L[:, :, 1]    = T - K[:, 1] * Z[:, :, 1]
        r0            = Z_transp_invF * v[1] + L[:, :, 1]' * r[1, :]
        N0            = Z_transp_invF * Z[:, :, 1] + L[:, :, 1]' * N[:, :, 1] * L[:, :, 1]
    end
    alpha[1, :]   = a[1, :] + P[:, :, 1] * r0
    V[:, :, 1]    = P[:, :, 1] - P[:, :, 1] * N0 * P[:, :, 1]
    # Return the smoothed state structure
    return Smoother(alpha, V)
end

function get_log_likelihood_params(model::StateSpaceModel{T}, filter_type::Type{UnivariateKalmanFilter{T}}) where T <: AbstractFloat

    # Obtain innovation v and its variance F
    kfilter = univariate_kalman_filter(model)

    # Return v and F
    return kfilter.v, kfilter.F
end

function kfas(model::StateSpaceModel{T}, filter_type::Type{UnivariateKalmanFilter{T}}) where T <: AbstractFloat

    # Run filter and smoother
    filtered_state = univariate_kalman_filter(model)
    smoothed_state = univariate_smoother(model, filtered_state)
    v = Matrix{T}(undef, length(filtered_state.v), 1)
    v[:, 1] = filtered_state.v
    F = Array{T}(undef, 1, 1, length(filtered_state.F))
    F[1, 1, :] = filtered_state.F.*ones(1, 1)

    return FilterOutput(filtered_state.a, filtered_state.att, v,
                        filtered_state.P, filtered_state.Ptt, F,
                        filtered_state.steadystate, filtered_state.tsteady),
           SmoothedState(smoothed_state.alpha, smoothed_state.V)
end


text:
function LinearKalmanFilter(arg0::MatrixDecomposer, arg1::LinearProcess, arg2::ProcessEstimate)
    return LinearKalmanFilter((MatrixDecomposer, LinearProcess, ProcessEstimate), arg0, arg1, arg2)
end

function estimation_step(obj::LinearKalmanFilter, arg0::Measurement)
    return jcall(obj, "estimationStep", ProcessEstimate, (Measurement,), arg0)
end

function get_corrected(obj::AbstractKalmanFilter)
    return jcall(obj, "getCorrected", ProcessEstimate, ())
end

function get_predicted(obj::AbstractKalmanFilter)
    return jcall(obj, "getPredicted", ProcessEstimate, ())
end


text:
module TestEstimator
using Test
using JuMP, Ipopt

include("../src/SimplePF.jl")
import .SimplePF
export runtests

include("../src/Estimator.jl")
import .Estimator

include("cases_no_v.jl")
using .Cases

function test_build_DER_set_buses(sys)
    recived = Estimator.build_DER_set_buses(sys)
    println([dg.bus for dg in sys.dgs])
    @test recived == [dg.bus for dg in sys.dgs]

end

function test_build_DER_scenario(sys, name)
    if name == "case3_dist"
        scenario = [
            [1, 1], [1, 2], [1, 3], [2, 1], [2, 2], [2, 3],
        ]
        func = Estimator.build_DER_scenario(sys)
        for (k, K) in enumerate(scenario)
            for (d, expected) in enumerate(K)
                @test func(k, d) == sys.dgs[d].scenario[expected]
            end
        end

    end

end

function test_build_sets(sys)
    Ω = sys.buses
    sub_der = [Estimator.build_DER_set_buses(sys); sys.substation.bus]
    bΩ = filter(bus -> bus ∉ sub_der, Ω)
    L = collect(1:length(sys.m_load))
    D = collect(1:length(sys.dgs))
    if isempty(sys.dgs)
        K = [1]
    else
        K = collect(1:reduce(*, length(der.scenario) for der in sys.dgs))
    end

    S = collect(1:length(sys.m_new_dg))

    (Ωr, bΩr, Lr, Kr, Dr, Sr) = Estimator.build_sets(sys)
    @test Ω == Ωr
    @test bΩ == bΩr
    @test L == Lr
    @test K == Kr
    @test D == Dr
    @test S == Sr
end

function test_add_variables(sot, sys)
    (Ω, bΩ, L, K, D, S) = Estimator.build_sets(sys)
    n_L = length(L)
    n_D = length(D)
    n_K = length(K)
    n_S = length(S)
    n_bus = length(Ω)

    n_volts = 2 * n_bus * n_L * n_K * n_S
    n_p_der = n_D * n_L * n_K * n_S
    n_q_der = n_D * n_L * n_K * n_S
    n_p_hc = 1

    n_var = n_volts + n_p_der + n_q_der + n_p_hc

    sot = Estimator.add_variables(sot, sys)
    @testset "add_variables" begin
        @test sot isa Model
        V = sot[:V]
        pᴰᴱᴿ = sot[:pᴰᴱᴿ]
        qᴰᴱᴿ = sot[:qᴰᴱᴿ]
        pᴴᶜ = sot[:pᴴᶜ]
        @test num_variables(sot) == n_var

        @testset "V[$b, $l, $k, $s]" for b = Ω, l = L, k = K, s = S
            @testset "Start values" begin
                @test start_value(V[:Re, b, l, k, s]) == 1.0
                @test start_value(V[:Im, b, l, k, s]) == 0.0
            end

            @testset "Lower bound" begin
                @test lower_bound(V[:Re, b, l, k, s]) == -sys.VH
                @test lower_bound(V[:Im, b, l, k, s]) == -sys.VH
            end
            @testset "Upper bound" begin
                @test upper_bound(V[:Re, b, l, k, s]) == sys.VH
                @test upper_bound(V[:Im, b, l, k, s]) == sys.VH
            end
        end

        @testset "pᴰᴱᴿ[$d, $l, $k, $s]" for d = D, l = L, k = K, s = S
            @testset "Start values" begin
                @test start_value(pᴰᴱᴿ[d, l, k, s]) == 0.0
            end
            @testset "Bounds" begin
                @test lower_bound(pᴰᴱᴿ[d, l, k, s]) == sys.dgs[d].P_limit[1]
                @test upper_bound(pᴰᴱᴿ[d, l, k, s]) == sys.dgs[d].P_limit[2]
            end
        end

        @testset "qᴰᴱᴿ[$d, $l, $k, $s]" for d = D, l = L, k = K, s = S
            @testset "Start values" begin
                @test start_value(qᴰᴱᴿ[d, l, k, s]) == 0.0
            end
            @testset "Bounds" begin
                @test lower_bound(qᴰᴱᴿ[d, l, k, s]) == sys.dgs[d].Q_limit[1]
                @test upper_bound(qᴰᴱᴿ[d, l, k, s]) == sys.dgs[d].Q_limit[2]
            end
        end
        @testset "pᴴᶜ" begin
            @testset "Start values" begin
                @test start_value(pᴴᶜ) == 0.0
            end
            @testset "Bounds" begin
                @test lower_bound(pᴴᶜ) == 0.0
                @test has_upper_bound(pᴴᶜ) == false
            end
        end

    end

    return sot
end

function test_add_voltage_constraints(sot, sys)
    (Ω, bΩ, L, K, D, S) = Estimator.build_sets(sys)

    sot = Estimator.add_voltage_constraints(sot, sys)
    @test sot isa Model

    V = sot[:V]

    @test sot[:V²] isa JuMP.Containers.DenseAxisArray{JuMP.QuadExpr,4}
    V² = sot[:V²]

    @test sot[:V_module] isa JuMP.Containers.DenseAxisArray{JuMP.NonlinearExpression,4}
    V_module = sot[:V_module]

    voltage_constraint = sot[:voltage_constraint] #How to test it?

    for b = Ω, l = L, k = K, s = S
        @test V²[b, l, k, s] == V[:Re, b, l, k, s]^2 + V[:Im, b, l, k, s]^2
        obj = constraint_object(voltage_constraint[b, l, k, s])
        @test isequal_canonical(obj.func, V²[b, l, k, s])
        @test obj.set == MOI.Interval(sys.VL^2, sys.VH^2)

    end
    return sot
end

function test_add_I_V_relationship(sot, sys)
    (Ω, bΩ, L, K, D, S) = Estimator.build_sets(sys)
    G = real(sys.Y)
    B = imag(sys.Y)
    V = sot[:V]

    sot = Estimator.add_I_V_relationship(sot, sys)
    @test sot isa Model
    @test sot[:I] isa JuMP.Containers.DenseAxisArray{JuMP.AffExpr,5}
    I = sot[:I]

    for i = Ω, l = L, k = K, s = S
        @test I[:Re, i, l, k, s] == sum(
            G[i, j] * V[:Re, j, l, k, s] - B[i, j] * V[:Im, j, l, k, s]
            for j = Ω)

        @test I[:Im, i, l, k, s] == sum(
            B[i, j] * V[:Re, j, l, k, s] + G[i, j] * V[:Im, j, l, k, s]
            for j = Ω)
    end

    return sot
end

function test_add_S_VI_relationship(sot, sys)
    (Ω, bΩ, L, K, D, S) = Estimator.build_sets(sys)

    sot = Estimator.add_S_VI_relationship(sot, sys)
    V = sot[:V]
    I = sot[:I]

    @test sot[:P] isa JuMP.Containers.DenseAxisArray{JuMP.QuadExpr,4}
    @test sot[:Q] isa JuMP.Containers.DenseAxisArray{JuMP.QuadExpr,4}
    P = sot[:P]
    Q = sot[:Q]

    for b = Ω, l = L, k = K, s = S
        @test P[b, l, k, s] == SimplePF.mc_re(V[:Re, b, l, k, s], V[:Im, b, l, k, s], I[:Re, b, l, k, s], -I[:Im, b, l, k, s])

        @test Q[b, l, k, s] == SimplePF.mc_im(V[:Re, b, l, k, s], V[:Im, b, l, k, s], I[:Re, b, l, k, s], -I[:Im, b, l, k, s])
    end

    return sot
end

function test_add_substation_constraint(sot, sys)
    (Ω, bΩ, L, K, D, S) = Estimator.build_sets(sys)
    sot = Estimator.add_substation_constraint(sot, sys)
    @test sot isa Model

    V = sot[:V]
    I = sot[:I]
    P = sot[:P]
    Q = sot[:Q]
    i = sys.substation.bus
    Plimit = float(sys.substation.P_limit)
    Qlimit = float(sys.substation.Q_limit)

    for l = L, k = K, s = S
        @test !is_fixed(V[:Re, i, l, k, s])
        @test is_fixed(V[:Im, i, l, k, s])
        @test start_value(V[:Re, i, l, k, s]) == sys.substation.voltage
        @test fix_value(V[:Im, i, l, k, s]) == 0

        obj = constraint_object(sot[:sub_plimit][l, k, s])
        @test obj.set == MOI.Interval(0.0, Plimit)
        @test isequal_canonical(obj.func, P[i, l, k, s])

        obj = constraint_object(sot[:sub_qlimit][l, k, s])
        @test obj.set == MOI.Interval(0.0, Qlimit)
        @test isequal_canonical(obj.func, Q[i, l, k, s])

        obj = constraint_object(sot[:sub_current][l, k, s])
        @test isequal_canonical(obj.func, I[:Re, i, l, k, s]^2 + I[:Im, i, l, k, s]^2)
    end

    return sot
end

function test_add_line_constraint(sot, sys)

    G = real(sys.Y)
    B = imag(sys.Y)

    (Ω, bΩ, L, K, D, S) = Estimator.build_sets(sys)
    F = sys.amp
    Ωᴮ = keys(F)
    sot = Estimator.add_line_constraint(sot, sys)

    @test sot isa Model

    f = sot[:f]
    V = sot[:V]
    for (i, j) ∈ Ωᴮ, l = L, k = K, s = S
        vij_re = V[:Re, i, l, k, s] - V[:Re, j, l, k, s]
        vij_im = V[:Im, i, l, k, s] - V[:Im, j, l, k, s]
        @test f[:Re, (i, j), l, k, s] == vij_re * G[i, j] - vij_im * B[i, j]
        @test f[:Im, (i, j), l, k, s] == vij_im * G[i, j] + vij_re * B[i, j]

        obj = constraint_object(sot[:line_limit][(i, j), l, k, s])
        @test isequal_canonical(obj.func, f[:Re, (i, j), l, k, s]^2 + f[:Im, (i, j), l, k, s]^2)
        @test obj.set == MOI.LessThan(F[(i, j)]^2)
    end

    return sot

end

function test_add_power_injection_definition(sot, sys)
    (Ω, bΩ, L, K, D, S) = Estimator.build_sets(sys)
    B = Estimator.build_DER_set_buses(sys)
    sot = Estimator.add_power_injection_definition(sot, sys)
    PL = sys.PL
    QL = sys.QL
    subbus = sys.substation.bus

    P = sot[:P]
    Q = sot[:Q]
    pᴴᶜ = sot[:pᴴᶜ]
    pᴰᴱᴿ = sot[:pᴰᴱᴿ]
    qᴰᴱᴿ = sot[:qᴰᴱᴿ]

    @testset "add_power_injection_definition" begin
        p = sot[:p]
        q = sot[:q]
        @testset "Without DER" for b = bΩ, l = L, k = K, s = S
            μᴸ = sys.m_load[l]
            μᴴᶜ = sys.m_new_dg[s]
            @testset "No definition for substation bus" begin
                @test_throws KeyError p[subbus, l, k, s]
                @test_throws KeyError q[subbus, l, k, s]
            end
            @testset "Active" begin
                obj = constraint_object(p[b, l, k, s])
                @test isequal_canonical(obj.func, P[b, l, k, s] - μᴴᶜ * pᴴᶜ)
                @test obj.set == MOI.EqualTo(-μᴸ * PL[b])
            end
            @testset "Reactive" begin
                obj = constraint_object(q[b, l, k, s])
                @test isequal_canonical(obj.func, Q[b, l, k, s])
                @test obj.set == MOI.EqualTo(-μᴸ * QL[b])
            end
        end

        p_wder = sot[:p_wder]
        q_wder = sot[:q_wder]
        μᴰᴱᴿ = Estimator.build_DER_scenario(sys)
        @testset "With DER" for d = D, l = L, k = K, s = S
            b = B[d]
            Pᴰᴱᴿ = (1.0 - sys.dgs[d].alpha) * sys.dgs[d].S_limit
            μᴸ = sys.m_load[l]
            @testset "Active" begin
                obj = constraint_object(p_wder[d, l, k, s])
                @test isequal_canonical(obj.func, P[b, l, k, s] - pᴰᴱᴿ[d, l, k, s])
                @test obj.set == MOI.EqualTo(-μᴸ * PL[b] + Pᴰᴱᴿ * μᴰᴱᴿ(k, d))
            end
            @testset "Reactive" begin
                obj = constraint_object(q_wder[d, l, k, s])
                @test isequal_canonical(obj.func, Q[b, l, k, s] - qᴰᴱᴿ[d, l, k, s])
                @test obj.set == MOI.EqualTo(-μᴸ * QL[b])
            end
        end
    end

    return sot
end

function test_add_ders_limits(sot, sys)
    #P and Q limits was alredy tests in test_add_variables
    (Ω, bΩ, L, K, D, S) = Estimator.build_sets(sys)
    sot = Estimator.add_ders_limits(sot, sys)
    pᴰᴱᴿ = sot[:pᴰᴱᴿ]
    qᴰᴱᴿ = sot[:qᴰᴱᴿ]

    disco_der_limit = sot[:disco_der_limit]
    der_limit = sot[:der_limit]
    der_energy_limit = sot[:energy_limit]
    μᴰᴱᴿ = Estimator.build_DER_scenario(sys)
    for d = D, l = L, k = K, s = S
        obj = constraint_object(disco_der_limit[d, l, k, s])
        @test isequal_canonical(obj.func, pᴰᴱᴿ[d, l, k, s]^2 + qᴰᴱᴿ[d, l, k, s]^2)
        @test obj.set == MOI.LessThan((sys.dgs[d].alpha * sys.dgs[d].S_limit)^2)

        Pᴰᴱᴿ = (1.0 - sys.dgs[d].alpha) * sys.dgs[d].S_limit
        obj = constraint_object(der_limit[d, l, k, s])
        canonical = pᴰᴱᴿ[d, l, k, s]^2 + qᴰᴱᴿ[d, l, k, s]^2 + 2 * Pᴰᴱᴿ * μᴰᴱᴿ(k, d) * pᴰᴱᴿ[d, l, k, s]
        @test isequal_canonical(obj.func, canonical)
        left = (sys.dgs[d].S_limit)^2 - (Pᴰᴱᴿ * μᴰᴱᴿ(k, d))^2
        @test obj.set == MOI.LessThan(left)

        Eᴰᴱᴿ = sys.dgs[d].energy
        βᴰᴱᴿ = sys.dgs[d].beta
        obj = constraint_object(der_energy_limit[d, l, k, s])
        @test isequal_canonical(obj.func, pᴰᴱᴿ[d, l, k, s] * sys.time_curr)
        @test obj.set == MOI.Interval(-βᴰᴱᴿ * Eᴰᴱᴿ, βᴰᴱᴿ * Eᴰᴱᴿ)
    end

    return sot
end

function test_nl_pf(sot, sys, case)

    other_sys, _ = case()
    model = Estimator.nl_pf(Model(), other_sys)
    @testset "nl_pf" begin
        #Find a better way to test it
        @test sprint(print, sot) == sprint(print, model)
    end
    return sot
end

function test_model_should_be_feasible!(sot, sys, name)
    @testset "solve $(name)" begin
        set_optimizer(sot, Ipopt.Optimizer)
        set_silent(sot)
        optimize!(sot)

        @testset "Optimization" begin
            @testset "Solved: $(termination_status(sot))" begin
                @test termination_status(sot) == LOCALLY_SOLVED
            end
            @testset "Feasible: $(primal_status(sot))" begin
                @test primal_status(sot) == FEASIBLE_POINT
            end
        end
    end
    return sot
end

function test_objective(sot, sys)
    (Ω, bΩ, L, K, D, S) = Estimator.build_sets(sys)
    sot = Estimator.add_objective(sot, sys)
    Nᴮ = length(bΩ)
    pᴴᶜ = sot[:pᴴᶜ]
    @test MAX_SENSE == @inferred objective_sense(sot)
    @test objective_function_type(sot) == JuMP.AffExpr
    @test isequal_canonical(objective_function(sot), Nᴮ * pᴴᶜ)
    return sot
end

function runtests()

    @testset "Estimator" begin
        for case in [case3_dist, case33, case3_dist_no_dgs]
            sys, name = case()
            @testset "$name" begin
                test_build_DER_scenario(sys, name)
                test_build_sets(sys)
                sot = Model()
                sot = test_add_variables(sot, sys)
                sot = test_add_voltage_constraints(sot, sys)
                sot = test_add_I_V_relationship(sot, sys)
                sot = test_add_S_VI_relationship(sot, sys)
                sot = test_add_substation_constraint(sot, sys)
                sot = test_add_line_constraint(sot, sys)
                sot = test_add_power_injection_definition(sot, sys)
                sot = test_add_ders_limits(sot, sys)
                sot = test_nl_pf(sot, sys, case)
                sot = test_model_should_be_feasible!(sot, sys, name)
                sot = test_objective(sot, sys)
            end
        end
    end
end

end

TestEstimator.runtests()


text:
import Base: map, ==, hash, print

const Occupancy = Int; const EMPTY=0; const BLACK=1; const WHITE=2; const DISABLED=3
@enum Pieces L=5 Z=10 O=15 T=20
const Square = Tuple{Int,Int,Int} # one voxel
const Squares = Set{Square} # a set of voxels
const Board = Array{Occupancy, 3}
const Spot = Tuple{Int, Int} # just x/y, ignoring height
# struct Shape
#     blocks::Squares # shape locations (assuming 0,0,0 is the top-left-bottom of the cube) - these squares need to be empty to put shape
#     overhang::Squares # squares that need to be non-empty to put
# end
# ==(shape1::Shape, shape2::Shape) = shape1.blocks == shape2.blocks && shape1.overhang == shape2.overhang
# hash(shape::Shape) = hash(shape.blocks) + hash(shape.overhang)
const Shape = NamedTuple{(:blocks, :overhang), Tuple{Squares, Squares}}
const Move = NamedTuple{(:piece, :shape, :spot, :colour), Tuple{Pieces, Shape, Spot, Occupancy}}
EMPTY_MOVE = Move((O, Shape((Squares(), Squares())), (0,0), DISABLED))

struct StaticState

end

const PieceShapes = Dict(
    L => Shape((Set([(0,0,0), (1,0,0), (2,0,0), (2,1,0)]), Set{Squares}())), # L flat
    O => Shape((Set([(0,0,0), (0,1,0), (1,0,0), (1,1,0)]), Set{Squares}())), # O flat
    T => Shape((Set([(0,0,0), (1,0,0), (2,0,0), (1,1,0)]), Set{Squares}())), # T flat
    Z => Shape((Set([(0,0,0), (1,0,0), (1,1,0), (2,1,0)]), Set{Squares}())), # Z flat
    # Shape((Set([(0,0,0), (0,1,0), (0,0,1), (0,1,1)]), Set{Squares}())), # O horizontal
    # Shape((Set([(0,0,0), (1,0,0), (0,0,1), (1,0,1)]), Set{Squares}())), # O vertial
)

map(fun, dict::Dict) =  Dict(pair[1] => fun(pair[2]) for pair in dict)
map(fun, set::Set) =  Set([fun(x) for x in set])

mutable struct State
    board::Board
    white_pieces::Vector{Pieces}
    black_pieces::Vector{Pieces}
    whites_turn::Bool
    squares_to_connect::Vector{Tuple{Squares, Squares}}
    last_move::Move
end

function State()
    board = zeros(Occupancy,8,8,3)

    max_x, max_y, max_z = size(board)
    min_x, min_y, min_z = (1,1,1)
    left_side = [(x,y,z) for x = min_x for y = min_y:max_y for z = min_z:max_z]
    right_side = [(x,y,z) for x = max_x for y = min_y:max_y for z = min_z:max_z]
    top_side = [(x,y,z) for x = min_x:max_x for y = min_y for z = min_z:max_z]
    bottom_side = [(x,y,z) for x = min_x:max_x for y = max_y for z = min_z:max_z]
    squares_to_connect = [
        (Squares(left_side), Squares(right_side)),
        (Squares(top_side), Squares(bottom_side)),
        ]

    State(
        board,
        [L, L, L, L, Z, Z, Z, Z, T, T, T, T, O, O],
        [L, L, L, L, Z, Z, Z, Z, T, T, T, T, O, O],
        true,
        squares_to_connect,
        EMPTY_MOVE
    ) # default board
end

const piece_state = State(zeros(Occupancy,3,3,3), [], [], true, [], EMPTY_MOVE)
# const squares_to_connect = Vector{Tuple{Squares, Squares}}() # defines win condition

function HasSomeoneWon(state::State, last_move::Move)
    # if one of the BridgeComponents of a colour contains squares on opposing edges
    if IsBridgeBuilt(state, last_move)
        return (true, last_move.colour)
    end
    # game over, but nobody won
    if isempty(state.white_pieces) && isempty(state.black_pieces)
        return (true, DISABLED)
    end
    # nope, play on
    return (false, DISABLED)
end

function top_square(state::State, spot::Spot)
    spot_column = state.board[spot..., :]
    if spot_column[3] in [BLACK, WHITE]
        return (3, spot_column[3])
    elseif spot_column[2] in [BLACK, WHITE]
        return (2, spot_column[2])
    elseif spot_column[1] in [BLACK, WHITE]
        return (1, spot_column[1])
    else
        return (0, EMPTY)
    end
end

function BridgeNeighbours(board::Board, square::Square)
    # Neighbours of a square that can participate in a bridge
    # square needs to be on top, connections are North, West, East, South on the same level or one
    # level up/down
    #TODO two levels when all colours are the same :(
    #TODO column major bug?
    colour = board[square...]
    max_x, max_y, max_z = size(board)
    min_x, min_y, min_z = (1,1,1)
    possible_neighbour_spots = [
        (square[1] + 1, square[2], top_squares(square[1] + 1, square[2])),
        (square[1] - 1, square[2], top_squares(square[1] - 1, square[2])),
        (square[1], square[2] + 1, top_squares(square[1], square[2] + 1)),
        (square[1], square[2] - 1, top_squares(square[1], square[2] - 1)),
    ]
    possible_neighbours = [(x,y,z) for (x,y,z,col) in possible_neighbour_spots if col=colour]
    # TODO: special case with 1,3 and 3,1
    if square[3] == 3 && new_square[3] == 1 && board[new_square[1], new_square[2], 2] != colour ||
       square[3] == 1 && new_square[3] == 3 && board[square[1], square[2], 2] != colour
    # possible_neighbours = [(x,y,z) for (x,y,z)=possible_neighbours if min_x <= x <= max_x && min_y <= y <= max_y && min_z <= z <= max_z && on_top(board, x, y, z) && board[x,y,z] == colour]
    return Set(possible_neighbours)
end

on_top(board, x, y, z) = z == size(board)[3] || board[x,y,z+1] == EMPTY

function IsBridgeBuilt(state::State, last_move::Move)
    # has a bridge been formed?
    bollocks = squares_in_spot(last_move.shape.blocks, last_move.spot)
    bridge_to = Squares()
    for pairs in state.squares_to_connect
        # if we connected to a side this move, see if we can make it to the other side
        if !isempty(intersect(pairs[1], bollocks))
            bridge_to = union(bridge_to, pairs[2])
        end
        if !isempty(intersect(pairs[2], bollocks))
            bridge_to = union(bridge_to, pairs[1])
        end
    end

    if isempty(bridge_to) return (false) end

    starting_block = [(x,y,z) for (x,y,z)=bollocks if on_top(state.board, x, y, z)][1] # something on top - should we choose one on the side?
    @assert state.board[starting_block...] == last_move.colour
    reachable_blocks = connected_component(state.board, starting_block)
    bridge_built = !isempty(intersect(reachable_blocks, bridge_to))
    return bridge_built
end

function connected_component(board::Board, starting_block::Square)
    @assert on_top(board, starting_block...)
    component = Squares([starting_block])
    new_squares = Squares([starting_block])
    while true
        neighbouring_squares = [BridgeNeighbours(board, square) for square=new_squares]
        neighbouring_squares = union(neighbouring_squares...)
        new_squares = setdiff(neighbouring_squares, component)
        component = union(component, neighbouring_squares)
        if isempty(new_squares) break end
    end
    return component
end

function empty_and_occupied_sqares(board::Array{Occupancy, 3})
    empty_squares = Squares()
    occupied_squares = Squares()

    xlim, ylim, zlim = size(board)
    for z in 1:zlim
        for y in 1:ylim
            for x in 1:xlim
                if board[x,y,z] == EMPTY
                    push!(empty_squares, (x,y,z))
                elseif board[x,y,z] != DISABLED
                    push!(occupied_squares, (x,y,z))
                end
            end
        end
    end
    return empty_squares, occupied_squares
end

function can_put(board, shape::Shape, x_pos, y_pos)
    # TODO: make sure it fits... 2 disabled rows/columns?
    has_space = all([board[x_pos + x, y_pos + y , z + 1] == EMPTY for (x,y,z) = shape.blocks])
    no_overhang = all([board[x_pos + x, y_pos + y , z + 1] != EMPTY for (x,y,z) = shape.overhang])
    return has_space && no_overhang
end

function rotations(shape::Shape)
    all_shapes = Set{Shape}()
    for z=1:4
        for y=1:4
            for x=1:4
                push!(all_shapes, shape)
                shape = apply_rotate(x_rotate, shape)
            end
            shape = apply_rotate(y_rotate, shape)
        end
        shape = apply_rotate(z_rotate, shape)
    end
    return all_shapes
end

function apply_rotate(fun, shape)
    blocks = map(fun, shape.blocks)
    blocks = adjust_shape_vec(blocks)
    overhang = calculate_overhang(blocks)
    return Shape((blocks, overhang))
end

const ROT_MAX = 2
x_rotate(square::Square) = (square[1], square[3], ROT_MAX-square[2])
y_rotate(square::Square) = (ROT_MAX-square[3], square[2], square[1])
z_rotate(square::Square) = (square[2], ROT_MAX-square[1], square[3])

function adjust_shape_vec(squares::Squares)
    # make sure lowest coordinate is 0 (x,y and z)
    min_x = min([x for (x,y,z) = squares]...)
    min_y = min([y for (x,y,z) = squares]...)
    min_z = min([z for (x,y,z) = squares]...)
    new_blocks = map(square -> (square[1] - min_x, square[2] - min_y, square[3] - min_z), squares)
    return new_blocks
end

squares_in_spot(squares::Squares, spot::Spot) = map(square -> (square[1] + spot[1], square[2] + spot[2], square[3] + 1), squares)

# suares under this one
underlings(square::Square) = Set([(square[1], square[2], square[3]-i) for i = 1:square[3]])

function calculate_overhang(blocks::Squares)
    overhangs = union([underlings(square) for square = blocks]...)
    return setdiff(overhangs, blocks)
end

const PieceRotations = map(rotations, PieceShapes)

#TODO speedup: all_possible_moves in advance, and remove impossible ones as we go + can_put

function valid_moves(state::State)
    colour = state.whites_turn ? WHITE : BLACK
    possible_moves = Vector{Move}()
    # all_moves = Vector{Shape}()
    empty_squares, occupied_squares = empty_and_occupied_sqares(state.board)
    available_pieces = colour == WHITE ? unique(state.white_pieces) : unique(state.black_pieces)
    #uses PieceRotations and the board
    for x in 1:size(state.board)[1]
        for y in 1:size(state.board)[1]
            for piece in available_pieces
                for rotation in PieceRotations[piece]
                    # push!(all_moves, rotation)
                    if issubset(squares_in_spot(rotation.blocks, (x,y)), empty_squares) && issubset(squares_in_spot(rotation.overhang, (x,y)), occupied_squares)
                    # if can_put(state.board, rotation, x, y)
                        push!(possible_moves, Move((piece, rotation, (x, y), colour)))
                    end
                end
            end
        end
    end
    return possible_moves
end

function make_move!(state::State, move::Move)
    for square in move.shape.blocks
        @assert state.board[square[1] + move.spot[1], square[2] + move.spot[2], square[3] + 1] == EMPTY
        state.board[square[1] + move.spot[1], square[2] + move.spot[2], square[3] + 1] = move.colour
    end
    if state.whites_turn # TODO: dict of defaultdict
        deleteat!(state.white_pieces, findfirst(x->x==move.piece, state.white_pieces))
    else
        deleteat!(state.black_pieces, findfirst(x->x==move.piece, state.black_pieces))
    end
    state.whites_turn = !state.whites_turn
    state.last_move = move
end

function valid_moves_spot(state::State, spot::Spot, shape::Shape)
    #piece will go in a 3x3x3 cube centered at spot. With the heights known (0,1,2,X) this can
    #be precalculated for every board combination and every piece rotation, and pulled from a cache

    #to make it unique, piece has to occupy a spot on both the top and left side (cache might only contain these)
end

function print(shape::Shape)

end

const blocks = ["■□▣◾◼◻◽"]

const board_show =
"..1221.
...1..." # colors: black - black top and under, color indicates if bottom is different https://stackoverflow.com/questions/27929766/how-to-print-text-with-a-specified-rgb-value-in-julia

# Random policy
function next_move_random(state::State, possible_moves::Vector{Move})
    selected_move = rand(possible_moves)
    return selected_move
end

function play_game(white_policy, black_policy)
    state = State()
    game_over, winner = false, DISABLED
    while true
        print('.')
        possible_moves = valid_moves(state)
        if isempty(possible_moves)
            print("No valid moves")
            break
        end
        next_move = state.whites_turn ? white_policy(state, possible_moves) : black_policy(state, possible_moves)
        make_move!(state, next_move)
        game_over, winner = HasSomeoneWon(state, next_move)
        game_over && break
    end
    print("Game over, Winner: is $winner")
    return state, winner
end

state, winner = play_game(next_move_random, next_move_random)

max_stars_repo_name: robertsugar/bridget

max_issues_repo_name: robertsugar/bridget

max_forks_repo_name: robertsugar/bridget


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Drop <: Edit
        col::Symbol
        val::Any
        operation::String
    end

Remove information from the dataframe - either an entire column or rows containing specified values.

# Arguments
- `col::Symbol`: name of column containing data to remove
- `val::Any`: value to drop
- `operation::String`: how to determine what to drop
"""
mutable struct Drop <: Edit
    "name of column containing data to remove"
    col::Symbol
    "value to drop"
    val::Any
    "how to determine what to drop"
    operation::String
end

function Drop(; col, val, operation, )
    Drop(col, val, operation, )
end

"""Get [`Drop`](@ref) `col`."""
get_col(value::Drop) = value.col
"""Get [`Drop`](@ref) `val`."""
get_val(value::Drop) = value.val
"""Get [`Drop`](@ref) `operation`."""
get_operation(value::Drop) = value.operation

"""Set [`Drop`](@ref) `col`."""
set_col!(value::Drop, val) = value.col = val
"""Set [`Drop`](@ref) `val`."""
set_val!(value::Drop, val) = value.val = val
"""Set [`Drop`](@ref) `operation`."""
set_operation!(value::Drop, val) = value.operation = val

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
using PyPlot
using PyCall
using Base64
anim = pyimport( "matplotlib.animation")

function display_animation(name, make_axis, init, animate, frames=100, interval=20)

    #Construct Figure and do basic setup
    (fig, axis) = make_axis()

    # Create the animation object by calling the Python function FuncAnimaton
    myanim = anim.FuncAnimation(fig, (i) -> animate(axis,i), init_func= ()-> init(axis), frames=frames, interval=interval)

    # Convert it to an MP4 movie file and saved on disk in this format.
    #myanim[:save]("3Lines.mp4", bitrate=-1, extra_args=["-vcodec", "libx264", "-pix_fmt", "yuv420p"])
    myanim.save(name*".mp4", bitrate=-1, extra_args=["-vcodec", "libx264", "-pix_fmt", "yuv420p"])

    # Function for creating an embedded video given a filename
    function html_video(filename)
        open(filename) do f
            base64_video = Base64.base64encode(f)
            """<video controls src="data:video/x-m4v;base64,$base64_video">"""
        end
    end

    # Display the movie in a Julia cell as follows. Note it has animation controls for the user.
    #display("text/html", html_video("3Lines.mp4"))
    display("text/html", html_video(name*".mp4"))

    PyPlot.close(fig)
end

function test_animation()
    function make_axis()
        fig = figure("MyFigure",figsize=(5,5))
        ax = PyPlot.axes(xlim = (0,10),ylim=(0,10))
        (fig, ax)
    end

    # Define the init function, which draws the first frame (empty, in this case)
    function init(ax)
        xline1 = ax.plot([],[],"r-")[1]
        xline2 = ax.plot([],[],"g-")[1]
        xline3 = ax.plot([],[],"b-")[1]
        xline1.set_data([],[])
        xline2.set_data([],[])
        xline3.set_data([],[])
        return (xline1,xline2,xline3,Union{})  # Union{} is the new word for None
    end

    # Animate draws the i-th frame, where i starts at i=0 as in Python.
    function animate(ax,i)
        xline1 = ax.plot([],[],"r-")[1]
        xline2 = ax.plot([],[],"g-")[1]
        xline3 = ax.plot([],[],"b-")[1]
        x = (0:i)/10.0
        xline1.set_data(x,x)
        xline2.set_data(1 .+ x,x)
        xline3.set_data(2 .+ x,x)
        return (xline1,xline2,xline3,Union{})
    end

    display_animation("test1",make_axis,init,animate)
end

max_stars_repo_name: jamesonquinn/Gen

max_issues_repo_name: jamesonquinn/Gen

max_forks_repo_name: jamesonquinn/Gen


text:
module BirdPathTrigger

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Bird Path" => Ahorn.EntityPlacement(
        Maple.BirdPathTrigger,
        "rectangle"
    )
)

end


text:
# Clause set operation

struct RESDB
  cdb
  ldb
  pgr
end

#==
 CDB(clauses) -> LDB(Literals), PGR(Graph)
==#

include("common.jl")
include("reso.jl")

#==
 primitives for CDB
==#

function getvars(cdb, cid)
 cdb[cid][1]
end

function getvars(db::RESDB, cid)
 db.cdb[cid][1]
end

function getcls(cdb, cid)
 if cid == 0; return [[],EmptyClause] end
 cdb[cid]
end

function getcls(rdb::RESDB, cid)
 if cid == 0; return [[],EmptyClause] end
 rdb.cdb[cid]
end

function getlit(cdb, clp)
 cid,lid=clp
 return cdb[cid][2][lid].args[1], cdb[cid][2][lid].args[2]
end

function getlit(db::RESDB, clp)
 cid,lid=clp
 return db.cdb[cid][2][lid].args[1], db.cdb[cid][2][lid].args[2]
end

function litcounts(cdb)
 map(x->length(x[2]),cdb)
end

function litcounts(db::RESDB)
 map(x->length(x[2]),db.cdb)
end

#==
 maker for CDB,LDB,PGR
==#
makecent(cls) = (cls.args[1].args, cls.args[2].args[1].args)

function makecdb(cdb,cls)
  vars, cl = makecent(cls)
  push!(cdb, (vars, cl))
  return cdb
end

function addcdb(cdb, clss)
 for cls in clss
  makecdb(cdb,cls)
 end
 cdb
end

function makeldb(ldb, cid, cls)
@show cls
 for lid in 1:length(cls[2])
  lit=cls[2][lid]
@show lid
@show lit
  ldb[(cid,lid)]=lit
 end
 return ldb
end

function makepg(ldb)
 pgr=Dict()

 # ADHOC
 psyms = map(x->x.args[2].args[1],values(ldb))
 psyms = union(psyms,psyms)
 map(p->pgr[p]=[[],[]],psyms)

 for lid in keys(ldb)
  sign = ldb[lid].args[1]
  psym = ldb[lid].args[2].args[1]

  if sign == :+
    push!(pgr[psym][1],lid)
  else # sign == :-
    push!(pgr[psym][2],lid)
  end

 end
 return pgr
end

function makedb(clauses)
 cdb = []
 ldb = Dict()
 for cid in 1:length(clauses)
  cls = clauses[cid]
  cdb = makecdb(cdb, cls)
  ldb = makeldb(ldb, cid, cls)
 end
 pgr = makepg(ldb)
 return RESDB(cdb, ldb, pgr)
end

##
newcd(cdb) = length(cdb) + 1

function putcls(cls,cdb)::Int
 cid = 1
 for ci in cdb
  if equalclause(ci, cls)
   return cid
  end
  cid += 1
 end
 push!(cdb, cls)
 return cid +1
end

function putcdb(cdb,cls)
 push!(cdb, cls)
end

function putcdb(cls, cdb,ldb,pgr)
 cid = newcid(cdb)
 cdb = putcdb(cdb,cls)
 ldb = putldb(ldb,cid, cls)
 pgr = putpgr(pgr,cls)

 return cdb, ldb, pgr
end

function replacecdb(clsb, clsa, db)
@nyi :replacecdb
end

findunit(cdb) = find(x->x==1,litcounts(cdb))

#==
printing
==#

function printvars(vars)
 print("[")
 for i in 1:length(vars)
  if i < length(vars); print("$(vars[i]),") end
  if i == length(vars); print("$(vars[i])") end
 end
 print("]")
end

function printcdb(cdb)
 for i in 1:length(cdb)
  print("$i ")
  printvars(cdb[i][1])
  println(cdb[i][2])
 end
end

function printplist(plist)
 for i in 1:length(plist)
  if i < length(plist); print("$(plist[i]),") end
  if i == length(plist); print("$(plist[i])") end
 end
end

function printldb(ldb)
 for pr in keys(ldb)
  println("$pr => $(ldb[pr])")
 end
end

function printpgr(pgr)
 for p in keys(pgr)
  print("$p: ")
  print("\n +")
  printplist(pgr[p][1])
  print("\n -")
  printplist(pgr[p][2])
  println()
 end
end

## introduce a type for DB(cdb,ldb,pgr)

function printdb(resdb)
println("Clauses")
  printcdb(resdb.cdb)
println("\nLiterals")
  printldb(resdb.ldb)
println("\nGraph")
  printpgr(resdb.pgr)
end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
# Define utility for types Complex{D} where D<:Dual

# -------------------- #
# Complex{D} utilities #
# -------------------- #
const RealCplxDual = Union{D,Complex{D}} where D<:Dual
const RealCplxDualOrder2 = Union{D,Complex{D}} where D<:Dual{T,V,N} where {T,V<:Dual,N}

@inline value(z::Complex{D}) where D <: Dual = complex(value(real(z)),value(imag(z)))
@inline partials(z::Complex{D}, i::Int) where D <: Dual = complex(partials(real(z),i),partials(imag(z),i))
@inline partials(z::Complex{D}, i::Int, j::Int) where D <: Dual = complex(partials(real(z),i,j),partials(imag(z),i,j))

@inline npartials(z::Complex{D}) where D <: Dual = npartials(D)
@inline npartials(::Type{Complex{D}}) where D <: Dual = npartials(D)
@inline tagtype(z::Complex{D}) where D <: Dual = tagtype(D)
@inline tagtype(::Type{Complex{D}}) where D <: Dual = tagtype(D)
@inline valtype(z::Complex{D}) where D <: Dual = valtype(D)
@inline valtype(::Type{Complex{D}}) where D <: Dual = valtype(D)

@inline basetype(::V) where {V} = V # fallback
@inline basetype(::Type{V}) where {V} = V # fallback
@inline basetype(::Dual{T,V,N}) where {T,V,N} = V
@inline basetype(::Type{Dual{T,V,N}}) where {T,V,N} = V
@inline basetype(::Dual{T,V,N}) where {T,V<:Dual,N} = basetype(V)
@inline basetype(::Type{Dual{T,V,N}}) where {T,V<:Dual,N} = basetype(V)
@inline basetype(::Complex{D}) where {D<:Dual} = basetype(D)
@inline basetype(::Type{Complex{D}}) where {D<:Dual} = basetype(D)

@inline numbases(::Type{D}) where {D<:Dual} = div(sizeof(D), sizeof(basetype(D)))
@inline numbases(::Type{Complex{D}}) where {D<:Dual} = numbases(D)

# ------------------------- #
# dual2array and array2dual #
# ------------------------- #
function dual2array(X::Array{D,Dim}) where {D<:Dual, Dim}
    X = reinterpret(basetype(D), X, (numbases(D), size(X)...))
    Y = copy(X)
    return Y
end
function dual2array(X::Array{Complex{D},Dim}) where {D<:Dual, Dim}
    Xr = reinterpret(basetype(D), real(X), (numbases(D), size(X)...))
    Xi = reinterpret(basetype(D), imag(X), (numbases(D), size(X)...))
    Z = complex.(Xr, Xi)
    return Z
end

function array2dual(::Type{D}, X::Array{T,Dim}) where {D<:Dual, T, Dim}
    X = reinterpret(D, X, size(X)[2:end])
    Y = copy(X)
    return Y
end
function array2dual(::Type{Complex{D}}, X::Array{Complex{T},Dim}) where {D<:Dual, T, Dim}
    Xr = reinterpret(D, real(X), size(X)[2:end])
    Xi = reinterpret(D, imag(X), size(X)[2:end])
    Z = complex.(Xr, Xi)
    return Z
end

# ----------------------------------------------------- #
# randn(Complex{D}) where partials of D are random, too #
# ----------------------------------------------------- #
@generated function randn_tuple(rng::Random.AbstractRNG, ::Type{NTuple{N,V}}) where {N,V}
    return ForwardDiff.tupexpr(i -> :(randn(rng, V)), N)
end

@generated function randn_tuple(::Type{NTuple{N,V}}) where {N,V}
    return ForwardDiff.tupexpr(i -> :(randn(V)), N)
end

# randn(::Partials)
@inline Base.randn(partials::Partials) = randn(typeof(partials))
@inline Base.randn(::Type{Partials{N,V}}) where {N,V} = Partials{N,V}(randn_tuple(NTuple{N,V}))
@inline Base.randn(rng::Random.AbstractRNG, partials::Partials) = randn(rng, typeof(partials))
@inline Base.randn(rng::Random.AbstractRNG, ::Type{Partials{N,V}}) where {N,V} = Partials{N,V}(randn_tuple(rng, NTuple{N,V}))

# randn(::Dual)
@inline Base.randn(d::Dual) = randn(typeof(d))
@inline Base.randn(::Type{Dual{T,V,N}}) where {T,V,N} = Dual{T}(randn(V), randn(Partials{N,V}))
@inline Base.randn(rng::Random.AbstractRNG, d::Dual) = randn(rng, typeof(d))
@inline Base.randn(rng::Random.AbstractRNG, ::Type{Dual{T,V,N}}) where {T,V,N} = Dual{T}(randn(rng, V), randn(Partials{N,V}))

# randn(::Complex{Dual})
@inline Base.randn(d::Complex{D}) where D <: Dual = randn(typeof(d))
@inline Base.randn(::Type{Complex{D}}) where D <: Dual = Complex(randn(D), randn(D))
@inline Base.randn(rng::Random.AbstractRNG, d::Complex{D}) where D <: Dual = randn(rng, typeof(d))
@inline Base.randn(rng::Random.AbstractRNG, ::Type{Complex{D}}) where D <: Dual = Complex(randn(rng, D), randn(rng, D))


text:
using DDF

using LinearAlgebra
using StaticArrays
using Test

# Dmax-1 since these tests are expensive
const Dmax1 = max(0, Dmax - 1)
@testset "Op D=$D P2=$P2 R2=$R2 P1=$P1 R1=$R1" for D in 0:Dmax1,
P2 in (Pr, Dl),
R2 in 0:D,
P1 in (Pr, Dl),
R1 in 0:D

    S = Float64
    T = Rational{Int64}
    mfd = hypercube_manifold(Val(D), S)

    z = zero(Op{D,P2,R2,P1,R1,T}, mfd)
    e = one(Op{D,P1,R1,P1,R1,T}, mfd)
    e2 = one(Op{D,P2,R2,P2,R2,T}, mfd)
    A = rand(Op{D,P2,R2,P1,R1,T}, mfd)
    B = rand(Op{D,P2,R2,P1,R1,T}, mfd)
    C = rand(Op{D,P2,R2,P1,R1,T}, mfd)
    f = rand(Fun{D,P1,R1,D,S,T}, mfd)
    g = rand(Fun{D,P1,R1,D,S,T}, mfd)
    a = rand(T)
    b = rand(T)

    @test iszero(z)
    @test !iszero(e)
    @test isone(e)
    @test !isone(z)

    # Vector space

    @test +A == A
    @test (A + B) + C == A + (B + C)

    @test z + A == A
    @test A + z == A

    @test A + (-A) == z
    @test (-A) + A == z
    @test A - B == A + (-B)

    @test (a * b) * A == a * (b * A)
    @test A * (a * b) == (A * a) * b

    @test one(T) * a == a
    @test one(T) * A == A
    @test A * one(T) == A

    @test zero(T) * a == zero(T)
    @test zero(T) * A == z
    @test A * zero(T) == z

    if a != zero(T)
        @test a * inv(a) == one(T)

        @test inv(a) * (a * A) == A
        @test (A * a) * inv(a) == A

        @test inv(a) * A == a \ A
        @test A * inv(a) == A / a
    end

    @test (a + b) * A == a * A + b * A
    @test A * (a + b) == A * a + A * b

    @test a * (A + B) == a * A + a * B
    @test (A + B) * a == A * a + B * a

    # Category

    @test e2 * A == A
    @test A * e == A

    for P3 in (Pr, Dl), R3 in 0:D, P4 in (Pr, Dl), R4 in 0:D
        F = rand(Op{D,P3,R3,P2,R2,T}, mfd)
        G = rand(Op{D,P4,R4,P3,R3,T}, mfd)

        @test (G * F) * A == G * (F * A)
    end

    # # Groupoid ("Division ring")
    #
    # if A != z
    #     @test A * inv(A) == e
    #     @test inv(A) * A == e
    #     @test B / A == B * inv(A)
    #     @test A \ B == inv(A) * B
    # end

    # Adjoints

    @test A'' == A

    # Operators act on functions

    @test A * (f + g) == A * f + A * g
    @test (A + B) * f == A * f + B * f

    for P3 in (Pr, Dl), R3 in 0:D
        F = rand(Op{D,P3,R3,P2,R2,T}, mfd)
        @test (F * A) * f == F * (A * f)
    end

    maxabs(f) = norm(f.values, Inf)

    isgood = false
    for iter in 1:3
        # Add diagonal entries to help make F invertible
        F = one(Op{D,P1,R1,P1,R1,T}, mfd) + rand(Op{D,P1,R1,P1,R1,T}, mfd)

        # Note: \ converts rationals to Float64
        g = F * (F \ f)
        isgood = maxabs(g - f) ≤ 1.0e-12
        isgood && break
    end
    @test isgood

    isgood = false
    for iter in 1:3
        # Add diagonal entries to help make F and G invertible
        F = one(Op{D,P1,R1,P1,R1,T}, mfd) + rand(Op{D,P1,R1,P1,R1,T}, mfd)
        G = one(Op{D,P1,R1,P1,R1,T}, mfd) + rand(Op{D,P1,R1,P1,R1,T}, mfd)

        # Note: \ converts rationals to Float64
        g1 = (G * F) \ f
        g2 = F \ (G \ f)
        gscale = max(1, maxabs(g1), maxabs(g2))
        # isgood = maxabs(g1 - g2) ≤ 1.0e-12 * gscale
        isgood = maxabs(g1 - g2) ≤ 1.0e-8 * gscale
        isgood && break
    end
    @test isgood
end


text:
#!/usr/bin/env julia
using XAM,Statistics,GenomicFeatures

bName = ARGS[1]
cName = ARGS[2]

function histQuals(bamName,chrName)
	allquals=Int[]
	bamIn = open(BAM.Reader,bamName,index="$bamName.bai")
	count = 0

	#for record in eachoverlap(reader, "Chr2", 10000:11000)
	for a in eachoverlap(bamIn,chrName,1:1000000)
		count+=1
    	push!(allquals,BAM.mappingquality(a))
	end

	close(bamIn)

	println("Count: ",count)
	println("Mean quality: ",mean(allquals))
end

histQuals(bName,cName)

# time averagequal.jl /Users/james/projects/hla/capture/bam38/PC-HLA-EG4.bam
# Reading BAM...done.
# Average quality: 43.76016208794757
# averagequal.jl /Users/james/projects/hla/capture/bam38/PC-HLA-EG4.bam
# 15.48s user 0.94s system 99% cpu 16.422 total

# time averagequal.jl data/PC-HLA-EG4.bam
# Reading BAM...done.
# Average quality: 56.078356406954285
# averagequal.jl data/PC-HLA-EG4.bam
# 1.85s user 0.42s system 108% cpu 2.105 total

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
function part1(input)
    instructions = parse_input(input)
    position = navigate1(instructions)
    return prod(position)
end

function part2(input)
    instructions = parse_input(input)
    position = navigate2(instructions)
    return prod(position)
end

parse_input(input) = parse_line.(readlines(input))

function parse_line(line)
    direction, distance = split(line, " ")
    return direction, parse(Int, distance)
end

function navigate1(instructions)
    position = (0, 0)
    for (direction, distance) in instructions
        if direction == "forward"
            position = position .+ (distance, 0)
        elseif direction == "down"
            position = position .+ (0, distance)
        elseif direction == "up"
            position = position .- (0, distance)
        end
    end
    return position
end

function navigate2(instructions)
    position = (0, 0)
    aim = 0
    for (direction, distance) in instructions
        if direction == "forward"
            position = position .+ (distance, aim * distance)
        elseif direction == "down"
            aim += distance
        elseif direction == "up"
            aim -= distance
        end
    end
    return position
end


text:
"""
    FullFeaturedCPNN(;
        graphChain::Flux.Chain
        nodeChain::Flux.Chain
        globalChain::Flux.Chain
        outputChain::Flux.Dense
    )

This structure is here to provide a flexible way to create a nn model which respect this approach:
Making modification on the graph, then extract one node feature and modify it.
"""
Base.@kwdef struct FullFeaturedCPNN <: NNStructure
    graphChain::Flux.Chain = Flux.Chain()
    nodeChain::Flux.Chain = Flux.Chain()
    globalChain::Flux.Chain = Flux.Chain()
    outputChain::Union{Flux.Dense, Flux.Chain} = Flux.Chain()
end

# Enable the `|> gpu` syntax from Flux
Flux.@functor FullFeaturedCPNN

# TODO make it possible to use global features
function (nn::FullFeaturedCPNN)(states::BatchedDefaultTrajectoryState)
    @assert !isnothing(states.allValuesIdx)

    variableIdx = states.variableIdx
    batchSize = length(variableIdx)
    allValuesIdx = states.allValuesIdx
    actionSpaceSize = size(allValuesIdx, 1)
    mask = device(states) == Val(:gpu) ? CUDA.zeros(Float32, 1, actionSpaceSize, batchSize) : zeros(Float32, 1, actionSpaceSize, batchSize) # this mask will replace `reapeat` using broadcasted `+`

    # chain working on the graph(s) with the GNNs
    featuredGraph = nn.graphChain(states.fg)
    nodeFeatures = featuredGraph.nf # FxNxB
    globalFeatures = featuredGraph.gf # GxB

    # Extract the features corresponding to the varibales
    variableIndices = nothing
    Zygote.ignore() do
        variableIndices = Flux.unsqueeze(CartesianIndex.(variableIdx, 1:batchSize), 1)
    end
    variableFeatures = nodeFeatures[:, variableIndices] # Fx1xB
    variableFeatures = reshape(nn.nodeChain(RL.flatten_batch(variableFeatures)), :, 1, batchSize) # F'x1xB

    # Extract the features corresponding to the values
    valueIndices = nothing
    Zygote.ignore() do
        valueIndices = CartesianIndex.(allValuesIdx, repeat(transpose(1:batchSize); outer=(actionSpaceSize, 1)))
    end
    valueFeatures = nodeFeatures[:, valueIndices] # FxAxB
    valueFeatures = reshape(nn.nodeChain(RL.flatten_batch(valueFeatures)), :, actionSpaceSize, batchSize) # F'xAxB

    finalFeatures = nothing
    if sizeof(globalFeatures) != 0

        # Extract the global features
        globalFeatures = reshape(nn.globalChain(globalFeatures), :, 1, batchSize) # G'x1xB

        # Prepare the input of the outputChain
        finalFeatures = vcat(
            variableFeatures .+ mask, # F'xAxB
            globalFeatures .+ mask, # G'xAxB
            valueFeatures,
        ) # (F'+G'+F')xAxB
        finalFeatures = RL.flatten_batch(finalFeatures) # (F'+G'+F')x(A+B)
    else
        # Prepare the input of the outputChain
        finalFeatures = vcat(
            variableFeatures .+ mask, # F'xAxB
            valueFeatures,
        ) # (F'+F')xAxB
        finalFeatures = RL.flatten_batch(finalFeatures) # (F'+F')x(A+B)
    end

    # output layer
    predictions = nn.outputChain(finalFeatures) # Ox(A+B)
    output = reshape(predictions, actionSpaceSize, batchSize) # OxAxB

    return output
end


text:
using StaticArrays
using Test

using ClimateMachine
ClimateMachine.init()
using ClimateMachine.Atmos
using ClimateMachine.Mesh.Grids
using ClimateMachine.MoistThermodynamics
using ClimateMachine.VariableTemplates

using CLIMAParameters
using CLIMAParameters.Planet: grav, MSLP
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

function init_test!(bl, state, aux, (x, y, z), t)
    FT = eltype(state)

    z = FT(z)
    _grav::FT = grav(bl.param_set)
    _MSLP::FT = MSLP(bl.param_set)

    # These constants are those used by Stevens et al. (2005)
    qref = FT(9.0e-3)
    q_pt_sfc = PhasePartition(qref)
    Rm_sfc = FT(gas_constant_air(param_set, q_pt_sfc))
    T_sfc = FT(290.4)
    P_sfc = _MSLP

    # Specify moisture profiles
    q_liq = FT(0)
    q_ice = FT(0)

    θ_liq = FT(289.0)
    q_tot = qref

    ugeo = FT(7)
    vgeo = FT(-5.5)
    u, v, w = ugeo, vgeo, FT(0)

    # Pressure
    H = Rm_sfc * T_sfc / _grav
    p = P_sfc * exp(-z / H)

    # Density, Temperature
    ts = LiquidIcePotTempSHumEquil_given_pressure(bl.param_set, θ_liq, p, q_tot)
    ρ = air_density(ts)

    e_kin = FT(1 / 2) * FT((u^2 + v^2 + w^2))
    e_pot = _grav * z
    E = ρ * total_energy(e_kin, e_pot, ts)

    state.ρ = ρ
    state.ρu = SVector(ρ * u, ρ * v, ρ * w)
    state.ρe = E
    state.moisture.ρq_tot = ρ * q_tot

    return nothing
end

function main()
    FT = Float64

    # DG polynomial order
    N = 4

    # Domain resolution and size
    Δh = FT(40)
    Δv = FT(40)
    resolution = (Δh, Δh, Δv)

    xmax = FT(320)
    ymax = FT(320)
    zmax = FT(400)

    t0 = FT(0)
    timeend = FT(10)
    CFL = FT(0.4)

    ode_solver = ClimateMachine.ExplicitSolverType()
    driver_config = ClimateMachine.AtmosLESConfiguration(
        "Driver test",
        N,
        resolution,
        xmax,
        ymax,
        zmax,
        param_set,
        init_test!,
        solver_type = ode_solver,
    )
    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        Courant_number = CFL,
    )

    # Test the courant wrapper
    # by default the CFL should be less than what asked for
    CFL_nondiff = ClimateMachine.DGmethods.courant(
        ClimateMachine.Courant.nondiffusive_courant,
        solver_config,
    )
    @test CFL_nondiff < CFL
    CFL_adv = ClimateMachine.DGmethods.courant(
        ClimateMachine.Courant.advective_courant,
        solver_config,
    )
    CFL_adv_v = ClimateMachine.DGmethods.courant(
        ClimateMachine.Courant.advective_courant,
        solver_config;
        direction = VerticalDirection(),
    )
    CFL_adv_h = ClimateMachine.DGmethods.courant(
        ClimateMachine.Courant.advective_courant,
        solver_config;
        direction = HorizontalDirection(),
    )

    # compute known advective Courant number (based on initial conditions)
    ugeo_abs = FT(7)
    vgeo_abs = FT(5.5)
    Δt = solver_config.dt
    ca_h = ugeo_abs * (Δt / Δh) + vgeo_abs * (Δt / Δh)
    # vertical velocity is 0
    caᵥ = FT(0.0)
    @test isapprox(CFL_adv_v, caᵥ)
    @test isapprox(CFL_adv_h, ca_h, atol = 0.0005)
    @test isapprox(CFL_adv, ca_h, atol = 0.0005)

    cb_test = 0
    result = ClimateMachine.invoke!(solver_config)
    # cb_test should be zero since user_info_callback not specified
    @test cb_test == 0

    result = ClimateMachine.invoke!(
        solver_config,
        user_info_callback = (init) -> cb_test += 1,
    )
    # cb_test should be greater than one if the user_info_callback got called
    @test cb_test > 0

    # Test that if dt is not adjusted based on final time the CFL is correct
    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        Courant_number = CFL,
        timeend_dt_adjust = false,
    )

    CFL_nondiff = ClimateMachine.DGmethods.courant(
        ClimateMachine.Courant.nondiffusive_courant,
        solver_config,
    )
    @test CFL_nondiff ≈ CFL
end

main()