text:
#
# Immutable dictionary implemented as a bitmap trie.
#

# Trie nodes.

abstract ImmDictNode{K,V}

immutable ValueNode{K,V} <: ImmDictNode{K,V}
    key::K
    val::V
    key_hc::UInt
end

immutable CollisionNode{K,V} <: ImmDictNode{K,V}
    keys::Vector{K}
    vals::Vector{V}
    key_hc::UInt
end

immutable BitmapNode{K,V} <: ImmDictNode{K,V}
    branches::Vector{ImmDictNode{K,V}}
    bitmap::UInt32

    BitmapNode(branches, bitmap) = new(branches, bitmap)
    BitmapNode() = new(ImmDictNode{K,V}[], 0x00000000)
end

# Lookup in trie nodes.

find{K,V}(node::ValueNode{K,V}, key, key_hc, shift) =
    key_hc == node.key_hc && isequal(key, node.key) ?
        Nullable{V}(node.val) : Nullable{V}()

function find{K,V}(node::CollisionNode{K,V}, key, key_hc, shift)
    if key_hc == node.key_hc
        for j = eachindex(node.keys)
            if isequal(key, node.keys[j])
                return Nullable{V}(node.vals[j])
            end
        end
    end
    return Nullable{V}()
end

mask(hc::UInt, shift::Int) = (hc >>> shift) & 0b11111
bitpos(mask::UInt) = 0x00000001 << mask

function find{K,V}(node::BitmapNode{K,V}, key, key_hc, shift)
    pos = bitpos(mask(key_hc, shift))
    if (node.bitmap & pos) != 0
        j = 1 + count_ones(node.bitmap & (pos-1))
        return find(node.branches[j], key, key_hc, shift+5)::Nullable{V}
    end
    return Nullable{V}()
end

# Iteration.

start{K,V}(node::ValueNode{K,V}) = true
next{K,V}(node::ValueNode{K,V}, flag) = (Pair{K,V}(node.key, node.val), false)
done{K,V}(node::ValueNode{K,V}, flag) = !flag

start{K,V}(node::CollisionNode{K,V}) = 1
next{K,V}(node::CollisionNode{K,V}, k) =
    (Pair{K,V}(node.keys[k], node.vals[k]), k+1)
done{K,V}(node::CollisionNode{K,V}, k) = k > length(node.keys)

function start{K,V}(node::BitmapNode{K,V})
    for j = eachindex(node.branches)
        b = node.branches[j]
        bst = start(b)
        if !done(b, bst)
            return (j, bst)
        end
    end
    return ()
end

function next{K,V}(node::BitmapNode{K,V}, state)
    j, bst = state
    b = node.branches[j]
    pair, bst = next(b, bst)
    if !done(b, bst)
        return (pair, (j, bst))
    end
    while j < length(node.branches)
        j += 1
        b = node.branches[j]
        bst = start(b)
        if !done(b, bst)
            return (pair, (j, bst))
        end
    end
    return (pair, ())
end

done{K,V}(node::BitmapNode{K,V}, state) = state == ()

# Adding an element.

type PutResult{V}
    old::Nullable{V}
    changed::Bool

    PutResult() = new(Nullable{V}(), false)
end

function put{K,V}(node::ValueNode{K,V}, key, val, key_hc, shift, result)
    if key_hc == node.key_hc
        if isequal(key, node.key)
            if isequal(val, node.val)
                return node
            else
                result.old = Nullable{V}(node.val)
                result.changed = true
                return ValueNode{K,V}(key, val, key_hc)
            end
        else
            result.changed = true
            return CollisionNode{K,V}(K[node.key, key], [node.val, val], key_hc)
        end
    else
        branches = ImmDictNode{K,V}[node]
        bitmap = bitpos(mask(node.key_hc, shift))
        return put(BitmapNode{K,V}(branches, bitmap), key, val, key_hc, shift, result)
    end
end

function put{K,V}(node::CollisionNode{K,V}, key, val, key_hc, shift, result)
    if key_hc == node.key_hc
        N = length(node.keys)
        for j = 1:N
            if isequal(key, node.keys[j])
                if isequal(val, node.vals[j])
                    return node
                else
                    result.old = Nullable{V}(node.vals[j])
                    result.changed = true
                    vals = copy(node.vals)
                    vals[j] = val
                    return CollisionNode{K,V}(node.keys, vals, node.key_hc)
                end
            end
        end
        result.changed = true
        keys = Vector{K}(N+1)
        vals = Vector{K}(N+1)
        for j = 1:N
            keys[j] = node.keys[j]
            vals[j] = node.vals[j]
        end
        keys[N+1] = key
        vals[N+1] = val
        return CollisionNode{K,V}(keys, vals, node.key_hc)
    else
        nodes = ImmDictNode{K,V}[node]
        bitmap = bitpos(mask(node.key_hc, shift))
        return put(BitmapNode{K,V}(branches, bitmap), key, val, key_hc, shift, result)
    end
end

function put{K,V}(node::BitmapNode{K,V}, key, val, key_hc, shift, result)
    pos = bitpos(mask(key_hc, shift))
    j = 1 + count_ones(node.bitmap & (pos-1))
    if (node.bitmap & pos) != 0
        j = 1 + count_ones(node.bitmap & (pos-1))
        branch = put(node.branches[j], key, val, key_hc, shift+5, result)
        if result.changed
            branches = copy(node.branches)
            branches[j] = branch
            return BitmapNode{K,V}(branches, node.bitmap)
        else
            return node
        end
    end
    result.changed = true
    branch = ValueNode{K,V}(key, val, key_hc)
    N = length(node.branches)
    branches = Vector{ImmDictNode{K,V}}(N+1)
    for k = 1:j-1
        branches[k] = node.branches[k]
    end
    branches[j] = branch
    for k = j:N
        branches[k+1] = node.branches[k]
    end
    bitmap = node.bitmap | pos
    return BitmapNode{K,V}(branches, bitmap)
end

# Immutable dictionary.

immutable ImmutableDict{K,V} <: Associative{K,V}
    root::ImmDictNode{K,V}
    sz::Int
    hc::UInt

    ImmutableDict() = new(BitmapNode{K,V}(), 0, 0)
    ImmutableDict(root::ImmDictNode{K,V}, sz::Int, hc::UInt) = new(root, sz, hc)
end

ImmutableDict{K,V}(kvs::Pair{K,V}...) =
    let d = ImmutableDict{K,V}()
        for (k, v) in kvs
            d = assoc(d, k, v)
        end
        d
    end

length(idict::ImmutableDict) = idict.sz
isempty(idict::ImmutableDict) = idict.sz == 0

start(idict::ImmutableDict) = start(idict.root)
next(idict::ImmutableDict, state) = next(idict.root, state)
done(idict::ImmutableDict, state) = done(idict.root, state)

get{K,V}(idict::ImmutableDict{K,V}, key) =
    find(idict.root, key, hash(key), 0)::Nullable{V}

get{K,V}(idict::ImmutableDict{K,V}, key, default) =
    let maybe_val = get(idict, key)
        isnull(maybe_val) ? default : get(maybe_val)
    end

getindex{K,V}(idict::ImmutableDict{K,V}, key) =
    let maybe_val = get(idict, key)
        isnull(maybe_val) ? throw(KeyError(key)) : get(maybe_val)
    end

function assoc{K,V}(idict::ImmutableDict{K,V}, key, val)
    result = PutResult{V}()
    new_root = put(idict.root, key, val, hash(key), 0, result)
    if result.changed
        entry_hc = hash(key, hash(val))
        if !isnull(result.old)
            old_entry_hc = hash(key, hash(get(result.old)))
            hc = idict.hc $ entry_hc $ old_entry_hc
            return ImmutableDict{K,V}(new_root, idict.sz, hc)
        else
            sz = idict.sz + 1
            hc = idict.hc $ entry_hc
            return ImmutableDict{K,V}(new_root, sz, hc)
        end
    end
    return idict
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
using Test

ProjDir = joinpath(dirname(@__FILE__), "..", "examples", "BernoulliInitTheta")
cd(ProjDir) do

  isdir("tmp") &&
    rm("tmp", recursive=true);

  bernoullimodel = "
  data {
    int<lower=1> N;
    int<lower=0,upper=1> y[N];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    theta ~ beta(1,1);
    y ~ bernoulli(theta);
  }
  "

  bernoullidata = [
    Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1]),
    Dict("N" => 10, "y" => [0, 1, 0, 0, 0, 0, 1, 0, 0, 1]),
    Dict("N" => 10, "y" => [0, 0, 0, 0, 0, 0, 1, 0, 1, 1]),
    Dict("N" => 10, "y" => [0, 0, 0, 1, 0, 0, 0, 1, 0, 1])
  ]

  theta_init = [0.1, 0.4, 0.5, 0.9]
  bernoulliinit = Dict[
    Dict("theta" => theta_init[1]),
    Dict("theta" => theta_init[2]),
    Dict("theta" => theta_init[3]),
    Dict("theta" => theta_init[4]),
  ]

  monitor = ["theta", "lp__", "accept_stat__"]

  stanmodel = Stanmodel(name="bernoulli",
    model=bernoullimodel,
    init=Stan.Init(init=bernoulliinit),
    adapt=1);

  rc, a3d, cnames = stan(stanmodel, bernoullidata, CmdStanDir=CMDSTAN_HOME)

  if rc == 0
    sdf  = read_summary(stanmodel)
    @test sdf[sdf.parameters .== :theta, :mean][1] ≈ 0.34 rtol=0.1
  end

  isdir("tmp") &&
    rm("tmp", recursive=true);

end # cd

ext: jl


text:
module SwarmIntelligence

end


text:
# FFSSolver: (F)orward (S)ingle (S)hooting (Solver)
# Solves indirect trajectory optimization problem using a forward
# single shooting based approach.

# NOTE: Likely only supports spacecraft trajectory optimization problems
# involving a single spacecraft using a 6 element state representation
# with mass (7 total elements).

# Abstraction layer to allow for disptaching based on ShootingDataManager
abstract type AbstractFSSSolver{SDMT <: AbstractShootingDataManager} end

struct FSSSolver{SDMT,NLEF} <: AbstractFSSSolver{SDMT}

    # Shooting Data Manager
    sdm::SDMT

    # Non-linear equations function
    nlEqs::NLEF

end

# Constructor
function FSSSolver(initGuess, ics, fcs, bvpFuncNoSTM, bvpFuncWSTM;
    homotopy = false, homotopyParamVec = [1.0])

    # Get size of variables vector
    m = length(initGuess)
    n = 2*m

    # Set up ShootingDataManager
    if homotopy == false
        sdm = ShootingDataManager()
    else
        sdm = ShootingHomotopyDataManager()
        SetHomotopyParams!(sdm, homotopyParamVec)
    end
    SetInitGuessData!(sdm, initGuess)
    SetInitConditions!(sdm, ics)
    SetFinConditions!(sdm, fcs)
    SetPreallocatedVecs!(sdm, n + n*n)
    Initialize!(sdm)

    # Configuring input vector
    inVec = GetInputVec(sdm)
    inVec[1:m]      .= ics
    inVec[m+1:n]    .= initGuess
    inVec[n+1:end]  .= 0.0
    for i in 1:n
        inVec[n + i + (i-1)*n] = 1.0
    end

    # Configure NLsolve.jl input function
    if homotopy == false
        nlEqs = (F,J,λ) -> fssNLEqFunc!(F, J, λ, bvpFuncNoSTM, bvpFuncWSTM, sdm)
    else
        nlEqs = (F,J,λ,ϵ) -> fssHomotopyNLEqFunc!(F, J, λ, ϵ, bvpFuncNoSTM, bvpFuncWSTM, sdm)
    end

    FSSSolver(sdm, nlEqs)
end

function FSSSolver(ics, fcs, bvpFuncNoSTM, bvpFuncWSTM;
    homotopy = false, homotopyParamVec = [1.0])

    # Set up ShootingDataManager
    if homotopy == false
        sdm = ShootingDataManager()
    else
        sdm = ShootingHomotopyDataManager()
        SetHomotopyParams!(sdm, homotopyParamVec)
    end
    SetInitConditions!(sdm, ics)
    SetFinConditions!(sdm, fcs)

    # Configure NLsolve.jl input function
    if homotopy == false
        nlEqs = (F,J,λ) -> fssNLEqFunc!(F, J, λ, bvpFuncNoSTM, bvpFuncWSTM, sdm)
    else
        nlEqs = (F,J,λ,ϵ) -> fssHomotopyNLEqFunc!(F, J, λ, ϵ, bvpFuncNoSTM, bvpFuncWSTM, sdm)
    end

    FSSSolver(sdm, nlEqs)
end

function initializeData!(solver::FSSSolver, initGuess::AbstractVector)
    m = length(initGuess)
    n = 2*m
    SetInitGuessData!(solver.sdm, initGuess)
    SetPreallocatedVecs!(solver.sdm, n + n*n)
    Initialize!(solver.sdm)
    inVec = GetInputVec(solver.sdm)
    inVec[1:m]      .= solver.sdm.ics
    inVec[m+1:n]    .= initGuess
    inVec[n+1:end]  .= 0.0
    for i in 1:n
        inVec[n + i + (i-1)*n] = 1.0
    end

    return nothing
end

function solve!(solver::AbstractFSSSolver{ShootingDataManager}; factor = 1.0, ftol = 1e-8, showTrace = true,convergenceAttempts = 4)
    sdm = solver.sdm
    sol = nlsolve(only_fj!(solver.nlEqs), GetInitGuessData(sdm);
        show_trace = showTrace, factor = factor, ftol =ftol)
    solver.sol .= sol.zero
    sdm.cflag = (sol.residual_norm < ftol*100.0)
    return nothing
end

function solve!(solver::AbstractFSSSolver{ShootingHomotopyDataManager}; factor = 1.0, ftol = 1e-8, showTrace = true, convergenceAttempts = 4)
    sdm         = solver.sdm
    ϵs          = GetHomotopyParams(sdm)
    sols        = GetHomotopySolutionVector(sdm)
    cflags      = GetHomotopyConvergenceFlags(sdm)
    lastConvIdx = 1
    cFlag       = false
    @inbounds for i in 1:length(ϵs)
        # Get homotopy parameter and solve
        ϵ = ϵs[i]
        if i == 1
            sol = nlsolve(only_fj!((F,J,λ)->solver.nlEqs(F,J,λ,ϵ)),
                GetInitGuessData(sdm); show_trace = showTrace, factor = factor, ftol = ftol)
        else
            println("Shooting for ϵ = " * string(ϵ))
            if cFlag == true
                lastConvIdx = i - 1
            end
            sol = nlsolve(only_fj!((F,J,λ)->solver.nlEqs(F,J,λ,ϵ)),
                sols[lastConvIdx]; show_trace = true, factor = factor, ftol = ftol)
        end

        # Save solution
        cFlag = (sol.residual_norm < ftol*100.0)
        cflags[i] = cFlag
        sols[i] .= sol.zero

        # Edge case catch: If final homotopy step fails, attempt to correct
        # by reducing steps.
        if i == length(ϵs) && cflags[end] == false && cflags[end-1] == true
            println("Last step failed. Attempting to correct...")
            attempts    = 0
            success     = false
            localSol    = Vector{Float64}(undef, length(sols[1]))
            while success == false && attempts <= convergenceAttempts
                attempts    += 1
                numSteps    = 2 + attempts
                localSol    .= sols[end-1]
                for ϵ in LinRange(ϵs[end-1], ϵs[end], numSteps)[2:end] # This is allocating twice. Might want to improve but hitting this edge case is rare...
                    println("Shooting for ϵ = " * string(ϵ))
                    sol = nlsolve(only_fj!((F,J,λ)->solver.nlEqs(F,J,λ,ϵ)),
                        localSol; show_trace = true, factor = factor, ftol = ftol)
                    cFlag = (sol.residual_norm < ftol*100.0)
                    if cFlag == false
                        break
                    else
                        localSol .= sol.zero
                    end
                end
                if cFlag == true
                    success = true
                    sols[end]  .= localSol
                    cflags[end] = true
                end
            end
        end

        # Check if convergence has stalled
        if cFlag == false && (i - lastConvIdx) >= convergenceAttempts
            println("Convergence failed during continuation.")
            break
        end
    end
end

# NLsolve.jl function
function fssNLEqFunc!(F::Union{Nothing,AbstractArray}, J::Union{Nothing,AbstractArray},
                      λ::AbstractArray, bvpFuncNoSTM, bvpFuncWSTM,
                      sdm::ShootingDataManager)

    # Get data from ShootingDataManager
    z0 = GetInputVec(sdm)
    zf = GetOutputVec(sdm)
    fc = GetFinConditions(sdm)

    # Fill J with zeros if necessary
    if !(J === nothing)
        fill!(J, 0.0)
    end

    # Added generality to hopfully make current revision applicable to
    # more problems
    m = length(λ)
    n = 2*m

    # Evaluate boundary value problem and fill data arrays
    if J === nothing && F !== nothing
        # Initialize full initial condition vector w/o STM
        z0[m+1:n] .= λ

        # Evaluate boundary value function
        @views zf[1:n] .= bvpFuncNoSTM(z0[1:n])

        # Fill F (This is not generalized!)
        @views F[1:6] .= zf[1:6] .- fc[1:6]
        F[7] = zf[14] - fc[7]
    elseif J !== nothing
        # Initialize full initial condition vector w/ STM
        z0[m+1:n] .= λ
n
        # Evaluate boundary value function with STM
        zf .= bvpFuncWSTM(z0)

        # Fill F (This is not generalized!)
        if F !== nothing
            @views F[1:6] .= zf[1:6] .- fc[1:6]
            F[7] = zf[14] - fc[7]
        end

        # Fill J
        @views STM = reshape(zf[n+1:end], (n,n))
        @views J .= STM[[1,2,3,4,5,6,14], 8:14] # Not generalized!
    end

    return nothing
end

# NLsolve.jl function
function fssHomotopyNLEqFunc!(F::Union{Nothing,AbstractArray}, J::Union{Nothing,AbstractArray},
                              λ::AbstractArray, ϵ, bvpFuncNoSTM, bvpFuncWSTM,
                              sdm::ShootingHomotopyDataManager)

    # Get data from ShootingDataManager
    z0 = GetInputVec(sdm)
    zf = GetOutputVec(sdm)
    fc = GetFinConditions(sdm)

    # Fill J with zeros if necessary
    if !(J === nothing)
        fill!(J, 0.0)
    end

    # Added generality to hopfully make current revision applicable to
    # more problems
    m = length(λ)
    n = 2*m

    # Evaluate boundary value problem and fill data arrays
    if J === nothing && F !== nothing
        # Initialize full initial condition vector w/o STM
        z0[m+1:n] .= λ

        # Evaluate boundary value function
        @views zf[1:n] .= bvpFuncNoSTM(z0[1:n], ϵ)

        # Fill F (This is not generalized!)
        @views F[1:6] .= zf[1:6] .- fc[1:6]
        F[7] = zf[14] - fc[7]
    elseif J !== nothing
        # Initialize full initial condition vector w/ STM
        z0[m+1:n] .= λ
n
        # Evaluate boundary value function with STM
        zf .= bvpFuncWSTM(z0, ϵ)

        # Fill F (This is not generalized!)
        if F !== nothing
            @views F[1:6] .= zf[1:6] .- fc[1:6]
            F[7] = zf[14] - fc[7]
        end

        # Fill J
        @views STM = reshape(zf[n+1:end], (n,n))
        @views J .= STM[[1,2,3,4,5,6,14], 8:14] # Not generalized!
    end

    return nothing
end

function GetInitialGuessConverged(fsss::FSSSolver)
    return GetInitialGuessConverged(fsss.sdm)
end

function GetHomotopyConverged(fsss::AbstractFSSSolver{ShootingHomotopyDataManager})
    return GetHomotopyConverged(fsss.sdm)
end

function GetSolution(fsss::AbstractFSSSolver{ShootingDataManager})
    return GetSolution(fsss.sdm)
end

function GetHomotopySolutionVector(fsss::AbstractFSSSolver{ShootingHomotopyDataManager})
    return GetHomotopySolutionVector(fsss.sdm)
end

function GetHomotopyParams(fsss::AbstractFSSSolver{ShootingHomotopyDataManager})
    return GetHomotopyParams(fsss.sdm)
end

function GetHomotopyConvergenceFlags(fsss::AbstractFSSSolver{ShootingHomotopyDataManager})
    return GetHomotopyConvergenceFlags(fsss.sdm)
end


text:
struct Analysis{S, F, NT<:NamedTuple}
    f::F
    kwargs::NT
    Analysis{S}(f::F, kwargs::NT) where {S, F, NT<:NamedTuple} =
        new{S, F, NT}(f, kwargs)
end

Analysis{S}(f; kwargs...) where {S} = Analysis{S}(f, values(kwargs))
Analysis{S}(a::Analysis; kwargs...) where {S} = Analysis{S}(a.f, merge(a.kwargs, values(kwargs)))
Analysis(args...; kwargs...) = Analysis{:automatic}(args...; kwargs...)

getfunction(funcs, S) = funcs
getfunction(funcs::NamedTuple, S::Symbol) = getfield(funcs, S)
getfunction(a::Analysis{S}) where {S} = getfunction(a.f, S)

(a::Analysis{S})(; kwargs...) where {S} = Analysis{S}(a; kwargs...)
(a::Analysis{:automatic})(; kwargs...) = Analysis{:automatic}(a; kwargs...)
(a::Analysis{S})(args...) where {S} = getfunction(a)(args...; a.kwargs...)
(a::Analysis{:automatic})(args...) = (infer_axis(args...)(a))(args...)

discrete(a::Analysis) = Analysis{:discrete}(a.f, a.kwargs)
continuous(a::Analysis) = Analysis{:continuous}(a.f, a.kwargs)
vectorial(a::Analysis) = Analysis{:vectorial}(a.f, a.kwargs)
discrete(::Nothing) = nothing
continuous(::Nothing) = nothing
vectorial(::Nothing) = nothing

Base.get(a::Analysis, s::Symbol, def) = get(a.kwargs, s, def)
Base.get(f::Function, a::Analysis, s::Symbol) = get(f, a.kwargs, s)
function set(f::Function, a::Analysis, s::Symbol)
    val = get(f, a, s)
    nt = NamedTuple{(s,)}((val,))
    a(; nt...)
end

const FunctionOrAnalysis = Union{Function, Analysis}

# TODO compute axis if called standalone!
compute_axis(f::Function, args...) = compute_axis(Analysis(f), args...)

infer_axis(x::AbstractVector{T}, args...) where {T<:Union{Missing, Number}} = Analysis{:continuous}
infer_axis(x::AbstractVector{T}, args...) where {T<:Union{Missing, AbstractArray}} = Analysis{:vectorial}
infer_axis(x::AbstractVector{Missing}, args...) = error("All data is missing")
infer_axis(x, args...) = Analysis{:discrete}

get_axis(s::Analysis) = s.kwargs.axis

function compute_axis(a::Analysis, args...)
    a_inf = infer_axis(args...)(a)
    compute_axis(a_inf, args...)
end

continuous_axis(x, args...; npoints = 100) = range(extrema(x)...; length = npoints)

function compute_axis(a::Analysis{:continuous}, args...)
    set(a, :axis) do
        npoints = get(a, :npoints, 100)
        continuous_axis(args..., npoints = npoints)
    end
end

discrete_axis(x, args...) = unique(sort(x))

function compute_axis(a::Analysis{:discrete}, args...)
    set(a, :axis) do
        discrete_axis(args...)
    end
end

vectorial_axis(x, args...) = axes(x[1], 1)

function compute_axis(a::Analysis{:vectorial}, args...)
    set(a, :axis) do
        vectorial_axis(args...)
    end
end

has_error(a::Analysis) = has_error(getfunction(a))
has_error(a::Analysis, args...) = has_error(infer_axis(args...)(a))
has_error(args...) = false

_expectedvalue(x, y; axis = nothing, kwargs...) = lazy_summary(x, y; kwargs...)

has_error(::typeof(_expectedvalue)) = true

function _localregression(x, y; npoints = 100, axis = continuous_axis(x, y; npoints = npoints), kwargs...)
    min, max = extrema(x)
    model = loess(convert(Vector{Float64}, x), convert(Vector{Float64}, y); kwargs...)
    return ((val, predict(model, val)) for (ind, val) in enumerate(axis) if min < val < max)
end

function _alignedsummary(xs, ys; axis = vectorial_axis(xs, ys), min_nobs = 2, stats = summary, kwargs...)
    stat, func = initstat(stats; kwargs...)
    iter = (view(y, x) for (x, y) in zip(xs, ys))
    stats = OffsetArray([copy(stat) for _ in axis], axis)
    fitvecmany!(stats, iter)
    ((val, func(st)) for (val, st) in zip(axis, stats) if nobs(st) >= min_nobs)
end

has_error(::typeof(_alignedsummary)) = true

const prediction = Analysis((continuous = _localregression, discrete = _expectedvalue, vectorial = _alignedsummary))

function _density(x; npoints = 100, axis = continuous_axis(x, npoints = npoints), kwargs...)
    d = InterpKDE(kde(x; kwargs...))
    return ((val, pdf(d, val)) for val in axis)
end

function _frequency(x; axis = discrete_axis(x), npoints = 100)
    c = countmap(x)
    s = sum(values(c))
    return ((val, get(c, val, 0)/s) for val in axis)
end

const density = Analysis((continuous = _density, discrete = _frequency))

function _cumulative(x; npoints = 100, axis = continuous_axis(x, npoints = npoints), kwargs...)
    func = ecdf(x)
    return ((val, func(val)) for val in axis)
end

function _cumulative_frequency(x; axis = discrete_axis(x, npoints = npoints), kwargs...)
    func = ecdf(x)
    return ((val, func(val)) for val in axis)
end

const cumulative = Analysis((continuous = _cumulative, discrete = _cumulative_frequency))

function _hazard(t; npoints = 100, axis = continuous_axis(t, npoints = npoints), kwargs...)
    pdf_iter = _density(t; axis = axis, kwargs...)
    cdf_func = ecdf(t)
    ((val, pdf / (1 + step(axis) * pdf - cdf_func(val))) for (val, pdf) in pdf_iter)
end

function _hazard_frequency(t; axis = discrete_axis(t), kwargs...)
    pdf_iter = _frequency(t; axis = axis, kwargs...)
    cdf_func = ecdf(t)
    ((val, pdf / (1 + pdf - cdf_func(val))) for (val, pdf) in pdf_iter)
end

const hazard = Analysis((continuous = _hazard, discrete = _hazard_frequency))


text:
using Documenter, ProcessFlows

makedocs(sitename="Manual for ProcessFlows.jl")

max_stars_repo_name: pbayer/pflow

max_issues_repo_name: pbayer/pflow

max_forks_repo_name: pbayer/pflow


text:
if !@isdefined(Femshop)
    include("../Femshop.jl");
    using .Femshop
end
init_femshop("FVadvection1d");

useLog("FVadvection1dlog")

# Configuration setup
domain(1)
solverType(FV)
timeStepper(RK4)

# Mesh
n = 40 # number of elements
mesh(LINEMESH, elsperdim=n, bids=2)

# Variables and BCs
u = variable("u", SCALAR, CELL)
v = variable("v", SCALAR, CELL)
boundary(u, 1, FLUX, "t<0.2 ? 1 : 0")
boundary(u, 2, NO_BC)
boundary(v, 1, FLUX, "t<0.2 ? 1 : 0")
boundary(v, 2, NO_BC)

# Time interval and initial condition
T = 0.5;
timeInterval(T)
initial(u, "0")
initial(v, "0")

# The flux and source terms of the conservation equation
# F and S in the following equation:
# Dt(int(u dx)) = int(S dx) - int(F.n ds)
coefficient("a", 1) # advection velocity
# The "upwind" function applies upwinding to the term (a.n)*u with flow velocity a.
# The optional third parameter is for tuning. Default upwind = 0, central = 1. Choose something between these.
flux([u, v], ["upwind(a,u)", "upwind(a,v,0.75)"])
# Note that there is no source() for this problem.

#@exportCode("fvad1dcode") # uncomment to export generated code to a file
#@importCode("fvad1dcode") # uncomment to import code from a file

solve([u,v])

finalize_femshop()

##### Uncomment below to compare to exact solution

# # The exact solution with constant velocity v
# a = 1;
# exact = zeros(n);
# x = Femshop.fv_info.cellCenters[:]
# for i=1:n
#     xt = x[i] - a*T;
#     if xt < 0
#         exact[i] = xt > -0.2 ? 1 : 0
#     end
# end

# using Plots
# pyplot();
# display(plot([x x x], [exact u.values[:] v.values[:]], markershape=:circle, label=["exact" "upwind" "alpha=0.75"]))


text:
#=
mstep.jl

    Maximization step of the Expectation-Maximization (EM) algorithm to estimate
    the hyper parameters of a linear Gaussian State Space model.

@author: Quint Wiersma <q.wiersma@vu.nl>

@date: 2022/02/09
=#

"""
    update_model!(model, state, smoother, fixed)

Conditional update of model components of the linear Gaussian state space model.

#### Arguments
  - `state::EMOptimizerState` 	: state variables
  - `smoother::Smoother`		: Kalman smoother output
  - `fixed::NamedTuple` 		: fixed hyper parameters

#### Returns
  - `model::StateSpaceModel`: state space model
"""
function update_model! end

"""
    mstep!(state, model, smoother, fixed,; ϵ_abs=1e-7, ϵ_rel=1e-3, max_iter=1000)

Update state and model using the maximization step of the EM algorithm.

#### Arguments
  - `state::EMState`        : state variables
  - `model::StateSpaceModel`: state space model
  - `smoother::Smoother`    : Kalman smoother output
  - `fixed::NamedTuple`     : fixed hyper parameters
  - `ϵ_abs::Real`           : absolute tolerance
  - `ϵ_rel::Real`           : relative tolerance
  - `max_iter::Integer`     : max number of iterations
"""
function mstep!(state::EMState, model::StateSpaceModel, smoother::Smoother,
                fixed::NamedTuple;
                ϵ_abs::Real=1e-7, ϵ_rel::Real=1e-3, max_iter::Integer=1000)
    # Get dims
    T_len= size(model.y,2)

    # lag and lead smoothed states
    α_lag= view(smoother.α,:,1:T_len-1)
    α_lead= view(smoother.α,:,2:T_len)

    # Update buffer variables
    # sum
    sum!(state.V_sum, view(smoother.V, :, :, 1, :))
    # variance
    # lead
    state.V_0.= state.V_sum .- view(smoother.V,:,:,1,1)
    mul!(state.V_0, α_lead, transpose(α_lead), 1., 1.)
    # lag
    state.V_1.= state.V_sum .- view(smoother.V,:,:,1,T_len)
    mul!(state.V_1, α_lag, transpose(α_lag), 1., 1.)
    # covariance
    sum!(state.V_01, view(smoother.V, :, :, 2, 2:T_len))
    mul!(state.V_01, α_lead, transpose(α_lag), 1., 1.)

    # Initialize stopping flags
    abs_change= Inf
    rel_change= Inf
    # Initialize iteration counter
    iter= 1
    # M-step (maximization)
    while (abs_change > ϵ_abs && rel_change > ϵ_rel) && iter < max_iter
        # Store current parameters
        copyto!(state.ψ_prev_m, state.ψ)

        # Update model
        update_model!(model, state, smoother, fixed)

        # Update state
        get_params!(state.ψ, model)

        # Store change in state
        state.Δ.= state.ψ .- state.ψ_prev_m

        # Absolute change
        abs_change= maximum(abs, state.Δ)
        # Relative change
        rel_change= abs_change * inv(1 + maximum(abs, state.ψ))

        # Update iteration counter
        iter+=1
    end

    return nothing
end

"""
    mstep!(state, model, smoother, fixed)

Update state and model using the conditional maximization step of the ECM
algorithm.

#### Arguments
  - `state::ECMState`       : state variables
  - `model::StateSpaceModel`: state space model
  - `smoother::Smoother`    : Kalman smoother output
  - `fixed::NamedTuple`     : fixed hyper parameters
"""
function mstep!(state::ECMState, model::StateSpaceModel, smoother::Smoother,
                pen::NamedTuple)
    # Get dims
    T_len= size(model.y,2)

    # lag and lead smoothed states
    α_lag= view(smoother.α,:,1:T_len-1)
    α_lead= view(smoother.α,:,2:T_len)

    # Update buffer variables
    # sum
    sum!(state.V_sum, view(smoother.V, :, :, 1, :))
    # variance
    # lead
    state.V_0.= state.V_sum .- view(smoother.V,:,:,1,1)
    mul!(state.V_0, α_lead, transpose(α_lead), 1., 1.)
    # lag
    state.V_1.= state.V_sum .- view(smoother.V,:,:,1,T_len)
    mul!(state.V_1, α_lag, transpose(α_lag), 1., 1.)
    # covariance
    sum!(state.V_01, view(smoother.V, :, :, 2, 2:T_len))
    mul!(state.V_01, α_lead, transpose(α_lag), 1., 1.)

    # Update model
    update_model!(model, state, smoother, pen)

    # Update state
    get_params!(state.ψ, model)

    return nothing
end


text:
module DSAModels

export AbstractShockAccelerationEfficiency,
    # efficiency models
    Kang07, KR13, CS14, Ryu19, P16,
    # mach nummber dependent efficiency
    η_Ms, η_Ms_acc, η_Ms_reacc,
    # B-field dependent efficiency
    η_B,
    ηB_acc_e, ηB_acc_p, ηB_reacc_e, ηB_reacc_p

"""
    AbstractShockAccelerationEfficiency

Abstract type for shock acceleration efficiencies
"""
abstract type AbstractShockAccelerationEfficiency end

"""
    Helper functions
"""
# power functions
@inline p2(x::T) where {T} = x * x
@inline p3(x::T) where {T} = x * x * x

"""
    kr_fitting_function(x::Real, a0::Real, a1::Real, a2::Real, a3::Real, a4::Real)

Helper function to use the fitting function from KR07.
"""
@inline function kr_fitting_function(x::Real,
    a0::Real, a1::Real, a2::Real, a3::Real, a4::Real)
    mm = x - 1
    return (a0 + a1 * mm + a2 * p2(mm) + a3 * p3(mm) + a4 * p2(p2(mm))) / p2(p2(x))
end

"""
    η_Ms(η_model::AbstractShockAccelerationEfficiency, M::T, X_cr::T) where {T<:Real}

Calculates the sonic Mach Number dependent η_model as a linear interpolation between acceleration and reacceleration.
"""
function η_Ms(η_model::AbstractShockAccelerationEfficiency, M::T, X_cr::T) where {T<:Real}

    # check if pressure ratio is greater than the target ratio
    if X_cr > η_model.X_cr
        X_cr = η_model.X_cr
    end

    # initial acceleration
    η1 = η_Ms_acc(η_model, M)
    X1::T = 0

    # reacceleration
    η2 = η_Ms_reacc(η_model, M)
    X2 = η_model.X_cr

    # linear interpolation between the two cases:
    η_tot::T = η1 + (X_cr - X1) * (η2 - η1) / (X2 - X1)

    # check if η_tot is smaller than the maximum η_model
    if η_tot > η_model.η_max
        return η_model.η_max
    else
        if η_tot < 0
            return T(0)
        else
            return η_tot
        end
    end

end

include("mach_models/Kang07.jl")
include("mach_models/KR13.jl")
include("mach_models/CS14.jl")
include("mach_models/Ryu19.jl")
include("mach_models/Pfrommer16.jl")

include("B_models/pais.jl")

end # module


text:
function csdSum{T<:AbstractFloat}(x::Vector{T})
    xlen = length(x)
    y = zeros(T,xlen)
    y[1] = x[1]
    for i in 2:length(x)
        y[i], y[i-1] = eftAdd(x[i], x[i-1])
    end
    sum(y[1:(end-1)]) + y[end]              # as if rounded from 2x significant bits
end

function csdDot{T<:AbstractFloat}(x::Vector{T}, y::Vector{T})
    xlen = length(x)
    if xlen != length(y)
       throw(ErrorException("vector lengths differ")
    elseif xlen==0
       return(x)
    end

    sHi, compensation = eftMul(x[1], y[1])
    for i in 2:xlen
        pHi, pLo = eftMul(x[i], y[i])
        sHi, sLo = eftAdd(sHi, pHi)
        compensation += pLo + sLo
    end

    sHi+compensation                          # as if rounded from 2x significant bits
end


text:
immutable Binomial <: DiscreteUnivariateDistribution
    n::Int
    p::Float64

    function Binomial(n::Int, p::Float64)
        n >= 0 || error("n must be non-negative.")
        0.0 <= p <= 1.0 || error("p must be in [0, 1]")
        new(n, p)
    end

    Binomial(n::Integer, p::Real) = Binomial(int(n), float64(p))
    Binomial(n::Integer) = Binomial(int(n), 0.5)
    Binomial() = new(1, 0.5)
end

@distr_support Binomial 0 d.n

@_jl_dist_2p Binomial binom

#### Parameters

ntrials(d::Binomial) = d.n
succprob(d::Binomial) = d.p
failprob(d::Binomial) = 1.0 - d.p

params(d::Binomial) = (d.n, d.p)

#### Properties

mean(d::Binomial) = ntrials(d) * succprob(d)
var(d::Binomial) = ntrials(d) * succprob(d) * failprob(d)
mode(d::Binomial) = ((n, p) = params(d); n > 0 ? round(Int,(n + 1) * d.prob) : 0)
modes(d::Binomial) = Int[mode(d)]

median(d::Binomial) = round(Int,mean(d))

function skewness(d::Binomial)
    n, p1 = params(d)
    p0 = 1.0 - p1
    (p0 - p1) / sqrt(n * p0 * p1)
end

function kurtosis(d::Binomial)
    n, p = params(d)
    u = p * (1.0 - p)
    (1.0 - 6.0 * u) / (n * u)
end

function entropy(d::Binomial; approx::Bool=false)
    n, p1 = params(d)
    (p1 == 0.0 || p1 == 1.0 || n == 0) && return 0.0
    p0 = 1.0 - p1
    if approx
        return 0.5 * (log(twoπ * n * p0 * p1) + 1.0)
    else
        lg = log(p1 / p0)
        lp = n * log(p0)
        s = exp(lp) * lp
        for k = 1:n
           lp += log((n - k + 1) / k) + lg
           s += exp(lp) * lp
        end
        return -s
    end
end

#### Evaluation

immutable RecursiveBinomProbEvaluator <: RecursiveProbabilityEvaluator
    n::Int
    coef::Float64   # p / (1 - p)
end

RecursiveBinomProbEvaluator(d::Binomial) = RecursiveBinomProbEvaluator(d.n, d.p / (1.0 - d.p))
nextpdf(s::RecursiveBinomProbEvaluator, p::Float64, x::Integer) = ((s.n - x + 1) / x) * s.coef * p
_pdf!(r::AbstractArray, d::Binomial, rgn::UnitRange) = _pdf!(r, d, rgn, RecursiveBinomProbEvaluator(d))

function mgf(d::Binomial, t::Real)
    n, p = params(d)
    (1.0 - p + p * exp(t)) ^ n
end

function cf(d::Binomial, t::Real)
    n, p = params(d)
    (1.0 - p + p * exp(im * t)) ^ n
end

#### Fit model

immutable BinomialStats <: SufficientStats
    ns::Float64   # the total number of successes
    ne::Float64   # the number of experiments
    n::Int        # the number of trials in each experiment

    function BinomialStats(ns::Real, ne::Real, n::Integer)
        new(float64(ns), float64(ne), int(n))
    end
end

function suffstats{T<:Integer}(::Type{Binomial}, n::Integer, x::AbstractArray{T})
    ns = zero(T)
    for i = 1:length(x)
        @inbounds xi = x[i]
        0 <= xi <= n || throw(DomainError())
        ns += xi
    end
    BinomialStats(ns, length(x), n)
end

function suffstats{T<:Integer}(::Type{Binomial}, n::Integer, x::AbstractArray{T}, w::AbstractArray{Float64})
    ns = 0.
    ne = 0.
    for i = 1:length(x)
        @inbounds xi = x[i]
        @inbounds wi = w[i]
        0 <= xi <= n || throw(DomainError())
        ns += xi * wi
        ne += wi
    end
    BinomialStats(ns, ne, n)
end

suffstats{T<:Integer}(::Type{Binomial}, data::(Int, Array{T})) = suffstats(Binomial, data...)
suffstats{T<:Integer}(::Type{Binomial}, data::(Int, Array{T}), w::Array{Float64}) = suffstats(Binomial, data..., w)

fit_mle(::Type{Binomial}, ss::BinomialStats) = Binomial(ss.n, ss.ns / (ss.ne * ss.n))

fit_mle{T<:Integer}(::Type{Binomial}, n::Integer, x::Array{T}) = fit_mle(Binomial, suffstats(Binomial, n, x))
fit_mle{T<:Integer}(::Type{Binomial}, n::Integer, x::Array{T}, w::Array{Float64}) = fit_mle(Binomial, suffstats(Binomial, n, x, w))
fit_mle{T<:Integer}(::Type{Binomial}, data::(Int, Array{T})) = fit_mle(Binomial, suffstats(Binomial, data))
fit_mle{T<:Integer}(::Type{Binomial}, data::(Int, Array{T}), w::Array{Float64}) = fit_mle(Binomial, suffstats(Binomial, data, w))


text:
# ---
# title: 790. Domino and Tromino Tiling
# id: problem790
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/domino-and-tromino-tiling/description/>
# hidden: true
# ---
#
# We have two types of tiles: a 2x1 domino shape, and an "L" tromino shape.
# These shapes may be rotated.
#
#
#
#     XX  <- domino
#
#     XX  <- "L" tromino
#     X
#
#
# Given N, how many ways are there to tile a 2 x N board? **Return your answer
# modulo 10^9 + 7**.
#
# (In a tiling, every square must be covered by a tile. Two tilings are
# different if and only if there are two 4-directionally adjacent cells on the
# board such that exactly one of the tilings has both squares occupied by a
# tile.)
#
#
#
#     **Example:**
#     Input: 3
#     Output: 5
#     Explanation:
#     The five different ways are listed below, different letters indicates different tiles:
#     XYZ XXZ XYY XXY XYY
#     XYZ YYZ XZZ XYY XXY
#
# **Note:**
#
#   * N  will be in range `[1, 1000]`.
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function TA_dijkstra_shortest_paths(graph, travel_time, origin, init_node, term_node, first_thru_node)
    no_node = nv(graph)
    no_arc = ne(graph)

    distmx = Inf*ones(no_node, no_node)
    for i in 1:no_arc
      if term_node[i] >= first_thru_node
          distmx[init_node[i], term_node[i]] = travel_time[i]
      end
    end

    state = dijkstra_shortest_paths(graph, origin, distmx)
    return state
end

function TA_dijkstra_shortest_paths(graph, travel_time, origin, init_node, term_node)
    no_node = nv(graph)
    no_arc = ne(graph)

    distmx = Inf*ones(no_node, no_node)
    for i in 1:no_arc
      distmx[init_node[i], term_node[i]] = travel_time[i]
    end

    state = dijkstra_shortest_paths(graph, origin, distmx)
    return state
end

function create_graph(init_node, term_node)
    @assert Base.length(init_node)==Base.length(term_node)

    no_node = max(maximum(init_node), maximum(term_node))
    no_arc = Base.length(init_node)

    graph = DiGraph(no_node)
    for i=1:no_arc
        add_edge!(graph, init_node[i], term_node[i])
    end
    return graph
end

function get_vector(state, origin, destination, link_dic)
    current = destination
    parent = -1
    x = zeros(Int, maximum(link_dic))

    while parent != origin && origin != destination && current != 0
        parent = state.parents[current]

        # println("origin=$origin, destination=$destination, parent=$parent, current=$current")

        if parent != 0
            link_idx = link_dic[parent,current]
            if link_idx != 0
                x[link_idx] = 1
            end
        end

        current = parent
    end

    return x
end

function add_demand_vector!(x, demand, state, origin, destination, link_dic)
  current = destination
  parent = -1

  while parent != origin && origin != destination && current != 0
      parent = state.parents[current]

      if parent != 0
          link_idx = link_dic[parent,current]
          if link_idx != 0
              x[link_idx] += demand
          end
      end

      current = parent
  end
end

# function get_shortest_path(init_node, term_node, link_length, origin, destination)
#     @assert Base.length(init_node)==Base.length(term_node)
#     @assert Base.length(init_node)==Base.length(link_length)
#
#     graph = create_graph(init_node, term_node)
#
#     state = dijkstra_shortest_paths(graph, link_length, origin)
#
#     path = get_path(state, origin, destination)
#     x = get_vector(path, init_node, term_node)
#
#     return path, x
# end
#
# function get_vector(state, origin, destination, init_node, term_node)
#     current = destination
#     parent = -1
#     x = zeros(Int, Base.length(init_node))
#
#     while parent != origin
#         parent = state.parents[current]
#
#         for j=1:Base.length(init_node)
#             if init_node[j]==parent && term_node[j]==current
#                 x[j] = 1
#                 break
#             end
#         end
#
#         current = parent
#     end
#
#     return x
# end


text:
### A Pluto.jl notebook ###
# v0.16.1

using Markdown
using InteractiveUtils

# ╔═╡ 50a16b1b-f79f-43ad-a5b4-71c99afba8ca
# julia stuff -
using PlutoUI

# ╔═╡ 26a6730a-21dd-11ec-3bbc-430e71b09495
md"""

#### Analysis of the Throttle Valve in the Vapor Compression Refrigeration Cycle (VCRC)

The throttle valve is used in the VCRC to convert a high-pressure saturated liquid stream to a vapor liquid mixture at lower pressure that is then fed into the evaporator. The throttle value is an _irreversible_ adiabatic operation.

Saturated T-508B enters into a throttle value at P$_{1}$ = 0.8 MPa (operating point $\mathcal{O}_{1}$) and exits at P$_{2}$ = 0.1 MPa (operating point $\mathcal{O}_{2}$). The mass flow rate of R-508B is $\dot{m}$ = 10 kg/s. The throttle value operates at steady-state.

Compute the following using enthalpy, entropy and R-508B pressure-enthalpy diagram:

1. The temperatures at $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$
1. The enthapy values at $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$, and the difference across the throttle valve
1. The entropy values at $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$, and the difference across the throttle value
1. Is there entropy generation between $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$? If so, can we compute $\dot{S}_{G}$?
1. What is the _quality_ of the outlet of the throttle valve?

"""

# ╔═╡ 3e2e945e-f5ee-4e28-87db-681e33726cb8
mdot = 10.0 	# units: kg/s

# ╔═╡ effcdabc-95ce-4a6c-8f47-b9b59c2d08b7
md"""
##### Pressure-Enthaplpy diagram for E-508B
"""

# ╔═╡ ed264b61-aca9-42b3-ab2e-0ba94fd26e30
PlutoUI.LocalResource("./figs/Fig-PH-508B-F1.png")

# ╔═╡ 52557b82-e273-478c-8572-115f90a24062
md"""

##### Lookup: The temperatures at $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$

We can lookup the temperatures at each operating point from the R-508B PH diagram.
"""

# ╔═╡ cb43d69c-abbc-45af-8775-3be12f0da20a
# Temperatures from the PH diagram
begin

	T₁ = -41.0 # units C (orange T)
	T₂ = -89.0 # units C (blue T)

	# return -
	(T₁,T₂)
end

# ╔═╡ ee496dc4-c9dc-46ee-bbae-49b41602626a
md"""

##### Compute: The _enthapy_ values at $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$, and the difference across the throttle valve

Throttle valves are _isoenthalpic_ i.e., there is no _enthaply_ change across the valve. To see this, let's start with the steady-state open energy balance:

$$\dot{Q}+\dot{W}_{sh}+\dot{W}_{b}+\sum_{s=1}^{\mathcal{S}}v_{s}\dot{m}_{s}H_{s} = 0$$

There is a single input, and a single output so we can show from the total mass balance that $\dot{m}_{1}=\dot{m}_{2}\equiv\dot{m}$. Next, a throttle valve is _adiabatic_, thus, $\dot{Q}$ = 0. Finally, there are no moving parts in a throttle valve and the system is at steady-state so there is no time-dependent boundary expansion or shart work terms i.e., $\dot{W}_{b}$ = 0 and $\dot{W}_{sh}$ = 0. Putting all these toget gives the energy balance:

$$\dot{m}\left(\sum_{s=1}^{\mathcal{S}}v_{s}H_{s}\right) = 0$$

Since $\dot{m}\neq{0}$, the summation term must be zero, or $H_{1}=H_{2}$.

"""

# ╔═╡ 67efc9e2-8edc-4551-b49d-c787c920230d
# Enthalpy from the PH-diagram
begin
	H₁ = 140.0 		# units: kJ/kg
	H₂ = 140.0 		# units: kJ/kg
	ΔH = H₂ - H₁ 	# units: kJ/kg

	# return -
	(H₁,H₂, ΔH)
end

# ╔═╡ 2990c0d2-b160-4b7b-9821-36cf744440a0
md"""

##### Compute: The _entropy_ values at $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$, and the difference across the throttle valve. Is there entropy generation?

Yes, there will be an entropy difference between $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$ (and entropy generation) because throttle values are _irreversible_. We can lookup the entropy values at $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$ and use the open second law (entropy balance) to compute the rate of entropy generation $\dot{S}_{G}$.

"""

# ╔═╡ d6960e86-cb59-4a9f-81cb-a3a54ae272f1
# Entropy values from the PH-diagram
begin
	S₁ = (0.75 + 0.8)/2 # units: kJ/kg-K
	S₂ = 0.8 			# units: kJ/kg-K
	ΔS = (S₂ - S₁)		# units: kJ/kg-G

	# return -
	(S₁,S₂, ΔS)
end

# ╔═╡ c301e469-9bca-475e-8e0d-101d0214ba0c
md"""
###### Compute the rate of entropy generation $\dot{S}_{G}$

Start from the open second law (entropy balance):

$$-\left(\sum_{s=1}^{\mathcal{S}}v_{s}\dot{m}_{s}S_{s}+\frac{\dot{Q}}{T_{\sigma}}\right) = \dot{S}_{G}$$

There is a single input, and a single output so we can show from the total mass balance that $\dot{m}_{1}=\dot{m}_{2}\equiv\dot{m}$.
Next, a throttle valve is _adiabatic_, thus, $\dot{Q}$ = 0. Finally, we have an expression for the rate of entropy generation:

$$-\dot{m}\left(\sum_{s=1}^{\mathcal{S}}v_{s}S_{s}\right) = \dot{S}_{G}$$

or (after running the summation out):

$$\dot{m}\left(S_{2} - S_{1}\right) = \dot{S}_{G}$$

"""

# ╔═╡ 5c72e8ac-24d0-4212-b3e8-b20a9c5c0171
# compute the rate of entropy generation -
S_dot_gen = mdot*(S₂ - S₁)	# units: kJ/s-K -or- W/K

# ╔═╡ d1de9861-f6a0-487c-9ef7-35fd3f1b6a68
md"""
#### Compute: what is the quality of the stream exiting the throttle valve?

The term _quality_ (denoted by $\theta$) in this context (or in the context of the Rankine cycle) refers to the fraction of vapor in the stream. To compute this, let's start with the expression for the entropy $S_\star$ which is valid for operating point $\star$ under the saturation dome (at some constant T):

$$S_{\star} = \theta_\star{S}_{v}+\left(1-\theta_\star\right)S_{l}$$

where $S_{v}$ denotes the entropy of a saturated vapor at T, $S_{l}$ denotes the entropy of a saturated liquid at T, and $\theta_\star$ denotes the quality of operating point $\star$ at T. We can lookup values for $S_{v}$ and $S_{l}$ as a function of T from the PH-diagram (or the associated saturation tables).

To comput $\theta_\star$ we can rearrange the expression above to give:

$$\theta_\star = \frac{S_{\star}- S_{l}}{S_{v} - S_{l}}$$

"""

# ╔═╡ 96bfd2aa-0ffc-4ed5-be20-1f3ca8bb88bb
begin

	# Lookup from the R-508B saturation table for T = -89C
	Sₗ = 0.5232 # units: kJ/kg-K
	Sᵥ = 1.4124 # units: kJ/kg-K

	# return -
	nothing
end

# ╔═╡ a80a7162-8aee-4a54-9199-bf2e5a0a0769
θ = (S₂ - Sₗ)/(Sᵥ - Sₗ)

# ╔═╡ 8f92b67f-555a-4797-b8fb-3bea29951f25
html"""<style>
main {
    max-width: 1200px;
  	width: 80%;
	margin: auto;
	font-family: "Roboto, monospace";
}

a {
	color: blue;
	text-decoration: none;
}

.H1 {
    padding: 0px 30px;
}
</style>
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.14"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[HypertextLiteral]]
git-tree-sha1 = "72053798e1be56026b81d4e2682dbe58922e5ec9"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.0"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "a8709b968a1ea6abc2dc1967cb1db6ac9a00dfb6"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.5"

[[PlutoUI]]
deps = ["Base64", "Dates", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "d1fb76655a95bf6ea4348d7197b22e889a4375f4"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.14"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╠═50a16b1b-f79f-43ad-a5b4-71c99afba8ca
# ╟─26a6730a-21dd-11ec-3bbc-430e71b09495
# ╠═3e2e945e-f5ee-4e28-87db-681e33726cb8
# ╟─effcdabc-95ce-4a6c-8f47-b9b59c2d08b7
# ╠═ed264b61-aca9-42b3-ab2e-0ba94fd26e30
# ╟─52557b82-e273-478c-8572-115f90a24062
# ╠═cb43d69c-abbc-45af-8775-3be12f0da20a
# ╟─ee496dc4-c9dc-46ee-bbae-49b41602626a
# ╠═67efc9e2-8edc-4551-b49d-c787c920230d
# ╟─2990c0d2-b160-4b7b-9821-36cf744440a0
# ╠═d6960e86-cb59-4a9f-81cb-a3a54ae272f1
# ╟─c301e469-9bca-475e-8e0d-101d0214ba0c
# ╠═5c72e8ac-24d0-4212-b3e8-b20a9c5c0171
# ╟─d1de9861-f6a0-487c-9ef7-35fd3f1b6a68
# ╠═96bfd2aa-0ffc-4ed5-be20-1f3ca8bb88bb
# ╠═a80a7162-8aee-4a54-9199-bf2e5a0a0769
# ╟─8f92b67f-555a-4797-b8fb-3bea29951f25
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
x = [0, 0, 1, 1, 1]
y = [0, 1, 0, 1, 1]
cm = IntConfusionMatrix(x, y)

@test typeof(cm) <: InfoTheory.ConfusionMatrix

@test sort(InfoTheory.xvalues(cm)) == [0, 1]
@test sort(InfoTheory.yvalues(cm)) == [0, 1]

xi0 = InfoTheory.xindex(cm, 0)
xi1 = InfoTheory.xindex(cm, 1)
@test in(xi0, 1:2)
@test in(xi1, 1:2)
@test xi0 != xi1

yi0 = InfoTheory.yindex(cm, 0)
yi1 = InfoTheory.yindex(cm, 1)
@test in(yi0, 1:2)
@test in(yi1, 1:2)
@test yi0 != yi1

@test count(cm, 0, 0) == 1
@test count(cm, 0, 1) == 1
@test count(cm, 1, 0) == 1
@test count(cm, 1, 1) == 2

x2 = [0, 0, 1, 1, -1]
y2 = [0, 1, 0, 1, 1]
cm2 = IntConfusionMatrix(x2, y2)
@test sort(InfoTheory.xvalues(cm2)) == [-1, 0, 1]
@test sort(InfoTheory.yvalues(cm2)) == [0, 1]

@test count(cm2, 0, 0) == 1
@test count(cm2, 0, 1) == 1
@test count(cm2, 1, 0) == 1
@test count(cm2, 1, 1) == 1
@test count(cm2, -1, 1) == 1
@test count(cm2, -1, 0) == 0


text:
# test_problems.jl : definition of some test problems

g(x, y) = Float64[1 + sin(π * xᵢ) * sin(π * yᵢ) for xᵢ in x, yᵢ in y]
g(x, y, z) = Float64[1 + sin(π * xᵢ) * sin(π * yᵢ) * sin(π * zᵢ) for xᵢ in x, yᵢ in y, zᵢ in z]

f(n, m) = g(range(1/n, stop = 1 - 1/n, length = n - 1), range(1/m, stop = 1 - 1/m, length = m - 1))
f(n, m, l) = g(range(1/n, stop = 1 - 1/n, length = n - 1),range(1/m, stop = 1 - 1/m, length = m - 1), range(1/l, stop = 1 - 1/l, length = l - 1))

p₁₂(n) = laplace2d(n, n)
p₁₃(n) = laplace3d(n, n, n)

p₂₂(n) = elliptic2d(f(n, n))
p₂₃(n) = elliptic3d(f(n, n, n))

p₃₂(n) = anisotropic2d(1e-4, n, n)
p₃₃(n) = anisotropic3d(1e-4,1, n, n, n)

p₄₂(n) = anisotropic2d(1e-4, 10 * π/180, n, n)


text:
# test stencil operator

include("grids.jl")
# include("operators.jl")

abstract type Stencil end

# Internal cell type
struct InternalCell <: Stencil end

# 6 faces of the domain
struct WestBoundaryCell <: Stencil end
struct EastBoundaryCell <: Stencil end
struct SouthBoundaryCell <: Stencil end
struct NorthBoundaryCell <: Stencil end
struct BottomBoundaryCell <: Stencil end
struct TopBoundaryCell <: Stencil end

# 12 sides of the domain
struct SouthWestBoundaryCell <: Stencil end
struct SouthEastBoundaryCell <: Stencil end
struct NorthWestBoundaryCell <: Stencil end
struct NorthEastBoundaryCell <: Stencil end
struct BottomWestBoundaryCell <: Stencil end
struct BottomEastBoundaryCell <: Stencil end
struct TopWestBoundaryCell <: Stencil end
struct TopEastBoundaryCell <: Stencil end
struct BottomSouthBoundaryCell <: Stencil end
struct BottomNorthBoundaryCell <: Stencil end
struct TopSouthBoundaryCell <: Stencil end
struct TopNorthBoundaryCell <: Stencil end

# 8 corners of the domain
struct BottomSouthWestBoundaryCell <: Stencil end
struct BottomSouthEastBoundaryCell <: Stencil end
struct BottomNorthWestBoundaryCell <: Stencil end
struct BottomNorthEastBoundaryCell <: Stencil end
struct TopSouthWestBoundaryCell <: Stencil end
struct TopSouthEastBoundaryCell <: Stencil end
struct TopNorthWestBoundaryCell <: Stencil end
struct TopNorthEastBoundaryCell <: Stencil end

"""Stencil for internal cell"""
function get_stencil!(cell, stencil_type::InternalCell, stencil, nb)
    (i, j, k) = cell
    nb = 6
    # stencil = Array{Tuple,1}(undef,6)
    stencil[1] = (i - 1, j, k)
    stencil[2] = (i + 1, j, k)
    stencil[3] = (i, j - 1, k)
    stencil[4] = (i, j + 1, k)
    stencil[5] = (i, j, k - 1)
    stencil[6] = (i, j, k + 1)
    return stencil, nb
end

# Stencil for west boundary cell
function get_stencil!(cell, stencil_type::WestBoundaryCell, stencil, nb)
    (i, j, k) = cell
    nb = 5
    # stencil = Array{Tuple,1}(undef,5)
    stencil[1] = (i + 1, j, k)
    stencil[2] = (i, j - 1, k)
    stencil[3] = (i, j + 1, k)
    stencil[4] = (i, j, k - 1)
    stencil[5] = (i, j, k + 1)
    return stencil, nb
end

# Stencil for east boundary cell
function get_stencil!(cell, stencil_type::EastBoundaryCell, stencil, nb)
    (i, j, k) = cell
    nb = 5
    # stencil = Array{Tuple,1}(undef,5)
    stencil[1] = (i - 1, j, k)
    stencil[2] = (i, j - 1, k)
    stencil[3] = (i, j + 1, k)
    stencil[4] = (i, j, k - 1)
    stencil[5] = (i, j, k + 1)
    return stencil, nb
end

# Stencil for south boundary cell
function get_stencil!(cell, stencil_type::SouthBoundaryCell, stencil, nb)
    (i, j, k) = cell
    nb = 5
    # stencil = Array{Tuple,1}(undef,5)
    stencil[1] = (i - 1, j, k)
    stencil[2] = (i + 1, j, k)
    stencil[3] = (i, j + 1, k)
    stencil[4] = (i, j, k - 1)
    stencil[5] = (i, j, k + 1)
    return stencil, nb
end

# Stencil for north boundary cell
function get_stencil!(cell, stencil_type::NorthBoundaryCell, stencil, nb)
    (i, j, k) = cell
    nb = 5
    # stencil = Array{Tuple,1}(undef,5)
stencil[1] = (i - 1, j, k)
    stencil[2] = (i + 1, j, k)
    stencil[3] = (i, j - 1, k)
    stencil[4] = (i, j, k - 1)
    stencil[5] = (i, j, k + 1)
    return stencil, nb
end

    # Stencil for bottom boundary cell
function get_stencil!(cell, stencil_type::BottomBoundaryCell, stencil, nb)
    (i, j, k) = cell
    nb = 5
    # stencil = Array{Tuple,1}(undef,5)
    stencil[1] = (i - 1, j, k)
    stencil[2] = (i + 1, j, k)
    stencil[3] = (i, j - 1, k)
    stencil[4] = (i, j + 1, k)
    stencil[5] = (i, j, k + 1)
    return stencil, nb
end

# Stencil for top boundary cell
function get_stencil!(cell, stencil_type::TopBoundaryCell, stencil, nb)
    (i, j, k) = cell
    nb = 5
    # stencil = Array{Tuple,1}(undef,5)
    stencil[1] = (i - 1, j, k)
    stencil[2] = (i + 1, j, k)
    stencil[3] = (i, j - 1, k)
    stencil[4] = (i, j + 1, k)
    stencil[5] = (i, j, k - 1)
    return stencil, nb
end

# Stencil for south-west boundary cell
function get_stencil!(cell, stencil_type::SouthWestBoundaryCell, stencil, nb)
(i, j, k) = cell
    nb = 4
    # stencil = Array{Tuple,1}(undef,4)
    stencil[1] = (i + 1, j, k)
    stencil[2] = (i, j + 1, k)
    stencil[3] = (i, j, k - 1)
    stencil[4] = (i, j, k + 1)
    return stencil, nb
end

# Stencil for south-west boundary cell
function get_stencil!(cell, stencil_type::SouthEastBoundaryCell, stencil, nb)
(i, j, k) = cell
    nb = 4
    # stencil = Array{Tuple,1}(undef,4)
    stencil[1] = (i - 1, j, k)
    stencil[2] = (i, j + 1, k)
    stencil[3] = (i, j, k - 1)
    stencil[4] = (i, j, k + 1)
    return stencil, nb
end

# Stencil for south-west boundary cell
function get_stencil!(cell, stencil_type::NorthWestBoundaryCell, stencil, nb)
(i, j, k) = cell
    nb = 4
    # stencil = Array{Tuple,1}(undef,4)
    stencil[1] = (i + 1, j, k)
    stencil[2] = (i, j - 1, k)
    stencil[3] = (i, j, k - 1)
    stencil[4] = (i, j, k + 1)
    return stencil, nb
end

# Stencil for south-west boundary cell
function get_stencil!(cell, stencil_type::NorthEastBoundaryCell, stencil, nb)
(i, j, k) = cell
    nb = 4
    # stencil = Array{Tuple,1}(undef,4)
    stencil[1] = (i - 1, j, k)
    stencil[2] = (i, j - 1, k)
    stencil[3] = (i, j, k - 1)
    stencil[4] = (i, j, k + 1)
    return stencil, nb
end

# Stencil for bottom-west boundary cell
    function get_stencil!(cell, stencil_type::BottomWestBoundaryCell, stencil, nb)
(i, j, k) = cell
    nb = 4
    # stencil = Array{Tuple,1}(undef,4)
    stencil[1] = (i + 1, j, k)
    stencil[2] = (i, j - 1, k)
    stencil[3] = (i, j + 1, k)
    stencil[4] = (i, j, k + 1)
    return stencil, nb
end

# Stencil for bottom-east boundary cell
    function get_stencil!(cell, stencil_type::BottomEastBoundaryCell, stencil, nb)
(i, j, k) = cell
    nb = 4
    # stencil = Array{Tuple,1}(undef,4)
    stencil[1] = (i - 1, j, k)
    stencil[2] = (i, j - 1, k)
    stencil[3] = (i, j + 1, k)
    stencil[4] = (i, j, k + 1)
    return stencil, nb
end

# Stencil for bottom-south boundary cell
function get_stencil!(cell, stencil_type::BottomSouthBoundaryCell, stencil, nb)
(i, j, k) = cell
    nb = 4
    # stencil = Array{Tuple,1}(undef,4)
    stencil[1] = (i - 1, j, k)
    stencil[2] = (i + 1, j, k)
    stencil[3] = (i, j + 1, k)
    stencil[4] = (i, j, k + 1)
    return stencil, nb
end

# Stencil for bottom-north boundary cell
function get_stencil!(cell, stencil_type::BottomNorthBoundaryCell, stencil, nb)
(i, j, k) = cell
    nb = 4
    # stencil = Array{Tuple,1}(undef,4)
    stencil[1] = (i - 1, j, k)
    stencil[2] = (i + 1, j, k)
    stencil[3] = (i, j - 1, k)
    stencil[4] = (i, j, k + 1)
    return stencil, nb
end

# Stencil for top-west boundary cell
    function get_stencil!(cell, stencil_type::TopWestBoundaryCell, stencil, nb)
(i, j, k) = cell
    nb = 4
    # stencil = Array{Tuple,1}(undef,4)
    stencil[1] = (i + 1, j, k)
    stencil[2] = (i, j - 1, k)
    stencil[3] = (i, j + 1, k)
    stencil[4] = (i, j, k - 1)
    return stencil, nb
end

# Stencil for top-east boundary cell
    function get_stencil!(cell, stencil_type::TopEastBoundaryCell, stencil, nb)
(i, j, k) = cell
    nb = 4
    # stencil = Array{Tuple,1}(undef,4)
    stencil[1] = (i - 1, j, k)
    stencil[2] = (i, j - 1, k)
    stencil[3] = (i, j + 1, k)
    stencil[4] = (i, j, k - 1)
    return stencil, nb
end

# Stencil for top-south boundary cell
function get_stencil!(cell, stencil_type::TopSouthBoundaryCell, stencil, nb)
(i, j, k) = cell
    nb = 4
    # stencil = Array{Tuple,1}(undef,4)
    stencil[1] = (i - 1, j, k)
    stencil[2] = (i + 1, j, k)
    stencil[3] = (i, j + 1, k)
    stencil[4] = (i, j, k - 1)
    return stencil, nb
end

# Stencil for top-north boundary cell
function get_stencil!(cell, stencil_type::TopNorthBoundaryCell, stencil, nb)
(i, j, k) = cell
    nb = 4
    # stencil = Array{Tuple,1}(undef,4)
    stencil[1] = (i - 1, j, k)
    stencil[2] = (i + 1, j, k)
    stencil[3] = (i, j - 1, k)
    stencil[4] = (i, j, k - 1)
    return stencil, nb
end

# Stencil for bottom south-west boundary cell
function get_stencil!(cell, stencil_type::BottomSouthWestBoundaryCell, stencil, nb)
    (i, j, k) = cell
nb = 3
    # stencil = Array{Tuple,1}(undef,3)
    stencil[1] = (i + 1, j, k)
    stencil[2] = (i, j + 1, k)
    stencil[3] = (i, j, k + 1)
    return stencil, nb
end

# Stencil for bottom south-east boundary cell
function get_stencil!(cell, stencil_type::BottomSouthEastBoundaryCell, stencil, nb)
    (i, j, k) = cell
nb = 3
    # stencil = Array{Tuple,1}(undef,3)
    stencil[1] = (i - 1, j, k)
    stencil[2] = (i, j + 1, k)
    stencil[3] = (i, j, k + 1)
    return stencil, nb
end

# Stencil for bottom north-west boundary cell
function get_stencil!(cell, stencil_type::BottomNorthWestBoundaryCell, stencil, nb)
    (i, j, k) = cell
nb = 3
    # stencil = Array{Tuple,1}(undef,3)
    stencil[1] = (i + 1, j, k)
    stencil[2] = (i, j - 1, k)
    stencil[3] = (i, j, k + 1)
    return stencil, nb
end

# Stencil for bottom north-east boundary cell
function get_stencil!(cell, stencil_type::BottomNorthEastBoundaryCell, stencil, nb)
    (i, j, k) = cell
    nb = 3
    # stencil = Array{Tuple,1}(undef,3)
    stencil[1] = (i - 1, j, k)
    stencil[2] = (i, j - 1, k)
    stencil[3] = (i, j, k + 1)
    return stencil, nb
end

# Stencil for top south-west boundary cell
function get_stencil!(cell, stencil_type::TopSouthWestBoundaryCell, stencil, nb)
    (i, j, k) = cell
    nb = 3
    # stencil = Array{Tuple,1}(undef,3)
    stencil[1] = (i + 1, j, k)
    stencil[2] = (i, j + 1, k)
    stencil[3] = (i, j, k - 1)
    return stencil, nb
end

# Stencil for top south-east boundary cell
function get_stencil!(cell, stencil_type::TopSouthEastBoundaryCell, stencil, nb)
    (i, j, k) = cell
    nb = 3
    # stencil = Array{Tuple,1}(undef,3)
    stencil[1] = (i - 1, j, k)
    stencil[2] = (i, j + 1, k)
    stencil[3] = (i, j, k - 1)
    return stencil, nb
end

# Stencil for top north-west boundary cell
function get_stencil!(cell, stencil_type::TopNorthWestBoundaryCell, stencil, nb)
    (i, j, k) = cell
    nb = 3
    # stencil = Array{Tuple,1}(undef,3)
    stencil[1] = (i + 1, j, k)
    stencil[2] = (i, j - 1, k)
    stencil[3] = (i, j, k - 1)
    return stencil, nb
end

# Stencil for top north-east boundary cell
function get_stencil!(cell, stencil_type::TopNorthEastBoundaryCell, stencil, nb)
    (i, j, k) = cell
    nb = 3
    # stencil = Array{Tuple,1}(undef,3)
    stencil[1] = (i - 1, j, k)
    stencil[2] = (i, j - 1, k)
    stencil[3] = (i, j, k - 1)
    return stencil, nb
end

function set_stencil_type(grid) # , stencil_type)

    nx = grid.nx
    ny = grid.ny
    nz = grid.nz

    stencil_type = Array{Any,3}(undef, (nx, ny, nz))

    for k = 2:nz
        for j = 2:ny
            for i = 1:nx
                stencil_type[i,j,k] = InternalCell()
            end
        end
    end

    # Initialize all cells to be of Internal cell type
    # stencil_type[: ,: ,: ] .= InternalCell

    # Modify the 6 faces of the domain
    for k = 2:nz
        for j = 2:ny
            stencil_type[1 ,j ,k ] = WestBoundaryCell()
            stencil_type[nx,j ,k ] = EastBoundaryCell()
        end
    end

    for k = 2:nz
        for i = 1:nx
            stencil_type[i ,1 ,k ] = SouthBoundaryCell()
            stencil_type[i ,ny,k ] = NorthBoundaryCell()
        end
    end

    for j = 2:ny
        for i = 1:nx
        stencil_type[i ,j ,1 ] = BottomBoundaryCell()
            stencil_type[i ,j ,nz] = TopBoundaryCell()
        end
    end

    # Modify the 12 sides of the domain
    for k = 1:nz
        stencil_type[1 ,1 ,k ] = SouthWestBoundaryCell()
        stencil_type[nx,1 ,k ] = SouthEastBoundaryCell()
        stencil_type[1 ,ny,k ] = NorthWestBoundaryCell()
        stencil_type[nx,ny,k ] = NorthEastBoundaryCell()
    end

    for j = 1:ny
        stencil_type[1 ,j ,1 ] = BottomWestBoundaryCell()
        stencil_type[nx,j ,1 ] = BottomEastBoundaryCell()
        stencil_type[1 ,j ,nz] = TopWestBoundaryCell()
        stencil_type[nx,j ,nz] = TopEastBoundaryCell()
    end

    for i = 1:nx
        stencil_type[i ,1 ,1 ] = BottomSouthBoundaryCell()
        stencil_type[i ,ny,1 ] = BottomNorthBoundaryCell()
        stencil_type[i ,1 ,nz] = TopSouthBoundaryCell()
        stencil_type[i ,ny,nz] = TopNorthBoundaryCell()
    end

    # Finally modify the 4 domain corners
    stencil_type[1 ,1 ,1 ] = BottomSouthWestBoundaryCell()
    stencil_type[nx,1 ,1 ] = BottomSouthEastBoundaryCell()
    stencil_type[1 ,ny,1 ] = BottomNorthWestBoundaryCell()
    stencil_type[nx,ny,1 ] = BottomNorthEastBoundaryCell()
    stencil_type[1 ,1 ,nz] = TopSouthWestBoundaryCell()
    stencil_type[nx,1 ,nz] = TopSouthEastBoundaryCell()
    stencil_type[1 ,ny,nz] = TopNorthWestBoundaryCell()
    stencil_type[nx,ny,nz] = TopNorthEastBoundaryCell()

    return stencil_type
end

@inline function kernel!(cell, h, hnew, K, stencil, nb)
    (i, j, k) = cell
    # @show nb
    hnew[i,j,k] = h[i,j,k]
    for ib = 1:nb
        hnew[i,j,k] += add_flux(cell, h, K, stencil[ib])
    end
    return hnew
end

function kernel_no_stencil!(cell, h, hnew, K)
    (i, j, k) = cell

    hnew[i,j,k] = h[i,j,k] +
K[i,j,k] * (h[i + 1,j  ,k  ] - 2h[i,j,k] + h[i - 1,j  ,k  ] +
                              h[i  ,j + 1,k  ] - 2h[i,j,k] + h[i  ,j - 1,k  ] +
    h[i  ,j  ,k + 1] - 2h[i,j,k] + h[i  ,j  ,k - 1] )
    return hnew
end

function add_flux(cell, h, K, neighb)
    (i, j, k) = cell
    (in, jn, kn) = neighb
    # Take the direction of the neighbour into account
    s = Float64(sign(i - in) + sign(j - jn) + sign(k - kn))
    flux = s * K[i,j,k] * (h[i,j,k] - h[in,jn,kn])
end

function test_full_stencil_based!(grid, h, hnew, K, stencil_type)
                nx = grid.nx
    ny = grid.ny
    nz = grid.nz

    nb      = 6
    stencil = Array{Tuple,1}(undef, 6)

    @inbounds for k = 1:nz
        @inbounds for j = 1:ny
            @inbounds for i = 1:nx
                cell = (i, j, k)
                stencil, nb = get_stencil!(cell, stencil_type[i,j,k], stencil, nb)
                hnew = kernel!(cell, h, hnew, K, stencil, nb)
                # @show size(neighbours)
            end
        end
    end

end

function test_stencil_based!(grid, h, hnew, K, stencil_type)
    nx = grid.nx
    ny = grid.ny
    nz = grid.nz

    nb      = 6
    stencil = Array{Tuple,1}(undef, 6)

    @inbounds for k = 2:nz - 1
        @inbounds for j = 2:ny - 1
            @inbounds for i = 2:nx - 1
                cell = (i, j, k)
                stencil, nb = get_stencil!(cell, stencil_type[i,j,k], stencil, nb)
                hnew = kernel!(cell, h, hnew, K, stencil, nb)
            end
        end
    end

end

function test_no_stencil!(grid, h, hnew, K)
    nx = grid.nx
    ny = grid.ny
    nz = grid.nz

    @inbounds for k = 2:nz - 1
        @inbounds for j = 2:ny - 1
            @inbounds for i = 2:nx - 1
                cell = (i, j, k)
                hnew = kernel_no_stencil!(cell, h, hnew, K)
            end
        end
    end

end

function test_stencil()

    nx = 100
    ny = 80
    nz = 50

    @show total = nx * ny * nz
    @show internal = (nx - 1) * (ny - 1) * (nz - 1)
    @show boundary = total - internal
    @show boundary_β = 100 * (boundary / total)

    dx = 10.0
    dy = 10.0
    dz = 10.0
    useCUDA = false

    h0 = 10.0
    K0 = 1.0e-4

    h = fill(h0, (nx, ny, nz))
    hnew = copy(h)
    K = fill(K0, (nx, ny, nz))

    x, y, z = grid_coords(nx, ny, nz, dx, dy, dz, useCUDA)

    grid = Grid(nx, ny, nz, dx, dy, dz, x, y, z)

    # stencil_type = fill(nothing, (nx, ny, nz))

    stencil_type = set_stencil_type(grid) # , stencil_type)

    @time test_full_stencil_based!(grid, h, hnew, K, stencil_type)
    @time test_stencil_based!(grid, h, hnew, K, stencil_type)
    @time test_no_stencil!(grid, h, hnew, K)

end

test_stencil()


text:
#Create and initialize shaders
module ShadersAndVerticiesForText
using ModernGL, GeometryTypes, GLFW, ..ForDisplayStructs,  ..CustomFragShad,  ..ModernGlUtil

export getWordsVerticies

"""
creating VertexShader  so controlling structures like verticies, quads
gslString so version of GSLS we are using currently
  """
function createVertexShader(gslString::String)
vsh = """
$(gslString)
layout (location = 0) in vec3 aPos;
layout (location = 1) in vec3 aColor;
layout (location = 2) in vec2 aTexCoord;
out vec3 ourColor;
smooth out vec2 TexCoord0;
void main()
{
    gl_Position = vec4(aPos, 1.0);
    ourColor = aColor;
  // TexCoord0 = vec2(aTexCoord.y, aTexCoord.x);
   TexCoord0 = aTexCoord;

}
"""
return createShader(vsh, GL_VERTEX_SHADER)
end

"""
creating fragment Shader  so controlling colors and textures
gslString so version of GSLS we are using currently
  """
function createFragmentShader(gslString::String)
    fsh = """
    $(gslString)

    #extension GL_EXT_gpu_shader4 : enable    //Include support for this extension, which defines usampler2D

    out vec4 FragColor;
    in vec3 ourColor;
    smooth in vec2 TexCoord0;

    uniform usampler2D TextTexture1;
    void main() {

    uint text1Texel = texture2D(TextTexture1, TexCoord0).r ;

     if(text1Texel > 0){
      FragColor = vec4(0.0,0.0,0.0,1.0);  }
       else {
    FragColor = vec4(1.0,1.0,1.0,1.0);

    }
    }

    """
    return createShader(fsh, GL_FRAGMENT_SHADER)
    end

################### data to display verticies

# Specify how vertices are arranged into faces
# Face{N,T} type specifies a face with N vertices, with index type
# T (you should choose UInt32), and index-offset O. If you're
# specifying faces in terms of julia's 1-based indexing, you should set
# O=0. (If you instead number the vertices starting with 0, set
# O=-1.)
elements = Face{3,UInt32}[(0,1,2),          # the first triangle
(2,3,0)]          # the second triangle

"""
generetes verticies for quad used for displaying text
"""
function getWordsVerticies(fractionOfMainIm::Float32)::Vector{Float32}
  correctedWidthForTextAccounting = (-1+ fractionOfMainIm*2)

  return  Float32.([
  # positions                     // colors           // texture coords
   1.0,  1.0, 0.0,                 1.0, 0.0, 0.0,   1.0, 1.0,   # top right
   1.0, -1.0, 0.0,                0.0, 1.0, 0.0,   1.0, 0.0,   # bottom right
   correctedWidthForTextAccounting, -1.0, 0.0,   0.0, 0.0, 1.0,   0.0, 0.0,   # bottom left
   correctedWidthForTextAccounting,  1.0, 0.0,   1.0, 1.0, 0.0,   0.0, 1.0    # top left
   ])
end #getWordsVerticies

  end #..ShadersAndVerticies


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    triangles

function triangles(graph::UndirectedGraph)
    triads = Set()
    for n in nodeset(graph)
        nkeys = adjacencies(graph, n)
        length(nkeys) < 2 && continue
        for (u, v) in combinations(nkeys)
            if hasedge(graph, u, v)
                push!(triads, tuple(sort([n, u, v])...))
            end
        end
    end
    return triads
end


text:
#
# genmult.jl -
#
# Generalized dot product by grouping consecutive dimensions.
#
#-------------------------------------------------------------------------------
#
# This file is part of LazyAlgebra (https://github.com/emmt/LazyAlgebra.jl)
# released under the MIT "Expat" license.
#
# Copyright (c) 2017-2020 Éric Thiébaut.
#

module GenMult

export
    lgemm!,
    lgemm,
    lgemv!,
    lgemv

using ..LazyAlgebra
using ..LazyAlgebra: Complexes, Floats, Reals, axes, promote_multiplier,
    libblas, @blasfunc, BlasInt, BlasReal, BlasFloat, BlasComplex,
    bad_argument, bad_size
using ArrayTools  # for `cartesian_indices`, `is_flat_array`, etc.
using LinearAlgebra
using LinearAlgebra.BLAS

"""

```julia
Implementation(Val(:alg), args...)`
```

is used to quickly determine the most efficient implementation of the code to
use for algorithm `alg` with arguments `args...`.  The returned value is one of
four possible singletons:

- `Blas()` when highly optimized BLAS code can be used.  This is the preferred
  implementation as it is assumed to be the fastest.

- `Basic()` when *vector* and *matrix* arguments have respectively one and two
  dimensions.

- `Linear()` when *vector* and *matrix* arguments can be efficiently indexed
  by, respectively, one and two linear indices.

- `Generic()` to use generic implementation which can accommodate from any type
  of arguments and of multi-dimensional indices.  This implementation should be
  always safe to use and should provide the reference implementation of the
  algorithm `alg`.

Whenever possible, the best implementation is automatically determined at
compilation time by calling this method.

"""
abstract type Implementation end

for S in (:Basic, :Blas, :Linear, :Generic)
    @eval begin
        struct $S <: Implementation end
        @doc @doc(Implementation) $S
    end
end

incompatible_dimensions() =
    bad_size("incompatible dimensions")

invalid_transpose_character() =
    bad_argument("invalid transpose character")

include("lgemv.jl")
include("lgemm.jl")

end # module


text:
using AstroLib: trueanom, kepler_solver
using Unitful, UnitfulAstro
using Unitful: AbstractQuantity
using KeywordCalls
using Rotations

# Domain specific unit conversions / constants/ fallbacks
Unitful.preferunits(u"Msun,Rsun,d"...)
const G_unit = Unitful.G
const G_nom = ustrip(u"Rsun^3/Msun/d^2", G_unit)

# Helpers
include("kepler_helpers.jl")

"""
    KeplerianOrbit(; kwargs...)

Keplerian orbit parameterized by the basic observables of a transiting 2-body system.

# Parameters
* `period`/`P` -- The orbital period of the planet [d].
* `t0`/`t_0` -- The midpoint time of the reference transit [d].
* `tp`/`t_p` -- The time of periastron [d].
* `duration`/`τ`/`T` -- The transit duration [d].
* `a` -- The semi-major axis [R⊙].
* `aR_star`/`aRs` -- The ratio of the semi-major axis to star radius.
* `R_planet`/`Rp` -- The radius of the planet [R⊙].
* `R_star`/`Rs` -- The radius of the star [R⊙].
* `rho_star`/`ρ_star` -- The spherical star density [M⊙/R⊙³].
* `r`/`RpRs` -- The ratio of the planet radius to star radius.
* `b` -- The impact parameter, bounded between 0 ≤ b ≤ 1.
* `ecc`/`e` -- The eccentricity of the closed orbit, bounded between 0 ≤ ecc < 1.
* `incl` -- The inclination of the orbital plane relative to the axis perpendicular to the
           reference plane [rad]
* `omega`/`ω` -- The argument of periapsis [rad].
* `cos_omega`/`cos_ω` -- The cosine of the argument of periapsis.
* `sin_omega`/`sin_ω` -- The sine of the argument of periapsis.
* `Omega`/`Ω` -- The longitude of the ascending node [rad].
* `M_planet`/`Mp` -- The mass of the planet [M⊙].
* `M_star`/`Ms` -- The mass of the star [M⊙].

# Valid combinations
The following flowchart can be used to determine which parameters can define a `KeplerianOrbit`:
1. The `period` or `a` must be given. If both given, then neither `M_star` or `rho_star` can be defined because the stellar density is now implied.
2. Only `incl` or `b` can be given.
3. If `ecc` is given, then `omega` must also be given.
4. If no stellar parameters are given, the central body is assumed to be the Sun. If only `rho_star` is given, then `R_star` is defined to be 1 solar radius. Otherwise, at most two of `M_star`, `R_star`, and `rho_star` can be given.
5. Either `t0` or `tp` must be given, but not both.
"""
@concrete struct KeplerianOrbit <: AbstractOrbit
    period
    t0
    tp
    t_ref
    duration
    a
    a_planet
    a_star
    R_planet
    R_star
    rho_planet
    rho_star
    r
    aR_star
    b
    ecc
    M0
    cos_incl
    sin_incl
    cos_omega
    sin_omega
    cos_Omega
    sin_Omega
    incl
    omega
    Omega
    n
    M_planet
    M_star
end

function KeplerianOrbit(nt::NamedTuple{(
        :period, :t0, :tp, :duration,
        :a, :R_planet, :R_star,
        :rho_star,
        :r, :aR_star, :b, :ecc, :cos_omega, :sin_omega,
        :incl, :omega, :Omega,
        :M_planet, :M_star,
    )})
    if nt.period isa Real || nt.a isa Real
        G = G_nom
    else
        G = G_unit
    end

    if (isnothing(nt.ecc) || iszero(nt.ecc)) && !isnothing(nt.duration)
        isnothing(nt.b) && throw(ArgumentError(
            "`b` must also be provided for a circular orbit if `duration given`"
        ))
        isnothing(nt.r) && throw(ArgumentError(
            "`r` must also be provided if `duration` given"
        ))
    end

    a, aR_star, period, rho_star, R_star, M_star, M_planet, duration = compute_consistent_inputs(
        nt.a, nt.aR_star, nt.period, nt.rho_star, nt.R_star, nt.M_star, nt.M_planet, G,
        nt.ecc, nt.duration, nt.b, nt.r,
    )
    r = isnothing(nt.r) ? zero(aR_star) : nt.r
    M_tot = compute_M_tot(M_star, M_planet)
    R_planet = compute_R_planet(R_star, r, nt.R_planet)
    rho_planet = compute_rho(M_planet, R_planet)

    n = 2.0 * π / period
    a_star = compute_a_X(a, M_planet, M_tot)
    a_planet = -compute_a_X(a, M_star, M_tot)

    # Omega
    Omega = nt.Omega
    sin_Omega, cos_Omega = sincos(Omega)

    # Eccentricity
    ecc = nt.ecc
    if iszero(ecc)
        M0 = 0.5 * π
        incl_factor_inv = 1.0
        omega = zero(Omega)
        cos_omega, sin_omega = 1.0, 0.0
    else
        if !isnothing(nt.omega)
            all(!isnothing, (nt.cos_omega, nt.sin_omega)) && throw(ArgumentError(
                "Only `ω`, or `cos_ω` and `sin_ω` can be provided"
            ))
            omega = nt.omega
            sin_omega, cos_omega = sincos(nt.omega)
        else
            cos_omega, sin_omega = nt.cos_omega, nt.sin_omega
            omega = atan(sin_omega, cos_omega)
        end
        E0 = compute_E0(ecc, cos_omega, sin_omega)
        M0 = compute_M0(ecc, E0)

        incl_factor_inv = compute_incl_factor_inv(ecc, sin_omega)
    end

    # Jacobian for cos(i) -> b
    dcosi_db = compute_dcosi_db(a, R_star, incl_factor_inv)

    if !isnothing(nt.b)
        if any(!isnothing, (nt.incl, duration))
            throw(ArgumentError("Only `incl`, `b`, or `duration` can be given"))
        end
        b = nt.b
        cos_incl = dcosi_db * b
        incl = acos(cos_incl)
        sin_incl = sin(incl)
        duration = nt.duration
    elseif !isnothing(nt.incl)
        !isnothing(nt.duration) && throw(ArgumentError(
            "Only `incl`, `b`, or `duration` can be given"
        ))
        incl = nt.incl
        sin_incl, cos_incl = sincos(incl)
        b = compute_b(cos_incl, dcosi_db)
        duration = nt.duration
    elseif !isnothing(nt.duration)
        duration = nt.duration
        b = compute_b(a_planet, R_star, duration, period, incl_factor_inv, ecc, sin_omega)
        cos_incl = dcosi_db * b
        incl = acos(cos_incl)
        sin_incl = sin(incl)
    else
        incl = 0.5 * π
        cos_incl = 0.0
        sin_incl = 1.0
        b = 0.0
        duration = nt.duration
    end

    # Compute remaining system parameters
    t0, tp = compute_t0_tp(nt.t0, nt.tp; M0=M0, n=n)
    t_ref = tp - t0

    # Sanitize dimensionless units
    if period isa AbstractQuantity
        r, aR_star, b = NoUnits.((r, aR_star, b))
    end

    return KeplerianOrbit(
        period, t0, tp, t_ref, duration,
        a, a_planet, a_star, R_planet, R_star,
        rho_planet, rho_star,
        r, aR_star, b, ecc, M0, cos_incl, sin_incl, cos_omega, sin_omega, cos_Omega, sin_Omega,
        incl, omega, Omega,
        n,
        M_planet, M_star,
    )
end

@kwcall KeplerianOrbit(
    period=nothing, t0=nothing, tp=nothing, duration=nothing,
    a=nothing, R_planet=nothing, R_star=nothing,
    rho_star=nothing,
    r=nothing, aR_star=nothing, b=nothing, ecc=0.0, cos_omega=nothing, sin_omega=nothing,
    incl=nothing, omega=0.0, Omega=0.0,
    M_planet=nothing, M_star=nothing,
)

@kwalias KeplerianOrbit [
    P => period,
    t_0 => t0,
    t_p => tp,
    τ => duration,
    T => duration,
    aRs => aR_star,
    Rp => R_planet,
    Rs => R_star,
    ρ_star => rho_star,
    RpRs => r,
    e => ecc,
    cos_ω => cos_omega,
    sin_ω => sin_omega,
    ω => omega,
    Ω => Omega,
    Mp => M_planet,
    Ms => M_star,
]

# Finds the position `r` of the planet along its orbit after rotating
# through the true anomaly `ν`, then transforms this from the
# orbital plan to the equatorial plane
# separation: aR_star, a/R_star, a_star / R_star, or a_planet / R_star
# TODO: consider moving this to a separate orbital calculations package in the future
function _position(orbit, separation, t)
    sin_ν, cos_ν = compute_true_anomaly(orbit, t)
    if isnothing(orbit.ecc) || iszero(orbit.ecc)
        r = separation
    else
        r = separation * (1 - orbit.ecc^2) / (1 + orbit.ecc * cos_ν)
    end
    # Transform from orbital plane to equatorial plane
    X = SA[r * cos_ν, r * sin_ν, zero(r)]
    R = RotZXZ(orbit.Omega, -orbit.incl, orbit.omega)
    return R * X
end
_star_position(orb, R_star, t) = _position.(orb, orb.a_star / R_star, t)
_planet_position(orb, R_star, t) = _position.(orb, orb.a_planet / R_star, t)
relative_position(orbit::KeplerianOrbit, t) = _position(orbit, -orbit.aR_star, t)

# Returns sin(ν), cos(ν)
function compute_true_anomaly(orbit::KeplerianOrbit, t)
    M = compute_M(t, orbit.t0, orbit.t_ref, orbit.n)
    if isnothing(orbit.ecc) || iszero(orbit.ecc)
        return sincos(M)
    else
        E = kepler_solver(M, orbit.ecc)
        return sincos(trueanom(E, orbit.ecc))
    end
end

function flip(orbit::KeplerianOrbit, R_planet)
    M_planet = orbit.M_planet
    if isnothing(orbit.ecc) || iszero(orbit.ecc)
        return KeplerianOrbit(
            period = orbit.period,
            tp = orbit.tp + 0.5*orbit.period,
            incl = orbit.incl,
            Omega = orbit.Omega,
            omega = orbit.omega,
            M_star = orbit.M_planet,
            M_planet = orbit.M_star,
            R_star = R_planet,
            R_planet = orbit.R_star,
            ecc = orbit.ecc,
        )
    else
        return KeplerianOrbit(
            period = orbit.period,
            tp = orbit.tp,
            incl = orbit.incl,
            Omega = orbit.Omega,
            omega = orbit.omega - π,
            M_star = orbit.M_planet,
            M_planet = orbit.M_star,
            R_star = R_planet,
            R_planet = orbit.R_star,
            ecc = orbit.ecc,
        )
    end
end

function compute_consistent_inputs(a, aR_star, period, rho_star, R_star, M_star, M_planet, G, ecc, duration, b, r)
    if isnothing(a) && isnothing(period)
        throw(ArgumentError("At least `a` or `P` must be specified"))
    end

    if (isnothing(ecc) || iszero(ecc)) && !isnothing(duration)
        isnothing(R_star) && (R_star = compute_R_star_nom(G))
        aR_star = compute_aor(duration, period, b, r=r)
        a = compute_a(aR_star, R_star)
        duration = nothing
    end

    if isnothing(M_planet) && (!isnothing(a) || !isnothing(period))
        M_planet = compute_M_planet_nom(G)
    end

    # Compute implied stellar density
    implied_rho_star = false
    if all(!isnothing, (a, period))
        if any(!isnothing, (rho_star, M_star))
            throw(ArgumentError(
                "If both `a` and `P` are given, `rho_star` or `M_star` cannot be defined"
            ))
        end

        # Default to R_star = 1 R⊙ if not provided
        isnothing(R_star) && (R_star = oneunit(a))

        # Compute implied mass
        M_tot = compute_M_tot(a, G, period)

        # Compute implied density
        M_star = M_tot - M_planet
        rho_star = compute_rho_star(M_star, R_star)
        implied_rho_star = true
    end

    # Check combination of stellar params are valid
    if all(isnothing, (R_star, M_star))
        R_star = compute_R_star_nom(G)
        isnothing(rho_star) && (M_star = oneunit(M_planet))
    end

    if !implied_rho_star && sum(isnothing, (rho_star, R_star, M_star)) ≠ 1
        throw(ArgumentError(
            "Must provide exactly two of: `rho_star`, `R_star`, or `M_star` if rho_star not implied"
        ))
    end

    # Compute stellar params
    if isnothing(rho_star)
        rho_star = compute_rho_star(M_star, R_star)
    elseif isnothing(R_star)
        R_star = compute_R_star(rho_star, M_star)
    else
        M_star = compute_M_star(rho_star, R_star)
    end

    # Compute planet params
    M_tot = compute_M_tot(M_star, M_planet)
    if isnothing(a)
        if isnothing(aR_star)
            a = compute_a(M_tot, period, G)
            aR_star = a / R_star
        else
            a = compute_a(aR_star, R_star)
        end
    else
        period = compute_period(M_tot, a, G)
        aR_star = a / R_star
    end

    return a, aR_star, period, rho_star, R_star, M_star, M_planet, duration
end

stringify_units(value, unit_str) = @sprintf "%.4f %s" value unit_str
function stringify_units(value::Unitful.AbstractQuantity, unit_str)
    u = upreferred(value)
    return stringify_units(ustrip(u), string(unit(u)))
end
stringify_units(value::Nothing, unit) = "$(value)"
stringify_units(value) = @sprintf "%.4f"  value
function Base.show(io::IO, ::MIME"text/plain", orbit::KeplerianOrbit)
    print(
        IOContext(io, :fancy_exponent => false),
        """
        Keplerian Orbit
         P:      $(stringify_units(orbit.period, "d"))
         t₀:     $(stringify_units(orbit.t0, "d"))
         tₚ:     $(stringify_units(orbit.tp, "d"))
         t_ref:  $(stringify_units(orbit.t_ref, "d"))
         τ:      $(stringify_units(orbit.duration, "d"))
         a:      $(stringify_units(orbit.a, "R⊙"))
         aₚ:     $(stringify_units(orbit.a_planet, "R⊙"))
         aₛ:     $(stringify_units(orbit.a_star, "R⊙"))
         Rₚ:     $(stringify_units(orbit.R_planet, "R⊙"))
         Rₛ:     $(stringify_units(orbit.R_star, "R⊙"))
         ρₚ:     $(stringify_units(orbit.rho_planet, "M⊙/R⊙³"))
         ρₛ:     $(stringify_units(orbit.rho_star, "M⊙/R⊙³"))
         r:      $(stringify_units(orbit.r))
         aRₛ:    $(stringify_units(orbit.aR_star))
         b:      $(stringify_units(orbit.b))
         ecc:    $(stringify_units(orbit.ecc))
         cos(i): $(stringify_units(orbit.cos_incl))
         sin(i): $(stringify_units(orbit.sin_incl))
         cos(ω): $(stringify_units(orbit.cos_omega))
         sin(ω): $(stringify_units(orbit.sin_omega))
         cos(Ω): $(stringify_units(orbit.cos_Omega))
         sin(Ω): $(stringify_units(orbit.sin_Omega))
         i:      $(stringify_units(orbit.incl, "rad"))
         ω:      $(stringify_units(orbit.omega, "rad"))
         Ω:      $(stringify_units(orbit.Omega, "rad"))
         Mₚ:     $(stringify_units(orbit.M_planet, "M⊙"))
         Mₛ:     $(stringify_units(orbit.M_star, "M⊙"))"""
    )
end


text:
using ArrayLayouts: MemoryLayout
using BandedMatrices

import BandedMatrices:
    AbstractBandedMatrix, bandwidths, bandeddata

using LinearAlgebra
import LinearAlgebra:
    LU, ZeroPivotException,
    ldiv!, lu!, lu, factorize

import Base: @propagate_inbounds

"""
    CollocationMatrix{T} <: AbstractBandedMatrix{T}

B-spline collocation matrix, defined by

```math
C_{ij} = b_j(x_i),
```

where ``\\bm{x}`` is a set of collocation points.

Wraps a `BandedMatrix`, providing an efficient LU factorisation without pivoting
adapted from de Boor (1978). The factorisation takes advantage of the total
positivity of spline collocation matrices (de Boor 2001, p. 175).

# Factorisation

`CollocationMatrix` supports in-place LU factorisation using [`lu!`](@ref), as
well as out-of-place factorisation using [`lu`](@ref). LU decomposition may also
be performed via `factorize`.

"""
struct CollocationMatrix{
        T,
        M <: BandedMatrix{T},
    } <: AbstractBandedMatrix{T}
    data :: M
    CollocationMatrix(x::BandedMatrix{T}) where {T} = new{T, typeof(x)}(x)
end

# This affects printing e.g. in the REPL.
MemoryLayout(::Type{<:CollocationMatrix{T,M}}) where {T,M} = MemoryLayout(M)

factorize(A::CollocationMatrix) = lu(A)
bandwidths(A::CollocationMatrix) = bandwidths(parent(A))
bandeddata(A::CollocationMatrix) = bandeddata(parent(A))

Base.parent(A::CollocationMatrix) = A.data
Base.size(A::CollocationMatrix) = size(parent(A))
Base.copy(A::CollocationMatrix) = CollocationMatrix(copy(parent(A)))

# Adapted from BandedMatrices
function Base.array_summary(
        io::IO, A::CollocationMatrix, inds::Tuple{Vararg{Base.OneTo}},
    )
    T = eltype(A)
    print(io, Base.dims2string(length.(inds)), " CollocationMatrix{$T} with bandwidths $(bandwidths(A))")
end

@inline @propagate_inbounds Base.getindex(A::CollocationMatrix, i...) =
    parent(A)[i...]

@inline @propagate_inbounds Base.setindex!(A::CollocationMatrix, v, i...) =
    parent(A)[i...] = v

const CollocationLU{T} = LU{T, <:CollocationMatrix{T}} where {T}

"""
    LinearAlgebra.lu!(C::CollocationMatrix, pivot = Val(false); check = true)

Perform in-place LU factorisation of collocation matrix without pivoting.

Takes advantage of the totally positive property of collocation matrices
appearing in spline calculations (de Boor 1978).

The code is ported from Carl de Boor's BANFAC routine in FORTRAN77, via its
[FORTRAN90 version by John Burkardt](https://people.sc.fsu.edu/~jburkardt/f_src/pppack/pppack.f90).
"""
function lu!(C::CollocationMatrix, _pivot_IGNORED = Val(false); check = true)
    check || throw(ArgumentError("`check = false` not yet supported"))
    if size(C, 1) != size(C, 2)
        throw(DimensionMismatch(
            "factorisation of non-square collocation matrices not supported"
        ))
    end
    w = bandeddata(C)
    nbandl, nbandu = bandwidths(C)
    nrow = size(C, 1)
    nroww = size(w, 1)
    @assert nrow == size(w, 2)
    @assert nroww == nbandl + nbandu + 1
    isempty(C) && error("matrix is empty")
    middle = nbandu + 1  # w[middle, :] contains the main diagonal of A
    Cfact = LU(C, Int[], 0) :: CollocationLU  # factors, ipiv, info

    if nrow == 1
        @inbounds iszero(w[middle, nrow]) && throw(ZeroPivotException(1))
        return Cfact
    end

    if nbandl == 0
        # A is upper triangular. Check that the diagonal is nonzero.
        for i = 1:nrow
            @inbounds iszero(w[middle, i]) && throw(ZeroPivotException(i))
        end
        return Cfact
    end

    if nbandu == 0
        # A is lower triangular. Check that the diagonal is nonzero and
        # divide each column by its diagonal.
        @inbounds for i = 1:(nrow - 1)
            pivot = w[middle, i]
            iszero(pivot) && throw(ZeroPivotException(i))
            ipiv = inv(pivot)
            for j = 1:min(nbandl, nrow - i)
                w[middle + j, i] *= ipiv
            end
        end
        return Cfact
    end

    # A is not just a triangular matrix.
    # Construct the LU factorization.
    @inbounds for i = 1:(nrow - 1)
        pivot = w[middle, i]  # pivot for the i-th step
        iszero(pivot) && throw(ZeroPivotException(i))
        ipiv = inv(pivot)

        # Divide each entry in column `i` below the diagonal by `pivot`.
        for j = 1:min(nbandl, nrow - i)
            w[middle + j, i] *= ipiv
        end

        # Subtract A[i, i+k] * (i-th column) from (i+k)-th column (below row `i`).
        for k = 1:min(nbandu, nrow - i)
            factor = w[middle - k, i + k]
            for j = 1:min(nbandl, nrow - i)
                w[middle - k + j, i + k] -= factor * w[middle + j, i]
            end
        end
    end

    # Check the last diagonal entry.
    @inbounds iszero(w[middle, nrow]) && throw(ZeroPivotException(nrow))

    Cfact
end

"""
    LinearAlgebra.lu(C::CollocationMatrix, pivot = Val(false); check = true)

Returns LU factorisation of collocation matrix.

See also [`lu!`](@ref).
"""
lu(C::CollocationMatrix, args...; kws...) = lu!(copy(C), args...; kws...)

ldiv!(F::CollocationLU, y::AbstractVector) = ldiv!(y, F, y)

# Solution of banded linear system A * x = y.
# The code is adapted from Carl de Boor's BANSLV routine in FORTRAN77, via its
# FORTRAN90 version by John Burkardt.
# Note that `x` and `y` may be the same vector.
function ldiv!(
        x::AbstractVector,
        F::CollocationLU,
        y::AbstractVector,
    )
    A = parent(F.factors) :: BandedMatrix
    w = BandedMatrices.bandeddata(A)
    nrow = size(A, 1)
    nbandl, nbandu = bandwidths(A)
    middle = nbandu + 1
    @assert size(A, 1) == size(A, 2)

    if !(length(x) == length(y) == nrow)
        throw(DimensionMismatch("vectors `x` and `y` must have length $nrow"))
    end

    if x !== y
        copy!(x, y)
    end

    if nrow == 1
        @inbounds x[1] /= w[middle, 1]
        return x
    end

    # Forward pass:
    #
    # For i = 1:(nrow-1), subtract RHS[i] * (i-th column of L) from the
    # right hand side, below the i-th row.
    if nbandl != 0
        for i = 1:(nrow - 1)
            jmax = min(nbandl, nrow - i)
            for j = 1:jmax
                @inbounds x[i + j] -= x[i] * w[middle + j, i]
            end
        end
    end

    # Backward pass:
    #
    # For i = nrow:-1:1, divide RHS[i] by the i-th diagonal entry of
    # U, then subtract RHS[i]*(i-th column of U) from right hand side, above the
    # i-th row.
    @inbounds for i = nrow:-1:2
        x[i] /= w[middle, i]
        for j = 1:min(nbandu, i - 1)
            x[i - j] -= x[i] * w[middle - j, i]
        end
    end

    @inbounds x[1] /= w[middle, 1]

    x
end

function Base.getproperty(F::CollocationLU, d::Symbol)
    factors = getfield(F, :factors) :: CollocationMatrix
    A = parent(factors) :: BandedMatrix
    data = bandeddata(A)
    l, u = bandwidths(A)
    middle = u + 1
    if d === :L
        L = similar(A, size(A)..., l, 0)
        bandeddata(L) .= @view data[middle:(middle + l), :]
        L[Band(0)] .= 1
        LowerTriangular(L)
    elseif d === :U
        U = similar(A, size(A)..., 0, u)
        bandeddata(U) .= @view data[1:middle, :]
        UpperTriangular(U)
    elseif d === :p
        Base.OneTo(size(F, 1))    # no row permutations (no pivoting)
    elseif d === :P
        LinearAlgebra.I  # no row permutations
    else
        getfield(F, d)
    end
end


text:
using DataStructures

import DataStructures.compare

import Base: push!, pop!, isempty, empty!, sizehint!
export push!, pop!, isempty!, empty!, sizehint!

struct Earlier end
compare(c::Earlier, x::Event, y::Event) = time(x) == time(y) ? kind(x) < kind(y) : time(x) < time(y)

#struct EventQueue
#  immediates::CircularDeque{Event}
#  laters::BinaryHeap{Event, Earlier}
#  EventQueue(immediate_size::Integer=128) = new(CircularDeque{Event}(immediate_size), BinaryHeap{Event, Earlier}())
#end

#push!(queue::EventQueue, event::Event; immediate::Bool=false) = begin immediate ? push!(queue.immediates, event) : push!(queue.laters, event); queue end
#pop!(queue::EventQueue) = isempty(queue.immediates) ? pop!(queue.laters) : popfirst!(queue.immediates)
#isempty(queue::EventQueue) = isempty(queue.immediates) & isempty(queue.laters)

struct EventQueue
  queue::BinaryHeap{Event, Earlier}
  EventQueue() = new( BinaryHeap{Event, Earlier}())
end

empty!(queue::EventQueue) = empty!(queue.queue.valtree)
push!(queue::EventQueue, event::Event; immediate::Bool=false) = push!(queue.queue, event)
pop!(queue::EventQueue) = pop!(queue.queue)
isempty(queue::EventQueue) = isempty(queue.queue)
sizehint!(queue::EventQueue, sz::Integer) = sizehint!(queue.queue, sz)