text:
import Base: REPL, LineEdit, Terminals
using RCall
using Compat

type FakeTerminal <: Base.Terminals.UnixTerminal
    in_stream::Base.IO
    out_stream::Base.IO
    err_stream::Base.IO
    hascolor::Bool
    raw::Bool
    FakeTerminal(stdin,stdout,stderr,hascolor=true) =
        new(stdin,stdout,stderr,hascolor,false)
end

Base.Terminals.hascolor(t::FakeTerminal) = t.hascolor
Base.Terminals.raw!(t::FakeTerminal, raw::Bool) = t.raw = raw
Base.Terminals.size(t::FakeTerminal) = (24, 80)

# fake repl

stdin_read,stdin_write = (Base.PipeEndpoint(), Base.PipeEndpoint())
stdout_read,stdout_write = (Base.PipeEndpoint(), Base.PipeEndpoint())
stderr_read,stderr_write = (Base.PipeEndpoint(), Base.PipeEndpoint())
Base.link_pipe(stdin_read,true,stdin_write,true)
Base.link_pipe(stdout_read,true,stdout_write,true)
Base.link_pipe(stderr_read,true,stderr_write,true)

repl = Base.REPL.LineEditREPL(FakeTerminal(stdin_read, stdout_write, stderr_write,false))

repltask = @async begin
    Base.REPL.run_repl(repl)
end

send_repl(x, enter=true) = write(stdin_write, enter? "$x\n" : x)

function check_repl(io::IO, x)
    read_task = @task readuntil(io, x)
    t = Base.Timer((_) -> Base.throwto(read_task,
                ErrorException("Expect \"$x\", but wait too long.")), 5)
    schedule(read_task)
    wait(read_task)
    close(t)
end

check_repl_stdout(x) = check_repl(stdout_read, x)
check_repl_stderr(x) = check_repl(stderr_read, x)

# waiting for the repl
send_repl("using RCall")

RCall.repl_init(repl)

send_repl("\$", false)
check_repl_stdout("R> ")

send_repl("bar = 'apple'")
check_repl_stdout("bar = 'apple'")

send_repl("paste0('pine', bar)")
check_repl_stdout("pineapple")

send_repl("mtca\t", false)
check_repl_stdout("mtcars")

send_repl("")
check_repl_stdout("\n")

send_repl("\\alp\t", false)
check_repl_stdout("\\alpha")

send_repl("\t", false)
check_repl_stdout("α")

send_repl("")
check_repl_stdout("\n")

send_repl("foo]")
check_repl_stderr("unexpected")

send_repl("stop('something is wrong')")
check_repl_stderr("something is wrong")

send_repl("not_found")
check_repl_stderr("not found")

send_repl("\b", false)
check_repl_stdout("julia> ")

send_repl("x = \"orange\"")
check_repl_stdout("x = \"orange\"")

send_repl("\$", false)
check_repl_stdout("R> ")

send_repl("\$x")
check_repl_stdout("orange")

send_repl("\$not_found")
check_repl_stderr("UndefVarError")


text:
module CallGraphs

import Base: hash, isless
export FileFunction, combine, emacs_highlighting, findcallers, findgc, highlightgc, parsedot, parsedots

const rnode = r"Node(\w+)"     # for extracting the Node# from a .dot file
const rlabel = r"\"{(.*)}\""   # for extracting the function name from a "label" in a .dot file
const rfname = r"(\w+)\(.*\)"  # for stripping the signature from a function
const ralphanum = r"\w+"       # for extracting chunks of alphanumeric characters (e.g., function names)

type FileFunction
    file::ASCIIString
    func::ASCIIString
end

isless(cf1::FileFunction, cf2::FileFunction) = cf1.file < cf2.file ? true : (cf1.file == cf2.file && cf1.func < cf2.func)
==(cf1::FileFunction, cf2::FileFunction) = cf1.file == cf2.file && cf1.func == cf2.func
hash(cf::FileFunction, h::UInt) = hash(cf.file, hash(cf.func, h))

### Parsing .dot callgraph files
@doc """
`cgs = parsedots(dirname=".")` parses all callgraph `*.dot` files in the specified directory.
See `combine` for merging them into a single call graph, and `parsedot` for detailed information
about the encoding within `cgs`.
""" ->
function parsedots(dirname=".")
    isdir(dirname) || error("$dirname is not a directory")
    fls = readdir(dirname)
    cgs = Array(Any, 0)
    for fl in fls
        basename, ext = splitext(fl)
        ext == ".dot" || continue
        cg = parsedot(joinpath(dirname, fl))
        push!(cgs, (basename, cg))
    end
    cgs
end

@doc """
`nodes, labels, calls, calledby, externally_defined = parsedot(filename)`

Parses a single callgraph file, stored in GraphViz `*.dot` file, as generated by `llvm`'s
`opt -analyze -dot-callgraph`. Nodes in the graph (functions) are represented as `UInt32`
values, most of which have a `label` (the name, and sometimes the signature, of the function).
`nodes` is a dictionary convering `lbl=>node`, and `labels` performs
the opposite mapping. For a function represented by `node`, `calls[node]` returns the set
of nodes that `node` calls. Conversely, `calledby[node]` returns the list of functions
that call `node`.  Finally, some of the functions lack definitions in the current file, and
are all represented as calling the node `externally_defined`.

This is a low-level interface; usually you want to use `parsedots` and `combine`.
""" ->
function parsedot(io::IO)
    lines = split(readall(io), '\n')
    labels = Dict{UInt32,ASCIIString}()
    nodes = Dict{ASCIIString,UInt32}()
    calls = Dict{UInt32,Set{UInt32}}()
    calledby = Dict{UInt32,Set{UInt32}}()
    local external_caller
    for line in lines
        mnode = match(rnode, line)
        mnode == nothing && continue
        node = UInt32(parseint(mnode.captures[1]))
        if contains(line, "label")
            # A new node
            mlabel = match(rlabel, line)
            label = mlabel.captures[1]
            if label == "external node"
                external_caller = node
                continue
            end
            labels[node] = label
            nodes[label] = node
            calls[node] = Set{UInt32}()
            if !haskey(calledby, node)
                calledby[node] = Set{UInt32}()
            end
        else
            node == external_caller && continue
            mnode = match(rnode, line, mnode.offsets[1])
            mnode == nothing && error("Don't understand format of ", line)
            isempty(mnode.captures) && error("Don't understand format of ", line)
            callee = UInt32(parseint(mnode.captures[1]))
            push!(calls[node], callee)
            if haskey(calledby, callee)
                push!(calledby[callee], node)
            else
                calledby[callee] = Set{UInt32}(node)
            end
        end
    end
    for (key, val) in calledby
        delete!(val, external_caller)
    end
    # Find a node with no label---this corresponds to externally-defined functions
    local externally_defined
    for (nd, s) in calledby
        if !haskey(labels, nd)
            externally_defined = nd
            break
        end
    end
    nodes, labels, calls, calledby, externally_defined
end

function parsedot(filename::AbstractString)
    results = open(filename) do file
        parsedot(file)
    end
    results
end

### Combining callgraphs from many files into a single callgraph
@doc """
`calls, calledby = combine(cgs...)` merges the callgraphs from multiple files into a single
large callgraph. Both `calls` and `calledby` are dictionaries, where the `key` is
a `FileFunction`, a type which encodes the `file` and `func` as strings.

`calls[FileFunction("file1.c", "function8")]` will return the set of functions that `"function8"` calls.
`calledby[FileFunction("file1.c", "function8")]` will return the set of functions calling `"function8"`.

In cases where `opt -analyze -dot-callgraph` retained the function argument signature,
`combine` will discard it.

See also `parsedots`, which generates `cgs`.
""" ->
function combine(cgs...)
    functions = Array(FileFunction, 0)
    lookup = Dict{ASCIIString, FileFunction}()
    for cg in cgs
        filename, results = cg
        nodes, labels, calls, calledby, externally_defined = results
        for (nd,lbl) in labels
            olbl = lbl
            # Discard argument signatures and just keep the function name
            m = match(rfname, lbl)
            if m != nothing
                lbl = m.captures[1]
            end
            func = FileFunction(filename, lbl)
            push!(functions, func)
            lookup[olbl] = func
        end
    end
    calls = Dict{FileFunction, Set{FileFunction}}()
    calledby = Dict{FileFunction, Set{FileFunction}}()
    for cg in cgs
        filename, results = cg
        nodes, labels, calls1, calledby1, externally_defined = results
        for (nd,lbl) in labels
            func = lookup[lbl]
            if !haskey(calls, func)
                calls[func] = Set{FileFunction}()
            end
            cls = calls[func]
            for c in calls1[nd]
                c == externally_defined && continue
                callee = lookup[labels[c]]
                push!(cls, callee)
                if !haskey(calledby, callee)
                    calledby[callee] = Set{FileFunction}()
                end
                push!(calledby[callee], func)
            end
        end
    end
    calls, calledby
end

### Finding the list of functions that call funcname
@doc """
`callers = findcallers(calledby, funcname; indirect=true)` returns a set of `FileFunction`s
that call the function with name `funcname`. `calledby` is the dictionary
returned by `combine`. If `indirect` is true, then `callers` includes functions that might
call `funcname` even indirectly, through other functions.
""" ->
function findcallers(calledby, funcname; indirect::Bool=true)
    local funckey
    found = false
    for (key, callers) in calledby
        if key.func == funcname
            found = true
            funckey = key
        end
    end
    found || error("$funcname not found")
    callers = calledby[funckey]
    push!(callers, funckey)
    # Handle recursive/indirect calls by iterating the call chain until
    # we reach steady-state
    if indirect
        while true
            newcallers = Set{FileFunction}()
            for func in callers
                push!(newcallers, func)
                if haskey(calledby, func)
                    union!(newcallers, calledby[func])
                end
            end
            if isempty(setdiff(newcallers, callers))
                break
            end
            callers = newcallers
        end
    end
    callers
end

@doc """
`names = uniquenames(filefuncs)` takes a vector of `FileFunction`s and returns
a set of unique function names.
""" ->
function uniquenames(filefuncs::Vector{FileFunction})
    names = Set{ASCIIString}()
    for f in filefuncs
        push!(names, f.func)
    end
    names
end

include("gc.jl")

end

ext: jl


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ bcc82d12-6566-11eb-325f-8b6c4ac362a0
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add("PlutoUI")
	using PlutoUI
end

# ╔═╡ 7f2c6b8a-6be9-4c64-b0b5-7fc4435153ee
@bind x Slider(50:100)

# ╔═╡ 2995d591-0f74-44e8-9c06-c42c2f9c68f8
@bind y Slider(x:200)

# ╔═╡ 6bc11e12-3bdb-4ca4-a36d-f8067af95ca5
x

# ╔═╡ 80789650-d01f-4d75-8091-6117a66402cb
y

# ╔═╡ Cell order:
# ╠═bcc82d12-6566-11eb-325f-8b6c4ac362a0
# ╠═7f2c6b8a-6be9-4c64-b0b5-7fc4435153ee
# ╠═2995d591-0f74-44e8-9c06-c42c2f9c68f8
# ╠═6bc11e12-3bdb-4ca4-a36d-f8067af95ca5
# ╠═80789650-d01f-4d75-8091-6117a66402cb


text:
module SurfaceQG

export
  Problem,
  set_b!,
  updatevars!,

  kinetic_energy,
  buoyancy_variance,
  buoyancy_dissipation,
  buoyancy_work,
  buoyancy_advection,
  kinetic_energy_advection

  using
    CUDA,
    Reexport

@reexport using FourierFlows

using LinearAlgebra: mul!, ldiv!
using FourierFlows: parsevalsum

nothingfunction(args...) = nothing

"""
    Problem(; parameters...)

Construct a Surface QG turbulence problem.
"""
function Problem(dev::Device=CPU();
  # Numerical parameters
          nx = 256,
          Lx = 2π,
          ny = nx,
          Ly = Lx,
          dt = 0.01,
  # Hyper-viscosity parameters
           ν = 0,
          nν = 1,
  # Timestepper and equation options
     stepper = "RK4",
       calcF = nothingfunction,
  stochastic = false,
           T = Float64)

  grid = TwoDGrid(dev, nx, Lx, ny, Ly; T=T)

  params = Params{T}(ν, nν, calcF)

  vars = calcF == nothingfunction ? Vars(dev, grid) : (stochastic ? StochasticForcedVars(dev, grid) : ForcedVars(dev, grid))

  equation = Equation(params, grid)

  return FourierFlows.Problem(equation, stepper, dt, grid, vars, params, dev)
end

# ----------
# Parameters
# ----------

"""
    Params(ν, nν, calcF!)

Returns the params for Surface QG turbulence.
"""
struct Params{T} <: AbstractParams
       ν :: T         # Buoyancy viscosity coefficient
      nν :: Int       # Buoyancy hyperviscous order
  calcF! :: Function  # Calculates the Fourier transform of the buoyancy forcing `Fh`
end

Params(ν, nν) = Params(ν, nν, nothingfunction)

# ---------
# Equations
# ---------

"""
    Equation(params, grid)

Returns the equation for Surface QG turbulence with `params` and `grid`.
"""
function Equation(params::Params, grid::AbstractGrid)
  L = @. - params.ν * grid.Krsq^params.nν
  CUDA.@allowscalar L[1, 1] = 0
  return FourierFlows.Equation(L, calcN!, grid)
end

# ----
# Vars
# ----

abstract type SurfaceQGVars <: AbstractVars end

struct Vars{Aphys, Atrans, F, P} <: SurfaceQGVars
        b :: Aphys
        u :: Aphys
        v :: Aphys
       bh :: Atrans
       uh :: Atrans
       vh :: Atrans
       Fh :: F
  prevsol :: P
end

const ForcedVars = Vars{<:AbstractArray, <:AbstractArray, <:AbstractArray, Nothing}
const StochasticForcedVars = Vars{<:AbstractArray, <:AbstractArray, <:AbstractArray, <:AbstractArray}

"""
    Vars(dev, grid)

Returns the vars for unforced surface QG turbulence on device dev and with `grid`.
"""
function Vars(::Dev, grid::AbstractGrid) where Dev
  T = eltype(grid)
  @devzeros Dev T (grid.nx, grid.ny) b u v
  @devzeros Dev Complex{T} (grid.nkr, grid.nl) bh uh vh
  return Vars(b, u, v, bh, uh, vh, nothing, nothing)
end

"""
    ForcedVars(dev, grid)

Returns the vars for forced surface QG turbulence on device `dev` and with
`grid`.
"""
function ForcedVars(dev::Dev, grid) where Dev
  T = eltype(grid)
  @devzeros Dev T (grid.nx, grid.ny) b u v
  @devzeros Dev Complex{T} (grid.nkr, grid.nl) bh uh vh Fh
  return Vars(b, u, v, bh, uh, vh, Fh, nothing)
end

"""
    StochasticForcedVars(dev, grid)

Returns the `vars` for stochastically forced surface QG turbulence on device
`dev` and with `grid`.
"""
function StochasticForcedVars(dev::Dev, grid) where Dev
  T = eltype(grid)
  @devzeros Dev T (grid.nx, grid.ny) b u v
  @devzeros Dev Complex{T} (grid.nkr, grid.nl) bh uh vh Fh prevsol
  return Vars(b, u, v, bh, uh, vh, Fh, prevsol)
end

# -------
# Solvers
# -------

"""
    calcN_advection(N, sol, t, clock, vars, params, grid)

Calculates the advection term.
"""
function calcN_advection!(N, sol, t, clock, vars, params, grid)
  @. vars.bh = sol
  @. vars.uh =   im * grid.l  * sqrt(grid.invKrsq) * sol
  @. vars.vh = - im * grid.kr * sqrt(grid.invKrsq) * sol

  ldiv!(vars.u, grid.rfftplan, vars.uh)
  ldiv!(vars.v, grid.rfftplan, vars.vh)
  ldiv!(vars.b, grid.rfftplan, vars.bh)

  ub, ubh = vars.u, vars.uh         # use vars.u, vars.uh as scratch variables
  vb, vbh = vars.v, vars.vh         # use vars.v, vars.vh as scratch variables

  @. ub *= vars.b # u*b
  @. vb *= vars.b # v*b

  mul!(ubh, grid.rfftplan, ub) # \hat{u*b}
  mul!(vbh, grid.rfftplan, vb) # \hat{v*b}

  @. N = - im * grid.kr * ubh - im * grid.l * vbh
  return nothing
end

function calcN!(N, sol, t, clock, vars, params, grid)
  calcN_advection!(N, sol, t, clock, vars, params, grid)
  addforcing!(N, sol, t, clock, vars, params, grid)
  return nothing
end

addforcing!(N, sol, t, clock, vars::Vars, params, grid) = nothing

function addforcing!(N, sol, t, clock, vars::ForcedVars, params, grid)
  params.calcF!(vars.Fh, sol, t, clock, vars, params, grid)
  @. N += vars.Fh
  return nothing
end

function addforcing!(N, sol, t, clock, vars::StochasticForcedVars, params, grid)
  if t == clock.t # not a substep
    @. vars.prevsol = sol # sol at previous time-step is needed to compute budgets for stochastic forcing
    params.calcF!(vars.Fh, sol, t, clock, vars, params, grid)
  end
  @. N += vars.Fh
  return nothing
end

# ----------------
# Helper functions
# ----------------

"""
    updatevars!(prob)

Update variables in `vars` with solution in `sol`.
"""
function updatevars!(prob)
  vars, grid, sol = prob.vars, prob.grid, prob.sol

  @. vars.bh = sol
  @. vars.uh =   im * grid.l  * sqrt(grid.invKrsq) * sol
  @. vars.vh = - im * grid.kr * sqrt(grid.invKrsq) * sol

  ldiv!(vars.b, grid.rfftplan, deepcopy(vars.bh))
  ldiv!(vars.u, grid.rfftplan, deepcopy(vars.uh))
  ldiv!(vars.v, grid.rfftplan, deepcopy(vars.vh))

  return nothing
end

"""
    set_b!(prob, b)

Set the solution `sol` as the transform of b and update all variables.
"""
function set_b!(prob, b)
  mul!(prob.sol, prob.grid.rfftplan, b)
  CUDA.@allowscalar prob.sol[1, 1] = 0 # zero domain average

  updatevars!(prob)

  return nothing
end

"""
    kinetic_energy(prob)

Returns the domain-averaged surface kinetic energy. In SQG, this is
identical to half the domain-averaged surface buoyancy variance.
"""
@inline function kinetic_energy(prob)
  sol, vars, grid = prob.sol, prob.vars, prob.grid

  @. vars.uh =   im * grid.l  * sqrt(grid.invKrsq) * sol
  @. vars.vh = - im * grid.kr * sqrt(grid.invKrsq) * sol

  kinetic_energyh = vars.bh         # use vars.bh as scratch variable

  @. kinetic_energyh = 0.5 * (abs2(vars.uh) + abs2(vars.vh)) # ½(|û|²+|v̂|²)

  return 1 / (grid.Lx * grid.Ly) * parsevalsum(kinetic_energyh, grid)
end

"""
    buoyancy_variance(prob)

Returns the domain-averaged buoyancy variance. In SQG flows this is identical to
the domain-averaged velocity variance (twice the kinetic energy)
"""
@inline function buoyancy_variance(prob)
  sol, grid = prob.sol, prob.grid

  return 1 / (grid.Lx * grid.Ly) * parsevalsum(abs2.(sol), grid)
end

"""
    buoyancy_dissipation(prob)

Returns the domain-averaged dissipation rate of surface buoyancy variance due
to small scale diffusion/viscosity. nν must be >= 1.

In SQG, this is identical to twice the rate of kinetic energy dissipation
"""
@inline function buoyancy_dissipation(prob)
  sol, vars, params, grid = prob.sol, prob.vars, prob.params, prob.grid
  buoyancy_dissipationh = vars.uh         # use vars.uh as scratch variable

  @. buoyancy_dissipationh = 2 * params.ν * grid.Krsq^params.nν * abs2(sol)
  CUDA.@allowscalar buoyancy_dissipationh[1, 1] = 0

  return 1 / (grid.Lx * grid.Ly) * parsevalsum(buoyancy_dissipationh, grid)
end

"""
    buoyancy_work(prob)
    buoyancy_work(sol, vars, grid)

Returns the domain-averaged rate of work of buoyancy variance by the forcing Fh.
"""
@inline function buoyancy_work(sol, vars::ForcedVars, grid)
  buoyancy_workh = vars.uh         # use vars.uh as scratch variable

  @. buoyancy_workh =  2 * sol * conj(vars.Fh)    # 2*b̂*conj(f̂)
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(buoyancy_workh, grid)
end

@inline function buoyancy_work(sol, vars::StochasticForcedVars, grid)
  buoyancy_workh = vars.uh         # use vars.uh as scratch variable

  @. buoyancy_workh =  (vars.prevsol + sol) * conj(vars.Fh) # Stratonovich
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(buoyancy_workh, grid)
end

@inline buoyancy_work(prob) = buoyancy_work(prob.sol, prob.vars, prob.grid)

end # module


text:
# Tests for empty Interval
I = Interval()
# J = [1,2]
J = Interval(0.0,1.0)
@test string(I) == "∅"
@test empty(I)
@test 0 ∉ I
@test 1.3 ∉ I
@test left(I) == 0
@test right(I) == 0
@test I ⊆ J
@test J ⊈ I
@test I ⊆ I
@test cardinal(I) == 0
@test !disjoint(I,J)
@test !disjoint(J,I)
@test I ∩ J == Interval()
@test J ∩ I == Interval()
@test I ∪ J == J
@test J ∪ I == J
@test compact(I) == Interval()
@test compactness(I) == 0
@test superset(I) == Interval()
@test jaccard(I,J) == 0
@test jaccard(J,I) == 0
@test jaccard(I,I) == 1
@test overlap_coefficient(I,J) == 0
@test overlap_coefficient(J,I) == 0
@test overlap_coefficient(I,I) == 1
@test dice_coefficient(I,J) == 0
@test dice_coefficient(J,I) == 0
@test dice_coefficient(I,I) == 1

# I = [0,1]
I = Interval(0.0,1.0)
@test string(I) == "[0.0,1.0]"
@test !empty(I)

# J = ]0,1]
J = Interval(0.0,true,1.0)
@test string(J) == "]0.0,1.0]"
@test !empty(J)

# K = ]1,2[
K = Interval(1,true,2.0,true)
@test string(K) == "]1,2.0["
@test !empty(K)

# L = [1,2[
L = Interval(1,2,true)
@test string(L) == "[1,2["
@test !empty(L)

# M = [1,7[
M = Interval(1,7,true)
@test string(M) == "[1,7["
@test !empty(M)

# N = ]-1.3,2.6[
N = Interval(-1.3,true,2.6,true)
@test string(N) == "]-1.3,2.6["
@test !empty(N)

# Tests for left and right limits
@test left(I) == 0.0
@test right(I) == 1
@test left(J) == 0.0
@test right(L) == 2.0
@test left(L) == right(I)

# Tests for open/close
@test !I.open_left
@test !I.open_right
@test J.open_left
@test !J.open_right
@test K.open_left
@test K.open_right

# Equality tests
@test I == I
@test J == Interval(0.0,true,1.0)
@test I != J
@test J != K

# Tests for ∈
@test -0.1 ∉ I
@test 0.0 ∈ I
@test 0.5 ∈ I
@test 1.0 ∈ I
@test 1.00001 ∉ I
@test 0.0 ∉ J
@test 0 ∉ J
@test 0.000000001 ∈ J
@test 0.9 ∈ J
@test 2.0 ∉ K
@test 2 ∉ K
@test 1.999999999 ∈ K

# Tests for cardinal
@test cardinal(I) == 1.0
@test cardinal(J) == 1.0
@test cardinal(K) == 1.0
@test cardinal(L) == 1.0
@test cardinal(M) == 6.0
@test isapprox(cardinal(N),3.9)
@test cardinal(Interval(0,Inf,true)) == Inf
@test cardinal(Interval(-Inf,true,2)) == Inf
@test cardinal(Interval(-Inf,true,Inf,true)) == Inf

# Tests for ⊆
@test J ⊆ I
@test I ⊆ I
@test I ⊈ J
@test I ⊆ N
@test M ⊈ N
@test Interval(1,true,2.5) ⊆ Interval(1,true,2.6)
# [2.2,3[ ⊆ [2,3[
@test Interval(2.2,3,true) ⊆ Interval(2,3,true)
# [2.2,3] ⊈ [2,3[
@test Interval(2.2,3) ⊈ Interval(2,3,true)

# Tests for disjoint
# [0,1] and [0.5,0.9] are not disjoint
@test !disjoint(Interval(0.0,1.0),Interval(0.5,0.9))
# [0.5,0.9] and [0,1] are not disjoint
@test !disjoint(Interval(0.5,0.9),Interval(0.0,1.0))
# [0.0,1.0] and [1.1,2.0] are disjoint
@test disjoint(Interval(0.0,1.0),Interval(1.1,2.0))
# [1.1,2.0] and [0.0,1.0] are disjoint
@test disjoint(Interval(1.1,2.0),Interval(0.0,1.0))
# [0,1] and [1,2] are not disjoint
@test !disjoint(Interval(0.0,1.0),Interval(1.0,2.0))
# [0,1[ and [1,2] are disjoint
@test disjoint(Interval(0.0,1.0,true),Interval(1.0,2.0))
# [0,1] and ]1,2] are disjoint
@test disjoint( Interval(0.0,1.0), Interval(1.0,true,2.0))
# ]1,2[ and [0,1[ are disjoint
@test disjoint(Interval(1.0,true,2.0,true), Interval(0.0,1.0,true))
# ]0,2[ and ]0,1[ are not disjoint
@test !disjoint(Interval(0.0,true,2.0,true), Interval(0.0,true,1.0,true))
# ]-1,1.5[ and [1.5,1.67[ are disjoint
@test disjoint(Interval(-1.0,true,1.5,true),Interval(1.5,false,1.67,true))

# Tests for ∩
# [0,1] ∩ [0.5,1.5] = [0.5,1]
@test Interval(0.0,1.0) ∩ Interval(0.5,1.5) == Interval(0.5,1.0)
# [0.5,1.5] ∩ [0,1] = [0.5,1]
@test Interval(0.5,1.5) ∩ Interval(0.0,1.0) == Interval(0.5,1.0)
# [0,1] ∩ [1,1.5] = [1,1]
@test Interval(0.0,1.0) ∩ Interval(1.0,1.5) == Interval(1.0,1.0)
# [0,3.5] ∩ [1.2,1.6] = [1.2,1.6]
@test Interval(0.0,3.5) ∩ Interval(1.2,1.6) == Interval(1.2,1.6)
# [1.2,1.6] ∩ [0.0,3.5] = [1.2,1.6]
@test Interval(1.2,1.6) ∩ Interval(0.0,3.5) == Interval(1.2,1.6)
# [0,1.5[ ∩ ]1.5,1.6] = ]0,0[
@test Interval(0.0,1.5,true) ∩ Interval(1.5,true,1.6) == Interval(0.0,true,0.0,true)
# [0,1] ∩ ]1,2] = ]0,0[
@test Interval(0,1) ∩ Interval(1,true,2) == Interval(0,true,0,true)

# Tests for ∪
# [0,1] ∪ [0.5,1.5] = [0,1.5]
@test Interval(0.0,1.0) ∪ Interval(0.5,1.5) == Interval(0.0,1.5)
# [0.5,1.5] ∪ [0,1] = [0,1.5]
@test Interval(0.5,1.5) ∪ Interval(0.0,1.0) == Interval(0.0,1.5)
# [0,1] ∪ [1.0,1.5] = [0,1.5]
@test Interval(0.0,1.0) ∪ Interval(1.0,1.5) == Interval(0.0,1.5)
# [0,3.5] ∪ [0.5,1] = [0,3.5]
@test Interval(0.0,3.5) ∪ Interval(0.5,1.0) == Interval(0.0,3.5)
# [0.5,1] ∪ [0,3.5] = [0,3.5]
@test Interval(0.5,1.0) ∪ Interval(0.0,3.5) == Interval(0.0,3.5)
# [0,1] ∪ ]1,2] = [0,2]
@test Interval(0,1) ∪ Interval(1,true,2) == Interval(0,2)
# [1,3] ∪ ]3,4[ = [1,4[
@test Interval(3,true,4,true) ∪ Interval(1,3) == Interval(1,4,true)

# Tests for sampling intervals
@test 0 <= sample(Interval(0,1)) <= 1
@test 0 < sample(Interval(0,true,1,true)) < 1
samples = sample(Interval(2,true,6),100)
@test length(samples) == 100
for s in samples
	@test 2 < s <= 6
end

# Tests for superset
@test superset(I) == I
@test superset(J) == J
@test superset(K) == K

# Tests for compact
@test compact(I) == Interval(0,1)
@test compact(J) == Interval(0,1)
@test compact(K) == Interval(1,2)
@test compact(L) == Interval(1,2)
@test compact(M) == Interval(1,7)
@test compact(N) == Interval(-1.3,2.6)

# Tests for compactness
@test compactness(I) == 1
@test compactness(J) == 1
@test compactness(K) == 1

# Tests for jaccard similarity
# jaccard([0,1],[0,1]) = 1
@test jaccard(Interval(0,1),Interval(0,1)) == 1.0
# jaccard(]0,1[,[0,1]) = 1
@test jaccard(Interval(0,true,1,true),Interval(0,1)) == 1.0
# jaccard([0,2],[0,1]) = 1/2
@test jaccard(Interval(0,2),Interval(0,1)) == 0.5
# jaccard([-1,2],[0,1]) = 1/3
@test isapprox(jaccard(Interval(-1,2),Interval(0,1)),0.333333333333)
# jaccard([0,1],[1,2]) = 0
@test jaccard(Interval(0,1),Interval(1,2)) == 0.0
# jaccard([0,1],]1,2]) = 0
@test jaccard(Interval(0,1),Interval(1,true,2)) == 0.0
# jaccard([0,1],[1,1]) = 0
@test jaccard(Interval(0,1),Interval(1,1)) == 0.0
# jaccard([0,1],]1,1[) = 0
@test jaccard(Interval(0,1),Interval(1,true,1,true)) == 0.0
# jaccard(]-Inf,0],[0,Inf[) = 0
@test jaccard(Interval(-Inf,true,0),Interval(0,Inf,true)) == 0.0
# jaccard(]-Inf,Inf[,]-Inf,Inf[) = 1
@test jaccard(Interval(-Inf,true,Inf,true),Interval(-Inf,true,Inf,true)) == 1.0

# Tests for overlap coefficient
# overlap_coefficient([0,1],[0,1]) = 1
@test overlap_coefficient(Interval(0,1),Interval(0,1)) == 1.0
# overlap_coefficient(]0,1[,[0,1]) = 1
@test overlap_coefficient(Interval(0,true,1,true),Interval(0,1)) == 1.0
# overlap_coefficient([0,2],[0,1]) = 1
@test overlap_coefficient(Interval(0,2),Interval(0,1)) == 1
# overlap_coefficient([-1,2],[0,1]) = 1
@test overlap_coefficient(Interval(-1,2),Interval(0,1)) == 1
# overlap_coefficient([0,1],[1,2]) = 0
@test overlap_coefficient(Interval(0,1),Interval(1,2)) == 0.0
# overlap_coefficient([0,1],]1,2]) = 0
@test overlap_coefficient(Interval(0,1),Interval(1,true,2)) == 0.0
# overlap_coefficient([0,1],[1,1]) = 0
@test overlap_coefficient(Interval(0,1),Interval(1,1)) == 0.0
# overlap_coefficient([0,1],]1,1[) = 0
@test overlap_coefficient(Interval(0,1),Interval(1,true,1,true)) == 0.0
# overlap_coefficient(]-Inf,0],[0,Inf[) = 0
@test overlap_coefficient(Interval(-Inf,true,0),Interval(0,Inf,true)) == 0.0
# overlap_coefficient(]-Inf,Inf[,]-Inf,Inf[) = 1
@test overlap_coefficient(Interval(-Inf,true,Inf,true),Interval(-Inf,true,Inf,true)) == 1.0

# Tests for dice coefficient
# dice_coefficient([0,1],[0,1]) = 1
@test dice_coefficient(Interval(0,1),Interval(0,1)) == 1.0
# dice_coefficient(]0,1[,[0,1]) = 1
@test dice_coefficient(Interval(0,true,1,true),Interval(0,1)) == 1.0
# dice_coefficient([0,2],[0,1]) = 2/3
@test dice_coefficient(Interval(0,2),Interval(0,1)) == 2/3
# dice_coefficient([-1,2],[0,1]) = 1/2
@test dice_coefficient(Interval(-1,2),Interval(0,1)) == 1/2
# dice_coefficient([0,1],[1,2]) = 0
@test dice_coefficient(Interval(0,1),Interval(1,2)) == 0.0
# dice_coefficient([0,1],]1,2]) = 0
@test dice_coefficient(Interval(0,1),Interval(1,true,2)) == 0.0
# dice_coefficient([0,1],[1,1]) = 0
@test dice_coefficient(Interval(0,1),Interval(1,1)) == 0.0
# dice_coefficient([0,1],]1,1[) = 0
@test dice_coefficient(Interval(0,1),Interval(1,true,1,true)) == 0.0
# dice_coefficient(]-Inf,0],[0,Inf[) = 0
@test dice_coefficient(Interval(-Inf,true,0),Interval(0,Inf,true)) == 0.0
# dice_coefficient(]-Inf,Inf[,]-Inf,Inf[) = 1
@test dice_coefficient(Interval(-Inf,true,Inf,true),Interval(-Inf,true,Inf,true)) == 1.0


text:
module CurrencyAmounts

export Currency, CurrencyAmount, ExchangeRate, @currencies

import Base: +, -, *, /, ≈, zero, isless, isnan, isfinite
import Statistics

using RecipesBase

"""
Currency type, e.g. EUR, USD.
"""
struct Currency{T} end

Currency(c:: Symbol) = Currency{c}()
Currency(c:: AbstractString) = Currency(Symbol(c))

Base.Broadcast.broadcastable(c:: Currency) = Ref(c) # treat it as a scalar in broadcasting
Base.show(io:: IO, c:: Currency{T}) where {T} = print(io, string(T))

"""
    @currencies(syms)

Creates one or more currency definitions, these are global constants of type `Currency{T}`.

Example:

    @currencies EUR, USD, GBP
    @currencies EUR

the latter is equivalent to

    const EUR = Currency(:EUR)

"""
macro currencies(syms)
    args = syms isa Expr ? syms.args : [syms]
    for nam in args
        ccy = CurrencyAmounts.Currency{nam}()
        @eval __module__ const $nam = $ccy
    end
end

"""
    CurrencyAmount(x:: Number, c:: Currency)

Amount of a specific currency, i.e. a combination of a numeric value and a currency type.

Examples:

    CurrencyAmount(100.5, EUR) # assumes `@currencies EUR` has been called before
    100.5EUR # equivalent constructor
"""
struct CurrencyAmount{T <: Number, C <: Currency}
    amount:: T
end

CurrencyAmount(x:: Number, ::C) where C <: Currency = CurrencyAmount{typeof(x), C}(x)

Base.Broadcast.broadcastable(c:: CurrencyAmount) = Ref(c) # treat it as a scalar in broadcasting
Base.show(io:: IO, c:: CurrencyAmount{<: Number, C}) where {C} = print(io, c.amount, " ", C())
≈(x:: CurrencyAmount{T1, C}, y:: CurrencyAmount{T2, C}) where {C <: Currency, T1 <: Number, T2 <: Number} = x.amount ≈ y.amount
zero(c ::CurrencyAmount) = zero(typeof(c))
zero(::Type{CurrencyAmount{T, C}}) where {C <: Currency, T <: Number} = CurrencyAmount(zero(T), C())
isless(c1:: CurrencyAmount{<: Number, C}, c2:: CurrencyAmount{<: Number, C}) where {C} = c1.amount < c2.amount
Statistics.middle(c1 ::CurrencyAmount, c2:: CurrencyAmount) = (c1 + c2)/2
isnan(c ::CurrencyAmount) = isnan(c.amount)
isfinite(c ::CurrencyAmount) = isfinite(c.amount)

# construction of CurrencyAmount using multiplication
*(x:: Number, c:: Currency) = CurrencyAmount(x, c)

## CurrencyAmount arithmetics

# Addition and subtraction are only defined between equal currencies
+(x:: CurrencyAmount{T1, C}, y:: CurrencyAmount{T2, C}) where {C <: Currency, T1 <: Number, T2 <: Number} = CurrencyAmount(x.amount + y.amount, C())
-(x:: CurrencyAmount{T1, C}, y:: CurrencyAmount{T2, C}) where {C <: Currency, T1 <: Number, T2 <: Number} = CurrencyAmount(x.amount - y.amount, C())
-(x:: CurrencyAmount{T1, C}) where {C <: Currency, T1 <: Number} = CurrencyAmount(-x.amount, C())

*(x:: CurrencyAmount{T1, C}, y:: Number) where {C <: Currency, T1 <: Number} = CurrencyAmount(x.amount * y, C())
*(x:: Number, y:: CurrencyAmount) = y*x

# only division of a CurrencyAmount by a scalar is defined, not vice versa because this would give dimension of 1/Currency.
/(x:: CurrencyAmount{T1, C}, y:: Number) where {C <: Currency, T1 <: Number} = CurrencyAmount(x.amount / y, C())

# division of 2 amounts with the same currency gives a dimensionless result
/(x:: CurrencyAmount{<: Number, C}, y:: CurrencyAmount{<: Number, C}) where {C <: Currency} = x.amount / y.amount

"""
    ExchangeRate(x:: Number, ::BaseCurrency, ::QuoteCurrency)

Exchange rate between two currencies.
"""
struct ExchangeRate{T <: Number, BaseCurrency <: Currency, QuoteCurrency <: Currency}
    rate:: T
end

Base.Broadcast.broadcastable(c:: ExchangeRate) = Ref(c) # treat it as a scalar in broadcasting
Base.show(io:: IO, c:: ExchangeRate{<: Number, C1, C2}) where {C1, C2} = print(io, c.rate, " ", C2(), "/", C1())

## Exchange Rates
ExchangeRate(x:: Number, ::BaseCurrency, ::QuoteCurrency) where {BaseCurrency <: Currency, QuoteCurrency <: Currency} = ExchangeRate{typeof(x), BaseCurrency, QuoteCurrency}(x)
ExchangeRate(x:: Number, ::C, ::C) where {C <: Currency} = error("currencies must be different")
isless(x1:: ExchangeRate{<: Number, C1, C2}, x2:: ExchangeRate{<: Number, C1, C2}) where {C1, C2} = x1.rate < x2.rate
Statistics.middle(c1 ::ExchangeRate, c2:: ExchangeRate) = (c1 + c2)/2
isnan(x ::ExchangeRate) = isnan(x.rate)
isfinite(x ::ExchangeRate) = isfinite(x.rate)

≈(x:: ExchangeRate{<: Number, C1, C2}, y:: ExchangeRate{<: Number, C1, C2}) where {C1 <: Currency, C2 <: Currency} = x.rate ≈ y.rate

# in a currency fraction, the QuoteCurrency is the numerator and the BaseCurrency the denominator
/(x:: CurrencyAmount{T1, C1}, y:: CurrencyAmount{T2, C2}) where {C1 <: Currency, C2 <: Currency, T1 <: Number, T2 <: Number} = ExchangeRate(x.amount/y.amount, C2(), C1())
/(x:: CurrencyAmount, y:: Currency) = x / CurrencyAmount(1, y) # allowing syntax 1.2EUR/USD

# exchange rate arithmetics
*(x:: ExchangeRate{T1, C1, C2}, y:: Number) where {C1 <: Currency, C2 <: Currency, T1 <: Number} = ExchangeRate(x.rate * y, C1(), C2())
*(x:: Number, y:: ExchangeRate) = y*x

/(x:: ExchangeRate{T1, C1, C2}, y:: Number) where {C1 <: Currency, C2 <: Currency, T1 <: Number} = ExchangeRate(x.rate / y, C1(), C2())
/(y:: Number, x:: ExchangeRate{T1, C1, C2}) where {C1 <: Currency, C2 <: Currency, T1 <: Number} = ExchangeRate(y / x.rate, C2(), C1())

+(x:: ExchangeRate{<: Number, C1, C2}, y:: ExchangeRate{<: Number, C1, C2}) where {C1 <: Currency, C2 <: Currency} = ExchangeRate(x.rate + y.rate, C1(), C2())
-(x:: ExchangeRate{<: Number, C1, C2}, y:: ExchangeRate{<: Number, C1, C2}) where {C1 <: Currency, C2 <: Currency} = ExchangeRate(x.rate - y.rate, C1(), C2())

# division of 2 consistent exchange rates gives a scalar
/(x:: ExchangeRate{<: Number, C1, C2}, y:: ExchangeRate{<: Number, C1, C2}) where {C1 <: Currency, C2 <: Currency} = x.rate / y.rate

# exchange rate triangulation
/(x:: ExchangeRate{<: Number, BaseCurrency1, QuoteCurrency}, y:: ExchangeRate{<: Number, BaseCurrency2, QuoteCurrency}) where {QuoteCurrency <: Currency, BaseCurrency1 <: Currency, BaseCurrency2 <: Currency} = ExchangeRate(x.rate / y.rate, BaseCurrency1(), BaseCurrency2())
/(x:: ExchangeRate{<: Number, BaseCurrency, QuoteCurrency1}, y:: ExchangeRate{<: Number, BaseCurrency, QuoteCurrency2}) where {QuoteCurrency1 <: Currency, QuoteCurrency2 <: Currency, BaseCurrency <: Currency} = ExchangeRate(x.rate / y.rate, QuoteCurrency2(), QuoteCurrency1())
*(x:: ExchangeRate{<: Number, C_both, C1}, y:: ExchangeRate{<: Number, C2, C_both}) where {C_both <: Currency, C1 <: Currency, C2 <: Currency} = ExchangeRate(x.rate * y.rate, C2(), C1())
*(x:: ExchangeRate{<: Number, C1, C_both}, y:: ExchangeRate{<: Number, C_both, C2}) where {C_both <: Currency, C1 <: Currency, C2 <: Currency} = ExchangeRate(x.rate * y.rate, C1(), C2())

## currency conversions

# conversions using arithmetics
*(x:: CurrencyAmount{<: Number, C1}, y:: ExchangeRate{<: Number, C1, C2}) where {C1 <: Currency, C2 <: Currency} = CurrencyAmount(x.amount * y.rate, C2())
*(x:: ExchangeRate, y:: CurrencyAmount) = y*x

/(x:: CurrencyAmount{<: Number, C2}, y:: ExchangeRate{<: Number, C1, C2}) where {C1 <: Currency, C2 <: Currency} = CurrencyAmount(x.amount / y.rate, C1())

# explicit "smart" conversions

"""
   convert(target_currency:: Currency, amount:: CurrencyAmount, exchange_rate:: ExchangeRate):: CurrencyAmount

Converts the amount` to target currency using the given `exchange_rate`. If a conversion is not possible using this exchange rate, `missing` is returned.
"""
Base.convert(:: C2, amount:: CurrencyAmount{<: Number, C1}, rate:: ExchangeRate{<: Number, C1, C2}) where {C1 <: Currency, C2 <: Currency} = amount * rate
Base.convert(:: C1, amount:: CurrencyAmount{<: Number, C2}, rate:: ExchangeRate{<: Number, C1, C2}) where {C1 <: Currency, C2 <: Currency} = amount / rate
Base.convert(:: C1, amount:: CurrencyAmount{<: Number, C1}, :: ExchangeRate) where {C1 <: Currency} = amount # source and target currency are equal, do nothing
Base.convert(::Currency, ::CurrencyAmount, ::ExchangeRate) = missing

"""
   convert(target_currency:: Currency, amount:: CurrencyAmount, exchange_rate_iterable):: CurrencyAmount

Converts the amount` to target currency using one of the given rates in the `exchange_rate_iterable`. If a conversion is not possible using one of the exchange rate, `missing` is returned.
"""
function Base.convert(:: C1, amount:: CurrencyAmount{<: Number, C2}, itr) where {C1 <: Currency, C2 <: Currency}
    for exchange_rate in itr
        exchange_rate isa ExchangeRate || error("iterable must contain exchange rates")
        res = convert(C1(), amount, exchange_rate)
        if !ismissing(res)
            return res
        end
    end
    return missing
end

# plotting
@recipe function plot(x, y:: AbstractVector{CurrencyAmount{T, C}}) where {T, C <: Currency}
    yguide --> C()
    x, getproperty.(y, :amount)
end

@recipe function plot(x, y:: AbstractVector{ExchangeRate{T, BaseCurrency, QuoteCurrency}}) where {T, BaseCurrency <: Currency, QuoteCurrency <: Currency}
    yguide --> string(QuoteCurrency(), " / ", BaseCurrency())
    x, getproperty.(y, :rate)
end

end


text:
# getting elements
function get_id(el::MyXMLElement)
    return el.id
end

function get_data(bx::BEASTXMLElement)
    return find_element(bx, DataXMLElement)
end

function get_newick(bx::BEASTXMLElement)
    return find_element(bx, NewickXMLElement)
end

function get_treeModel(bx::BEASTXMLElement)
    return find_element(bx, TreeModelXMLElement)
end

function get_mbd(bx::BEASTXMLElement)
    return find_element(bx, MBDXMLElement)
end

function get_extension(bx::BEASTXMLElement)
    return find_element(bx, ModelExtensionXMLElement)
end

function get_repeatedMeasures(bx::BEASTXMLElement)
    return find_element(bx, RepeatedMeasuresXMLElement)
end

function get_integratedFactorModel(bx::BEASTXMLElement)
    return find_element(bx, IntegratedFactorsXMLElement)
end

function get_latentFactorModel(bx::BEASTXMLElement)
    return find_element(bx, LatentFactorModelXMLElement)
end

function get_traitLikelihood(bx::BEASTXMLElement)
    return find_element(bx, TraitLikelihoodXMLElement)
end

function get_operators(bx::BEASTXMLElement)
    ops = find_element(bx, OperatorsXMLElement)
    return ops.els
end

function get_mcmc(bx::BEASTXMLElement)
    return find_element(bx, MCMCXMLElement)
end

function get_timer(bx::BEASTXMLElement)
    return find_element(bx, TimerXMLElement)
end

function get_loadings_op(bx::BEASTXMLElement; component::String = "")
    loadings = get_loadings(bx)
    if typeof(loadings) <: ScaledOrthogonalMatrix
        if component == "scale"
            loadings = loadings.scale
        elseif component == "matrix"
            loadings = loadings.U
        else
            throw(ArgumentError("The loadings is a " * name(loadings) *
                " with components \"scale\" and \"matrix\"." *
                " You have set component=\"$component\"."))
        end
    end

    ops = get_operators(bx)
    for op in ops
        if get_parameter(op) === loadings
            return op
        end
    end
    error("No loadings operator.")
end

function get_op_by_parameter(bx::BEASTXMLElement, parameter::MyXMLElement)
    ops = get_operators(bx)
    for op in ops
        if get_parameter(op) === parameter
            return op
        end
    end
    error("cannot find operator for parameter with id '$(get_id(parameter))'")
end

function get_multiplicative_gamma_op(bx::BEASTXMLElement)
    ops = get_operators(bx)
    for op in ops
        t = typeof(op)
        if t <: NormalGammaPrecisionOperatorXMLElement
            if typeof(op.ggp) <: MultipilcativeGammaGibbsProvider
                return op
            end
        end
    end
    error("No multiplicative gamma gibbs operator.")
end

function get_loadings_scale(bx::BEASTXMLElement)
    ifm = get_integratedFactorModel(bx)
    return get_loadings_scale(ifm)
end

function get_factor_model(bx::BEASTXMLElement)
    fac_models = [IntegratedFactorsXMLElement, LatentFactorModelXMLElement]
    for model in fac_models
        fac_model = find_element(bx, model)
        if !isnothing(fac_model)
            return fac_model
        end
    end

    return nothing
end

function get_loadings(bx::BEASTXMLElement)
    fac_model = get_factor_model(bx)
    loadings = get_loadings(fac_model)
    return loadings
end

function get_factor_precisions(bx::BEASTXMLElement)
    fac_model = get_factor_model(bx)
    precs = get_precision(fac_model)
    return precs
end


text:
## subs
"""
Substitute values into a symbolic expression.

Examples
```
@syms x y
ex = x^2 + y^2
subs(ex, x, 1) # 1 + y^2
subs(ex, (x, 1)) # 1 + y^2
subs(ex, (x, 1), (y,x)) # 1 + x^2, values are substituted left to right.
subs(ex, x=>1)  # alternate to subs(x, (x,1))
subs(ex, x=>1, y=>1) # ditto
```
"""
function subs{T<:SymbolicType, S<:SymbolicType}(ex::T, var::S, val)
    s = Basic()
    ccall((:basic_subs2, libsymengine), Void, (Ptr{Basic}, Ptr{Basic}, Ptr{Basic}, Ptr{Basic}), &s, &ex, &var, &val)
    return s
end
subs{T <: SymbolicType, S<:SymbolicType}(ex::T, y::Tuple{S, Any}) = subs(ex, y[1], y[2])
subs{T <: SymbolicType, S<:SymbolicType}(ex::T, y::Tuple{S, Any}, args...) = subs(subs(ex, y), args...)
subs{T <: SymbolicType}(ex::T, d::Pair...) = subs(ex, [(p.first, p.second) for p in d]...)

## Lambdify

## Mapping of SymEngine Constants into julia values
constant_map = Dict("pi" => :pi, "E" => :e, "EulerGamma" => :γ)

## Map symengine classes to function names
fn_map = Dict(
              :Add => :+,
              :Sub => :-,
              :Mul => :*,
              :Div => :/,
              :Pow => :^,
              :re  => :real,
              :im  => :imag,
              :Abs => :abs # not really needed as default in map_fn covers this
              )

map_fn(key, fn_map) = haskey(fn_map, key) ? fn_map[key] : Symbol(lowercase(string(key)))

function walk_expression(ex)
    fn = get_symengine_class(ex)

    if fn == :Symbol
        return Symbol(toString(ex))
    elseif fn in [:Integer , :Rational]
        return N(ex)
    elseif fn == :Complex
        ## hacky
        x = eval(parse(replace(toString(ex), "I", "im")))
        return Expr(:call, :complex, real(x), imag(x))
    elseif fn == :Constant
        return constant_map[toString(ex)]

    end

    as = get_args(ex)

    Expr(:call, map_fn(fn, fn_map), [walk_expression(a) for a in as]...)
end

## evaluate symbolless expression or return afunction
function lambdify(ex::Basic)
    vars = free_symbols(ex)
    if length(vars) == 0
        _lambdify(ex)
    else
        _lambdify(ex, vars)
    end
end
lambdify(ex::BasicType) = lambdify(Basic(ex))

## return a number
function _lambdify(ex)
    body = walk_expression(ex)
    eval(body)
end

## return a function
function _lambdify(ex::Basic, vars)
    body = walk_expression(ex)

    try
        eval(Expr(:function,
                  Expr(:call, gensym(), map(Symbol,vars)...),
                  body))
    catch err
        throw(ArgumentError("Expression does not lambdify"))
    end
end

# N
"""

Convert a SymEngine numeric value into a number

"""
N(b::Basic) = N(BasicType(b))
N(b::BasicType{Val{:Integer}}) = eval(parse(toString(b)))  ## HACKY
N(b::BasicType{Val{:Rational}}) = eval(parse(replace(toString(b), "/", "//")))
N(b::BasicType{Val{:Complex}}) =  eval(parse(replace(toString(b), "I", "im")))

## function N(b::BasicType{Val{:Rational}})
##     println("XXX")
## end
## need to test for free_symbols, if none then we need to evaluate
function N(b::BasicType)
    b = convert(Basic, b)
    fs = free_symbols(b)
    if length(fs) > 0
        throw(ArgumentError("Object can have no free symbols"))
    end
    eval(lambdify(b))
end

N(a::Integer) = a
N(a::Rational) = a
N(a::Complex) = a


text:
struct TreeData
    last_node_size::Int    # Number of points in the last node
    leafsize::Int          # Number of points in each leaf node (except last)
    n_leafs::Int           # Number of leafs
    n_internal_nodes::Int  # Number of non leaf nodes
    cross_node::Int
    offset::Int
    offset_cross::Int
    last_full_node::Int
end

function TreeData(data::AbstractVector{V}, leafsize) where V
    n_dim, n_p = length(V), length(data)

    # If number of points is zero
    n_p == 0 && return TreeData(0, 0, 0, 0, 0, 0, 0, 0)

    n_leafs =  ceil(Integer, n_p / leafsize)
    n_internal_nodes = n_leafs - 1
    leafrow = floor(Integer, log2(n_leafs))
    cross_node = 2^(leafrow + 1)
    last_node_size = n_p % leafsize
    if last_node_size == 0
        last_node_size = leafsize
    end

    # This only happens when n_p / leafsize is a power of 2?
    if cross_node >= n_internal_nodes + n_leafs
        cross_node = div(cross_node, 2)
    end

    offset = 2(n_leafs - 2^leafrow) - 1
    k1 = (offset - n_internal_nodes - 1) * leafsize + last_node_size + 1
    k2 = -cross_node * leafsize + 1
    last_full_node = n_leafs + n_internal_nodes

    TreeData(last_node_size, leafsize, n_leafs,
    n_internal_nodes, cross_node, k1, k2, last_full_node)
end


text:
using NumPyArrays
using PyCall
using Test

@testset "NumPyArrays.jl" begin
    np = pyimport("numpy")
    let A = Float64[1 2; 3 4]
        # Normal array
        B = copy(A)
        C = NumPyArray(B)
        D = PyArray(C)
        @test pytypeof(C) == np.ndarray
        @test C == B
        @test D == B
        B[1] = 3
        @test C == B && C[1] == B[1]
        @test D == B && D[1] == B[1]
        C[1] += 1
        @test C == B && C[1] == B[1]
        @test D == B && D[1] == B[1]

        # SubArray
        B = view(A, 1:2, 2:2)
        C = NumPyArray(B)
        D = PyArray(C)
        @test pytypeof(C) == np.ndarray
        @test C == B
        @test D == B
        A[3] = 5
        @test C == B && C[1] == A[3]
        @test D == B && D[1] == A[3]
        C[1] += 1
        @test C == B && C[1] == B[1]
        @test D == B && D[1] == B[1]

        # ReshapedArray
        B = Base.ReshapedArray( A, (1,4), () )
        C = NumPyArray(B)
        D = PyArray(C)
        @test pytypeof(C) == np.ndarray
        @test C == B
        @test D == B
        A[2] = 6
        @test C == B && C[2] == A[2]
        @test D == B && D[2] == A[2]
        C[1] += 1
        @test C == B && C[1] == B[1]
        @test D == B && D[1] == B[1]

        # PermutedDimsArray
        B = PermutedDimsArray(A, (2,1) )
        C = NumPyArray(B)
        D = PyArray(C)
        @test pytypeof(C) == np.ndarray
        @test C == B
        @test D == B
        A[1] == 7
        @test C == B && C[1] == A[1]
        @test D == B && D[1] == A[1]
        C[1] += 1
        @test C == B && C[1] == B[1]
        @test D == B && D[1] == B[1]

        # ReinterpretArray
        B = reinterpret(UInt64, A)
        C = NumPyArray(B)
        D = PyArray(C)
        @test pytypeof(C) == np.ndarray
        @test C == B
        @test D == B
        A[1] = 12
        @test C == B && C[1] == reinterpret(UInt64, A[1])
        @test D == B && D[1] == reinterpret(UInt64, A[1])
        C[1] += 1
        @test C == B && C[1] == B[1]
        @test D == B && D[1] == B[1]

        # Test display
        rA = reinterpret(UInt8, zeros(Int8, 4, 4))
        nprA = NumPyArray(rA)
        io = IOBuffer()
        @test pytypeof(nprA) == np.ndarray
        @test show(io, nprA)  === nothing
        @test String(take!(io)) == "UInt8[0x00 0x00 0x00 0x00; 0x00 0x00 0x00 0x00; 0x00 0x00 0x00 0x00; 0x00 0x00 0x00 0x00]"

        sA = @view collect(1:16)[5:9]
        npsA = NumPyArray(sA)
        @test pytypeof(npsA) == np.ndarray
        @test show(io, npsA)  === nothing
        @test String(take!(io)) == "[5, 6, 7, 8, 9]"
        println()

        # Test roundtrip
        @test nprA == NumPyArray(nprA)
        @test npsA == NumPyArray(npsA)

        # Test operations
        @test sum(nprA) == np.sum(nprA)
        @test sum(npsA) == np.sum(npsA)
        @test length(nprA) == np.size(nprA)
        @test length(npsA) == np.size(npsA)
    end
end


text:
# module SDPhierarchy

# using MathProgComplex
# using DataStructures
using OPFInstances

export RelaxationContext, Moment, MomentMatrix, SDPDual, SDPPrimal
export SDP_Instance, SDP_Block, SDP_Moment, SDP_Problem

###############################################################################
## Relaxation context, symmetries and cliques
###############################################################################
mutable struct RelaxationContext
    ismultiordered::Bool
    issparse::Bool
    symmetries::Set{DataType}       # ::SortedSet{DataType}
    hierarchykind::Symbol           # :Complex or :Real
    renamevars::Bool                # Replace variables with by shorter named ones
    di::Dict{String, Int}
    ki::Dict{String, Int}
    cstrtypes::Dict{String, Symbol}
    relaxparams::OrderedDict{Symbol, Any}

    RelaxationContext() = new(false,
                              false,
                              Set{DataType}(),
                              :Real,
                              false,
                              Dict{String, Int}(),
                              Dict{String, Int}(),
                              Dict{String, Symbol}(),
                              get_defaultparams())
end

include(joinpath("core", "build_relaxationcontext.jl"))

# include("build_relctx.jl")
include(joinpath("core", "build_decomposition.jl"))

abstract type AbstractSymmetry end
type PhaseInvariance <: AbstractSymmetry end
include("symmetries.jl")

###############################################################################
## Moment Problem
###############################################################################
struct Moment
    conj_part::Exponent
    expl_part::Exponent
    clique::String
end

include(joinpath("base_types", "moment.jl"))

"""
    MomentMatrix{T}(mm, vars, order, matrixkind)

    Store a moment or localizing matrix of size `order`, corresponding to the `vars` variables in the `mm` dictionnary.
    **Note** that the matrix is indexed by a tuple of exponents, *the first of which contains only conjugated variables*, et second only real ones.
"""
mutable struct MomentMatrix{T}
    mm::Dict{Tuple{Exponent, Exponent}, Dict{Moment, T}}
    vars::Set{Variable}
    order::Int
    matrixkind::Symbol            # Either :SDP or :Sym
end

include(joinpath("base_types", "momentmatrix.jl"))

"""
    momentrel = SDPDual(obj, cstrs, moment_overlap)

    Store a Moment Relaxation problem.
"""
struct SDPDual{T}
    objective::Dict{Moment, T}                                  # A linear comb. of moments, to be maximized
    constraints::Dict{Tuple{String, String}, MomentMatrix{T}}   # A set of moment matrices, either SDP or Null. A constraint (`key[1]`) can be split on several cliques (`key[2]`)
    moments_overlap::Dict{Exponent, Set{String}}                # A set of clique per exponent, describing coupling constraints
end

include(joinpath("core", "build_momentrelaxation.jl"))

###############################################################################
## SOS Problem
###############################################################################
mutable struct SDPPrimal{T}
    block_to_vartype::Dict{String, Symbol}                       # Either :SDP, :Sym, :SDPc, :SymC
    blocks::Dict{Tuple{Moment, String, Exponent, Exponent}, T}   # ((α, β), block_name, γ, δ) -> coeff
    linsym::Dict{Tuple{Moment, String, Exponent}, T}             # ((α, β), block_name, var) -> coeff
    lin::Dict{Tuple{Moment, Exponent}, T}                        # ((α, β), var) -> coeff
    cst::Dict{Moment, T}                                         # (α, β) -> coeff
end

include(joinpath("core", "build_SOSrelaxation.jl"))
include(joinpath("io", "export_SDPPrimal.jl"))
include("SDPPrimal_cplx2real.jl")

###############################################################################
## Mosek Structures
###############################################################################
type SDP_Instance
  VAR_TYPES
  BLOCKS
  LINEAR
  CONST
end

type SDP_Block
  id::Int64
  name::String
  var_to_id::SortedDict{String, Int64}

  SDP_Block(id::Int64, name::String) = new(id, name, SortedDict{String, Int64}())
end

const SDP_Moment = Tuple{String, String, String}

type SDP_Problem
  # SDP vars
  name_to_sdpblock::SortedDict{String, SDP_Block}
  id_to_sdpblock::SortedDict{Int64, SDP_Block}

  # Scalar variables
  scalvar_to_id::Dict{String, Int64}

  # Objective / constraints
  obj_keys::SortedSet{SDP_Moment}
  name_to_ctr::SortedDict{SDP_Moment, Tuple{Int64, String, Float64, Float64}} # Id, type et bornes des contraintes
  id_to_ctr::SortedDict{Int64, SDP_Moment}

  matrices::SortedDict{Tuple{SDP_Moment, String, String, String}, Float64} # Matrices SDP du corps des contraintes / objectif
  linear::SortedDict{Tuple{SDP_Moment, String}, Float64} # Matrice portant les parties linéaires des contraintes
  cst_ctr::SortedDict{SDP_Moment, Float64} # Constante du corps des contraintes

  SDP_Problem() = new(SortedDict{String, SDP_Block}(),
                      SortedDict{Int64, SDP_Block}(),
                      Dict{String, Int64}(),
                      SortedSet{SDP_Moment}(),
                      SortedDict{SDP_Moment, Tuple{Int64, String, Float64, Float64}}(),
                      SortedDict{Int64, SDP_Moment}(),
                      SortedDict{Tuple{SDP_Moment, String, String, String}, Float64}(),
                      SortedDict{Tuple{SDP_Moment, String}, Float64}(),
                      SortedDict{SDP_Moment, Float64}()
                      )
end

include(joinpath("Mosek", "run_mosek.jl"))
include(joinpath("io", "build_SDP_Problem.jl"))

###############################################################################
## Unsorted
###############################################################################

include(joinpath("core", "run_hierarchy.jl"))

include("example_problems.jl")
include("utils.jl")
include(joinpath("io", "momentsos_io.jl"))

# end


text:
using StanModels, Distributions

# Define the Stan language model

m8_5s = "
data{
  int N;
  vector[N] y;
}
parameters{
  real sigma;
  real a1;
  real a2;
}
model{
  real mu;
  a1 ~ normal(0, 10);
  a2 ~ normal(0, 10);
  mu = a1 + a2;
  y ~ normal( mu , sigma );
}
";

# Define the Stanmodel and set the output format to :mcmcchains.

sm = SampleModel("m8.5s", m8_5s);

# Input data for cmdstan

m8_5_data = Dict("N" => 100, "y" => rand(Normal(0, 1), 100));

# Sample using cmdstan

(sample_file, log_file) = stan_sample(sm, data=m8_5_data);

rethinking = "
       mean   sd   5.5% 94.5% n_eff Rhat
a1    -0.08 7.15 -11.34 11.25  1680    1
a2    -0.05 7.15 -11.37 11.19  1682    1
sigma  0.90 0.07   0.81  1.02  2186    1
";

# Describe the draws
if !(sample_file == nothing)
  chn = read_samples(sm)
  describe(chn)
end


text:
@testset "Husert" begin
  a = fmpz[0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 87, -64, 15, -80, 0, 0,
           0, 0, 0, 49]
  b = fmpz[135, -2036, 2457, -243, -22300, -45, 679, -819, 81, 7434, 1658,
           -6743, 23115, -2294, -189610, 14922, -42584, 200191, -19879,
           -1616804, 25, -350, 455, -45, -4081]

  A = matrix(ZZ, 5, 5, a)
  B = matrix(ZZ, 5, 5, b)
  fl, C = Hecke._issimilar_husert(A, B)
  @test fl
  @test C * A == B * C
end


text:
module NonLinearOperatorMocks

export NonLinearOperatorMock

using Gridap
import Gridap: residual!, jacobian, jacobian!
import Gridap: create_in_domain

struct NonLinearOperatorMock <: NonLinearOperator end

function residual!(b::AbstractVector,::NonLinearOperatorMock,x::AbstractVector)
  b[1] = (x[1]-2)*(x[1]-1)
  b[2] = x[2]-3
end

function jacobian!(A::AbstractMatrix,::NonLinearOperatorMock,x::AbstractVector)
  A[1,1] = (x[1]-1) + (x[1]-2)
  A[1,2] = 0
  A[2,1] = 0
  A[2,2] = 1
end

function jacobian(op::NonLinearOperatorMock,x::AbstractVector)
  A = zeros(2,2)
  jacobian!(A,op,x)
  A
end

function create_in_domain(::NonLinearOperatorMock)
  zeros(2)
end

end # module


text:
abstract type KValueModel <: ThermoModel end

function kvalues(model::KValueModel,st::ThermodynamicState)
    return kvalues(state_type(st),model,st)
end

abstract type RachfordRiceSolver end
struct RachfordRice <: RachfordRiceSolver end
struct MultiPhaseRachfordRice <: RachfordRiceSolver end
struct PolynomialRachfordRice <: RachfordRiceSolver end
struct LeiboviciNichita  <: RachfordRiceSolver end
struct LiJohnsAhmadi <: RachfordRiceSolver end

struct AnalyticalRachfordRiceSolution{N} <: RachfordRiceSolver end

const RR = RachfordRice
const MultiRR = MultiPhaseRachfordRice
const PolynomialRR = PolynomialRachfordRice
const LN2 = LeiboviciNichita
const LJA = LiJohnsAhmadi

function vf0(K,z)
    N = length(K)
    Kmin = minimum(K)
    Kmax,i_kmax = findmax(K)
    z_of_Kmax = z[i_kmax]
    Kminm1 = one(Kmin)-Kmin
    V_over_F_min = ((Kmax-Kmin)*z_of_Kmax - Kminm1)/(Kminm1*(Kmax - one(Kmin)))
    V_over_F_max = one(Kmin)/Kminm1
    return (V_over_F_min,V_over_F_max)
end

function rr_find_zero(f0,(vmin,vmax),method::Roots.AbstractBracketing)
    return Roots.find_zero(f0,(vmin,vmax),method)
end

function rr_find_zero(f0,(vmin,vmax),method::Roots.AbstractNonBracketing)
    β0 = (vmin+vmax)/2
    return Roots.find_zero(f0,β0,method)
end

function rr_find_zero(f0,(vmin,vmax),method::Roots.Newton)
    β0 = (vmin+vmax)/2
    AD_f0(x) = autonewton(f0,x)
    return Roots.find_zero(f0,β0,method)
end

function flash_vfrac(model::RR,K,Z,roots_method=Roots.Order0())
    #@show K
    Kmin0 = minimum(K)
    Kmax0 = maximum(K)
    #if (Kmin0 > 1.0) | (Kmax0 < 1.0)
    #    return ArgumentError("no positive composition solution with given K-values")
    #end
    vfmin,vfmax = vf0(K,Z)
    abs(vfmin - vfmax) < sqrt(eps(1.0)) && return 0.5*(vfmin + vfmax)
    #vfmin = max(vfmin,zero(vfmin))
    #vfmax = min(vfmax,one(vfmax))
    #vfmin > vfmax && return 1.0
    #@show vfmin,vfmax
    function rr(k,z,β)
        if iszero(β)
            rr0(k,z)
        elseif isone(β)
            rr1(k,z)
        else
            rry(k,z,β)
        end
    end
    f0(β) = mapreduce((k,z)->rr(k,z,β),+,K,Z)
    βres = rr_find_zero(f0,(vfmin,vfmax),roots_method)
    return βres
end

function flash_vfrac(model::LeiboviciNichita,K,Z,roots_method=Roots.Newton())
    Kmin0 = minimum(K)
    Kmax0 = maximum(K)
    if (Kmin0 > 1.0) | (Kmax0 < 1.0)
        return ArgumentError("no positive composition solution with given K-values")
    end
    vfmin,vfmax = vf0(K,Z)
    ck = vfmin
    ckm1 = vfmax
    _λ(y) = ck + (ckm1-ck)/(one(y)+exp(-y))
    _y(λ) = log((ckm1-ck)/(λ-ck)-one(λ))
    function rr(k,z,y)
        c = 1/(1-k)
        return z/_y((y)-c)
    end
    println((_y(vfmin),_y(vfmax)))
    f0(β) = mapreduce((k,z)->rr(k,z,β),+,K,Z)
    yres = rr_find_zero(f0,(_y(vfmin),_y(vfmax)),roots_method)
    return _λ(yres)
end

function flash_vfrac(model::PolynomialRachfordRice,K,Z,method=nothing)
    n = length(Z)
    if 2 <= n <= 4
        return flash_vfrac(AnalyticalRachfordRiceSolution{n}(),K,Z)
    else
        return DomainError("analytical solutions are only available for 2,3,or 4 compounds")
    end
end
#=
n0 = b*nv + (1-b)*nl
nv + nl = n0
z0 = nv/(nv+nl) + (1-b)nl/(nv+nl)

=#

function rr0(k,z)
    km1 = k-one(k)
    return z*km1
end

function rr1(k::T,z) where T
    _1 = one(T)
    return z*(_1-_1/k)
end

function rry(k::T,z,y) where T
    _1 = one(T)
    km1 = k-_1
    return z*km1/(_1+y*km1)
end

function flash_eval(K,Z,β)
    if iszero(β)
        return mapreduce((k,z)->rr0(k,z),+,K,Z)
    elseif isone(β)
        return mapreduce((k,z)->rr1(k,z),+,K,Z)
    else
        return mapreduce((k,z)->rry(k,z,β),+,K,Z)
    end
end
###K.value

##kᵢ = xvᵢ/xlᵢ

# if kᵢ = 0, there is not liquid phase (non-condensable)
# if kᵢ = inf there is not gas phase (only liquid)
#kᵢ = 0, ignore in liquid phase
#kᵢ = Inf or NaN, ignore in liquid phase

"""
    flash_vapor(k,z,β)
Returns the gas phase composition, given k-values `k`, the initial molar composition `z` and the molar vapor fraction ´β´.

Each gas phase composition is calculated acording to:

    xvᵢ = kᵢzᵢ/(1+ β(kᵢ-1))

For an implace-if-posible version, use `flash_vapor!!(res,k,z,β)`

"""
function flash_vapor(k,z,β)
    function f(ki,zi)
    _1 = one(ki)
        return ki*zi/(_1+β*(ki-_1))
    end

    return map(f,k,z)
end
"""
    flash_liquid(k,z,β)
Returns the liquid phase composition, given k-values `k`, the initial molar composition `z` and the molar vapor fraction ´β´.

Each gas phase composition is calculated acording to:

    xlᵢ = zᵢ/(1+ β(kᵢ-1))

For an implace-if-posible version, use `flash_liquid!!(res,k,z,β)`
"""
function flash_liquid(k,z,β)
    function f(ki,zi)
        _1 = one(ki)
        return zi/(_1+β*(ki-_1))
    end
    return map(f,k,z)
end

function flash_vapor!!(res,k,z,β)
    _1 = one(eltype(res))
    @! res +=  k .* z ./ (_1 .+ β .*(k .-_1)) .- res

    return res
end

function flash_liquid!!(res,k,z,β)
    _1 = one(eltype(k))
    @! res += z ./ (_1 .+ β .*(k .-_1)) .- res
    return res
end

struct WilsonK{T} <: ThermoModel
    pc::T
    tc::T
    ω::T
end

struct MollerupK{T} <: ThermoModel
    pc::T
    tc::T
    ω::T
end

WilsonK(;pc,tc,ω) = WilsonK(pc,tc,ω)
MollerupK(;pc,tc,ω) = MollerupK(pc,tc,ω)
temperature(model::WilsonK,st::CriticalPoint,unit = u"K") = convert_unit(u"K",unit,model.tc)
temperature(model::MollerupK,st::CriticalPoint,unit = u"K") = convert_unit(u"K",unit,model.tc)
pressure(model::WilsonK,st::CriticalPoint,unit = u"Pa") = convert_unit(u"Pa",unit,model.pc)
pressure(model::MollerupK,st::CriticalPoint,unit = u"Pa") = convert_unit(u"Pa",unit,model.pc)
acentric_factor(model::WilsonK) = model.ω
acentric_factor(model::MollerupK) = model.ω

function WilsonK(model::ThermoModel)
    _tc = temperature(model,CriticalPoint(),u"K")
    _pc = pressure(model,CriticalPoint(),u"Pa")
    _ω = acentric_factor(model)
    return WilsonK(_pc,_tc,_ω)
end

function MollerupK(model::ThermoModel)
    _tc = temperature(model,CriticalPoint(),u"K")
    _pc = pressure(model,CriticalPoint(),u"Pa")
    _ω = acentric_factor(model)
    return MollerupK(_pc,_tc,_ω)
end

function kvalues(mt::MultiPT,model::WilsonK,st::ThermodynamicState)
    p = pressure(FromState(),st)
    t = temperature(FromState(),st)
    return kvalues_impl(mt,model,p,t)
end

function kvalues_impl(mt::MultiPT,model::WilsonK,p,t,x)
    pc =model.pc
    tc = model.tc
    ω = model.ω
    k0 = copy(x)
    k0-= k0
    k0 += exp.(log.(pc./p).+5.373 .*(1.0 .+ ω).*(1.0 .-tc./t))
end

function kvalues(mt::MultiPT,model::MollerupK,st::ThermodynamicState)
    p = pressure(FromState(),st)
    t = temperature(FromState(),st)
    x = mol_fraction(FromState(),st)
    return kvalues_impl(mt,model,p,t,x)
end

function kvalues_impl(mt::MultiPT,model::MollerupK,p,t,x)
    pc = pressure(model,ThermoState.CriticalPoint())
    tc = temperature(model,ThermoState.CriticalPoint())
    ω = acentric_factor(model)
    k0 = copy(x)
    k0 -= k0
    k0 +=   (pc ./ p) .* exp.(5.42 .* (1.0 .- (tc ./ t)))
end


text:
module YaoLocations

export CtrlLocations, Locations, AbstractLocations, LocationError, merge_locations

export Locations, CtrlLocations, AbstractLocations, LocationError, merge_locations

const LocationStorageTypes = Union{Int,NTuple{N,Int} where N,UnitRange{Int},Vector{Int}}

abstract type AbstractLocations end

"""
    plain(locations)

Return plain old location types that represents locations as integer ids,
such as `Int`, `NTuple{N, Int}`, `UnitRange{Int}`.
"""
function plain end

"""
    merge_locations(locations...)

Construct a new `Locations` by merging two or more existing locations.
"""
merge_locations(x::AbstractLocations, y::AbstractLocations, locations::AbstractLocations...) =
    merge_locations(merge_locations(x, y), locations...)

merge_locations(x::AbstractLocations) = x

"""
    Locations <: AbstractLocations

Type to annotate locations in quantum circuit.

    Locations(x)

Create a `Locations` object from a raw location statement. Valid storage types are:

- `Int`: single position
- `NTuple{N, Int}`: a list of locations
- `UnitRange{Int}`: contiguous locations

Other types will be converted to the storage type via `Tuple`.
"""
struct Locations{T<:LocationStorageTypes} <: AbstractLocations
    storage::T

    Base.@pure function Locations(x::T) where {T<:LocationStorageTypes}
        new{T}(x)
    end
end

# skip it if x is a location
Locations(x::Locations) = x
Locations(xs::NTuple{N,<:Locations}) where {N} = merge_locations(xs...)

Base.@propagate_inbounds Base.getindex(l::Locations, idx...) = Locations(getindex(l.storage, idx...))
Base.length(l::Locations) = length(l.storage)
Base.eltype(::Type{T}) where {T<:Locations} = Locations{Int}
Base.eltype(::Locations) = Locations{Int}
Base.Tuple(x::Locations) = (x.storage...,)
plain(x::Locations) = x.storage

function Base.iterate(l::Locations, st::Int = 1)
    st > length(l) && return
    return l[st], st + 1
end

struct LocationError <: Exception
    msg::String
end

@noinline function merge_locations(l1::Locations, l2::Locations)
    Locations((l1.storage..., l2.storage...))
end

# location mapping
# TODO: preserve sign when indexing
# TODO: provide a @inlocation macro via Expr(:meta, :inlocation, true) so when we compile to Julia functions
#       we can use unsafe_mapping directly
@noinline unsafe_mapping(parent::Locations{Int}, sub::Locations{Int}) = parent
@noinline unsafe_mapping(parent::Locations{Int}, sub::Locations{NTuple{N,Int}}) where {N} = parent
@noinline unsafe_mapping(parent::Locations{Int}, sub::Locations{UnitRange{Int}}) = parent
@noinline unsafe_mapping(parent::Locations{NTuple{N,Int}}, sub::Locations{Int}) where {N} =
    Locations(@inbounds parent[sub.storage])
@noinline unsafe_mapping(
    parent::Locations{NTuple{N,Int}},
    sub::Locations{NTuple{M,Int}},
) where {N,M} = Locations(map(x -> @inbounds(parent[x]), sub.storage))
@noinline unsafe_mapping(parent::Locations{NTuple{N,Int}}, sub::Locations{UnitRange{Int}}) where {N} =
    Locations(@inbounds parent[sub.storage])
@noinline unsafe_mapping(parent::Locations{UnitRange{Int}}, sub::Locations{Int}) =
    Locations(@inbounds parent[sub.storage])
@noinline unsafe_mapping(parent::Locations{UnitRange{Int}}, sub::Locations{NTuple{N,Int}}) where {N} =
    Locations(map(x -> @inbounds(parent[x]), sub.storage))
@noinline unsafe_mapping(parent::Locations{UnitRange{Int}}, sub::Locations{UnitRange{Int}}) =
    Locations(@inbounds parent[sub.storage])

@inline map_error(parent, sub) = throw(LocationError("got $sub in parent space $parent"))

@inline function map_check(parent, sub)
    map_check_nothrow(parent, sub) || map_error(parent, sub)
end

@noinline function map_check_nothrow(parent::Locations{Int}, sub::Locations{Int})
    sub.storage == 1
end

@noinline function map_check_nothrow(parent::Locations{Int}, sub::Locations{Tuple{Int}})
    sub.storage[1] == 1
end

@noinline function map_check_nothrow(parent::Locations{Int}, sub::Locations{NTuple{N,Int}}) where {N}
    false
end

@noinline function map_check_nothrow(parent::Locations{Int}, sub::Locations{UnitRange{Int}})
    (length(sub) == 1) && (sub.storage.start == 1)
end

@noinline function map_check_nothrow(parent::Locations{NTuple{N,Int}}, sub::Locations{Int}) where {N}
    1 <= sub.storage <= N
end

@noinline function map_check_nothrow(
    parent::Locations{NTuple{N,Int}},
    sub::Locations{NTuple{M,Int}},
) where {N,M}
    all(x -> (1 <= x <= N), sub.storage)
end

@noinline function map_check_nothrow(
    parent::Locations{NTuple{N,Int}},
    sub::Locations{UnitRange{Int}},
) where {N}
    (1 <= sub.storage.start) && (sub.storage.stop <= N)
end

@noinline function map_check_nothrow(parent::Locations{UnitRange{Int}}, sub::Locations{Int})
    1 <= sub.storage <= length(parent)
end

@noinline function map_check_nothrow(
    parent::Locations{UnitRange{Int}},
    sub::Locations{NTuple{N,Int}},
) where {N}
    all(x -> (1 <= x <= length(parent)), sub.storage)
end

@noinline function map_check_nothrow(
    parent::Locations{UnitRange{Int}},
    sub::Locations{UnitRange{Int}},
)
    (1 <= sub.storage.start) && (sub.storage.stop <= length(parent))
end

# comparing
function Base.:(==)(l1::Locations, l2::Locations)
    length(l1) == length(l2) || return false
    return all(l1.storage .== l2.storage)
end

struct CtrlFlags{L,N}
    data::NTuple{N,UInt64}
end

function flags(bits::Vararg{UInt8,L}) where {L}
    N = Base.num_bit_chunks(L)
    data = ntuple(N) do k
        x = UInt64(0)
        p = UInt64(0)
        for i in (64k-63):min(64k, L)
            x += bits[i] << p
            p += UInt64(1)
        end
        return x
    end
    return CtrlFlags{L,N}(data)
end

function _default_flags(L::Int)
    N = Base.num_bit_chunks(L)
    return CtrlFlags{L,N}(ntuple(x -> UInt64(0), N))
end

function Base.all(flags::CtrlFlags)
    return all(iszero, flags.data)
end

function Base.getindex(flags::CtrlFlags, idx::Int)
    i1, i2 = Base.get_chunks_id(idx)
    u = UInt64(1) << i2
    r = (flags.data[i1] & u) == 0
    return r
end

function merge_flags(a::CtrlFlags{LA}, b::CtrlFlags{LB}) where {LA,LB}
    L = LA + LB
    N = Base.num_bit_chunks(L)
    last_chunk_len = rem(LA, 64)

    if iszero(last_chunk_len)
        return CtrlFlags{L,N}((a.data..., b.data...))
    else
        return flags(ntuple(k -> UInt8(!a[k]), LA)..., ntuple(k -> UInt8(!b[k]), LB)...)
    end
end

function Base.show(io::IO, x::CtrlFlags{L,N}) where {L,N}
    print(io, "CtrlFlags(")
    p = L - 1
    for d in x.data
        for k in 0:min(63, p)
            print(io, d >> k & UInt64(1))
        end
        p -= 64
    end
    print(io, ")")
end

# CtrlLocations
struct CtrlLocations{T<:LocationStorageTypes,L,N} <: AbstractLocations
    storage::Locations{T}
    flags::CtrlFlags{L,N}

    Base.@pure CtrlLocations(storage::Locations{T}, flags::CtrlFlags{L,N}) where {T,L,N} =
        new{T,L,N}(storage, flags)
end

# skip itself
CtrlLocations(x::CtrlLocations) = x
CtrlLocations(x::Locations) = CtrlLocations(x, _default_flags(length(x)))
CtrlLocations(x::LocationStorageTypes) = CtrlLocations(Locations(x), _default_flags(length(x)))
CtrlLocations(x::LocationStorageTypes, configs::NTuple{L,UInt8}) where {L} =
    CtrlLocations(Locations(x), flags(configs...))

Base.length(l::CtrlLocations) = length(l.storage)

function Base.iterate(l::CtrlLocations, st::Int = 1)
    st > length(l) && return
    return (l.storage[st], l.flags[st]), st + 1
end

plain(l::CtrlLocations) = plain(l.storage)

function Base.show(io::IO, x::Locations)
    print(io, "Locations(")
    print_locations(io, x)
    print(io, ")")
end

function Base.show(io::IO, x::CtrlLocations)
    print(io, "CtrlLocations(")
    print_locations(io, x)
    print(io, ")")
end

function print_locations(io::IO, x::Locations)
    if x.storage isa AbstractRange
        return printstyled(io, x.storage; color = :light_blue)
    end

    nlocations = length(x)
    for i in 1:nlocations
        printstyled(io, x.storage[i]; color = :light_blue)

        if i != nlocations
            print(io, ", ")
        end
    end
end

function print_locations(io::IO, x::CtrlLocations)
    if all(x.flags)
        print_locations(io, x.storage)
    else
        nlocations = length(x)
        for i in 1:nlocations
            l = x.storage.storage[i]
            if x.flags[i]
                printstyled(io, l; color = :light_blue)
            else
                printstyled(io, "!", l; color = :light_blue)
            end

            if i != nlocations
                print(io, ", ")
            end
        end
    end
end

@noinline function merge_locations(l1::CtrlLocations, l2::CtrlLocations)
    CtrlLocations(merge_locations(l1.storage, l2.storage), merge_flags(l1.flags, l2.flags))
end

function merge_locations(l1::CtrlLocations, l2::Locations)
    throw(
        LocationError(
            "cannot merge CtrlLocations and Locations, convert to CtrlLocations to Locations",
        ),
    )
end

function merge_locations(l1::Locations, l2::CtrlLocations)
    throw(
        LocationError(
            "cannot merge CtrlLocations and Locations, convert to CtrlLocations to Locations",
        ),
    )
end

# NOTE: CtrlLocations can not be mapped by Locations
@noinline unsafe_mapping(parent::Locations, sub::CtrlLocations) =
    CtrlLocations(unsafe_mapping(parent, sub.storage), sub.flags)
@noinline map_check_nothrow(parent::Locations, sub::CtrlLocations) =
    map_check_nothrow(parent, sub.storage)

Base.:(==)(l1::CtrlLocations, l2::CtrlLocations) =
    (l1.storage == l2.storage) && (l1.flags == l2.flags)

# parent location has to be Locations since CtrlLocations can't be parent
@inline function Base.getindex(parent::Locations, sub::AbstractLocations)
    map_check(parent, sub)
    return unsafe_mapping(parent, sub)
end

Base.issorted(locs::AbstractLocations) = issorted(plain(locs))

Base.sort(locs::Locations{Int}) = locs
Base.sort(locs::Locations{UnitRange{Int}}) = locs

function Base.sort(locs::Locations{<:Tuple})
    return Locations(sort!(collect(locs.storage)))
end

function Base.sort(locs::Locations{Vector{Int}})
    return Locations(sort!(locs.storage))
end

end # module


text:
module NameList
# See Table 2 of Cohen et al, 2020
# namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_factor"] = 0.13
# namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["detrainment_factor"] = 0.51
# namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_massflux_div_factor"] = 0.0
# namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["turbulent_entrainment_factor"] = 0.015
# namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_smin_tke_coeff"] = 0.3
# namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["updraft_mixing_frac"] = 0.25
# namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_sigma"] = 10.0
# namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_scale"] = 0.004
# namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["sorting_power"] = 2.0
# namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["aspect_ratio"] = 0.2

# See Table 1 of Lopez Gomez et al, 2020
# namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["tke_ed_coeff"] = 0.14
# namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["tke_diss_coeff"] = 0.22
# namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["static_stab_coeff"] = 0.4 # Square of value in the paper

# See Table ? of He et al, 2021
# namelist["turbulence"]["EDMF_PrognosticTKE"]["pressure_normalmode_buoy_coeff1"] = 0.12 ==> alpha_b (scaling constant for virtual mass term)
# namelist["turbulence"]["EDMF_PrognosticTKE"]["pressure_normalmode_buoy_coeff2"] = 0.0 ==> alpha_b (scaling constant for virtual mass term)
# namelist["turbulence"]["EDMF_PrognosticTKE"]["pressure_normalmode_adv_coeff"] = 0.1 alpha_a (scaling constant for advection term)
# namelist["turbulence"]["EDMF_PrognosticTKE"]["pressure_normalmode_drag_coeff"] = 10.0 ==> alpha_d (scaling constant for drag term)
# namelist["turbulence"]["EDMF_PrognosticTKE"]["pressure_plume_spacing"] ==> r_d (horizontal length scale of plume spacing)

#NB: except for Bomex and life_cycle_Tan2018 cases, the parameters listed have not been thoroughly tuned/tested
# and should be regarded as placeholders only. Optimal parameters may also depend on namelist options, such as
# entrainment/detrainment rate formulation, diagnostic vs. prognostic updrafts, and vertical resolution
export default_namelist

using ArgParse
using ArtifactWrappers

import StaticArrays
const SA = StaticArrays

#! format: off
LESDrivenSCM_output_dataset = ArtifactWrapper(
    @__DIR__,
    isempty(get(ENV, "CI", "")),
    "LESDrivenSCM_output_dataset",
    ArtifactFile[
        ArtifactFile(url = "https://caltech.box.com/shared/static/0hnf7nkttueraaqf9tpkqsx38gjqx41p.nc", filename = "Stats.cfsite23_HadGEM2-A_amip_2004-2008.07.nc",),
    ],
)
LESDrivenSCM_output_dataset_path = get_data_folder(LESDrivenSCM_output_dataset)
#! format: on

import JSON

function parse_commandline()
    s = ArgParseSettings(; description = "namelist Generator")

    @add_arg_table! s begin
        "case_name"
        help = "The case name"
        arg_type = String
        required = true
    end

    return parse_args(s)
end

function default_namelist(::Nothing)

    args = parse_commandline()
    case_name = args["case_name"]
    return default_namelist(case_name)
end

function default_namelist(case_name::String; root::String = ".", write::Bool = true)

    namelist_defaults = Dict()
    namelist_defaults["meta"] = Dict()
    namelist_defaults["meta"]["uuid"] = basename(tempname())

    namelist_defaults["turbulence"] = Dict()
    namelist_defaults["turbulence"]["Ri_bulk_crit"] = 0.2
    namelist_defaults["turbulence"]["prandtl_number_0"] = 0.74

    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"] = Dict()
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["surface_area"] = 0.1
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["max_area"] = 0.9
    # mixing_length
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["tke_ed_coeff"] = 0.14
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["tke_diss_coeff"] = 0.22
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["static_stab_coeff"] = 0.4
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["Prandtl_number_scale"] = 53.0 / 13.0
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["smin_ub"] = 0.1
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["smin_rm"] = 1.5
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["l_max"] = 1.0e6
    # entrainment
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_factor"] = 0.13
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["detrainment_factor"] = 0.51
    # 1-layer nn parameters
    #! format: off
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["general_ent_params"] =
        SA.SVector(0.3038, 0.719,-0.910,-0.483,
                   0.739, 0.0755, 0.178, 0.521,
                   0.0, 0.0, 0.843,-0.340,
                   0.655, 0.113, 0.0, 0.0)
    #! format: on

    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_massflux_div_factor"] = 0.0
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["turbulent_entrainment_factor"] = 0.015
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_smin_tke_coeff"] = 0.3
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["updraft_mixing_frac"] = 0.25
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_sigma"] = 10.0
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["area_limiter_scale"] = 10.0
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["area_limiter_power"] = 3.0
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment_scale"] = 0.004
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["sorting_power"] = 2.0
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["min_upd_velocity"] = 0.001
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["aspect_ratio"] = 0.2
    # pressure
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["min_updraft_top"] = 500.0
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["pressure_normalmode_buoy_coeff1"] = 0.12
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["pressure_normalmode_buoy_coeff2"] = 0.0
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["pressure_normalmode_adv_coeff"] = 0.1
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["pressure_normalmode_drag_coeff"] = 10.0

    # stochastic closures
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["stochastic"] = Dict()
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["stochastic"]["closure"] = "none"
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["stochastic"]["entr_lognormal_var"] = 0.0
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["stochastic"]["detr_lognormal_var"] = 0.0
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["stochastic"]["sde_entr_theta"] = 1.0
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["stochastic"]["sde_entr_std"] = 0.0
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["stochastic"]["sde_detr_theta"] = 1.0
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["stochastic"]["sde_detr_std"] = 0.0

    # From namelist
    namelist_defaults["grid"] = Dict()
    namelist_defaults["grid"]["dims"] = 1

    namelist_defaults["thermodynamics"] = Dict()
    namelist_defaults["thermodynamics"]["thermal_variable"] = "thetal"
    namelist_defaults["thermodynamics"]["sgs"] = "quadrature"
    namelist_defaults["thermodynamics"]["quadrature_order"] = 3
    namelist_defaults["thermodynamics"]["quadrature_type"] = "log-normal" #"gaussian" or "log-normal"

    namelist_defaults["time_stepping"] = Dict()
    namelist_defaults["time_stepping"]["dt_max"] = 12.0
    namelist_defaults["time_stepping"]["dt_min"] = 1.0
    namelist_defaults["time_stepping"]["adapt_dt"] = true
    namelist_defaults["time_stepping"]["cfl_limit"] = 0.5

    namelist_defaults["microphysics"] = Dict()
    namelist_defaults["microphysics"]["precipitation_model"] = "None"

    namelist_defaults["turbulence"]["scheme"] = "EDMF_PrognosticTKE"

    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["updraft_number"] = 1
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["entrainment"] = "moisture_deficit" # "moisture_deficit" or "NN"
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["extrapolate_buoyancy"] = true
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["use_local_micro"] = true
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["constant_area"] = false
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["calculate_tke"] = true
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["mixing_length"] = "sbtd_eq"
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["env_buoy_grad"] = "quadratures"

    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["pressure_closure_buoy"] = "normalmode"
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["pressure_closure_drag"] = "normalmode"
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["pressure_closure_asp_label"] = "const"

    namelist_defaults["output"] = Dict()
    namelist_defaults["output"]["output_root"] = "./"

    namelist_defaults["stats_io"] = Dict()
    namelist_defaults["stats_io"]["stats_dir"] = "stats"
    namelist_defaults["stats_io"]["frequency"] = 60.0
    namelist_defaults["stats_io"]["skip"] = false

    if case_name == "Soares"
        namelist = Soares(namelist_defaults)
    elseif case_name == "Nieuwstadt"
        namelist = Nieuwstadt(namelist_defaults)
    elseif case_name == "Bomex"
        namelist = Bomex(namelist_defaults)
    elseif case_name == "life_cycle_Tan2018"
        namelist = life_cycle_Tan2018(namelist_defaults)
    elseif case_name == "Rico"
        namelist = Rico(namelist_defaults)
    elseif case_name == "TRMM_LBA"
        namelist = TRMM_LBA(namelist_defaults)
    elseif case_name == "ARM_SGP"
        namelist = ARM_SGP(namelist_defaults)
    elseif case_name == "GATE_III"
        namelist = GATE_III(namelist_defaults)
    elseif case_name == "DYCOMS_RF01"
        namelist = DYCOMS_RF01(namelist_defaults)
    elseif case_name == "GABLS"
        namelist = GABLS(namelist_defaults)
    elseif case_name == "SP"
        namelist = SP(namelist_defaults)
    elseif case_name == "DryBubble"
        namelist = DryBubble(namelist_defaults)
    elseif case_name == "LES_driven_SCM"
        namelist = LES_driven_SCM(namelist_defaults)
    else
        error("Not a valid case name")
    end

    if write
        write_file(namelist, root)
    end
    return namelist
end
function Soares(namelist_defaults)

    namelist = deepcopy(namelist_defaults)
    namelist["meta"]["casename"] = "Soares"

    namelist["grid"]["nz"] = 75
    namelist["grid"]["dz"] = 50.0

    namelist["time_stepping"]["t_max"] = 8 * 3600.0
    namelist["time_stepping"]["dt_min"] = 1.0

    namelist["meta"]["simname"] = "Soares"
    namelist["meta"]["casename"] = "Soares"

    return namelist
end
function Nieuwstadt(namelist_defaults)

    namelist = deepcopy(namelist_defaults)
    namelist["meta"]["casename"] = "Nieuwstadt"

    namelist["grid"]["nz"] = 75
    namelist["grid"]["dz"] = 50.0

    namelist["time_stepping"]["t_max"] = 8 * 3600.0
    namelist["time_stepping"]["dt_min"] = 1.2

    namelist["meta"]["simname"] = "Nieuwstadt"
    namelist["meta"]["casename"] = "Nieuwstadt"

    return namelist
end
function Bomex(namelist_defaults)

    namelist = deepcopy(namelist_defaults)
    namelist["meta"]["casename"] = "Bomex"

    namelist["grid"]["nz"] = 60
    namelist["grid"]["dz"] = 50.0

    namelist["time_stepping"]["t_max"] = 21600.0
    namelist["time_stepping"]["dt_min"] = 6.0

    namelist["meta"]["simname"] = "Bomex"
    namelist["meta"]["casename"] = "Bomex"

    return namelist
end
function life_cycle_Tan2018(namelist_defaults)

    namelist = deepcopy(namelist_defaults)
    namelist["meta"]["casename"] = "life_cycle_Tan2018"

    namelist["grid"]["nz"] = 75
    namelist["grid"]["dz"] = 40.0

    namelist["time_stepping"]["t_max"] = 6 * 3600.0
    namelist["time_stepping"]["dt_min"] = 10.0
    namelist["meta"]["simname"] = "life_cycle_Tan2018"
    namelist["meta"]["casename"] = "life_cycle_Tan2018"

    return namelist
end
function Rico(namelist_defaults)

    namelist = deepcopy(namelist_defaults)

    namelist["meta"]["casename"] = "Rico"

    namelist["grid"]["nz"] = 80
    namelist["grid"]["dz"] = 50.0

    # Rico cannot run with adaptive dt for now
    namelist["time_stepping"]["adapt_dt"] = true
    namelist["time_stepping"]["t_max"] = 86400.0
    namelist["time_stepping"]["dt_max"] = 9.0
    namelist["time_stepping"]["dt_min"] = 1.5

    namelist["microphysics"]["precipitation_model"] = "clima_1m"
    namelist["microphysics"]["τ_acnv_rai"] = 2500.0

    namelist["meta"]["simname"] = "Rico"
    namelist["meta"]["casename"] = "Rico"
    return namelist
end
function TRMM_LBA(namelist_defaults)

    namelist = deepcopy(namelist_defaults)

    namelist["meta"]["casename"] = "TRMM_LBA"

    namelist["grid"]["nz"] = 320
    namelist["grid"]["dz"] = 50.0

    namelist["time_stepping"]["t_max"] = 60 * 60 * 6.0
    namelist["time_stepping"]["dt_min"] = 1.0

    namelist["microphysics"]["precipitation_model"] = "cutoff"
    namelist["microphysics"]["τ_precip"] = 50.0

    namelist["meta"]["simname"] = "TRMM_LBA"
    namelist["meta"]["casename"] = "TRMM_LBA"

    return namelist
end
function ARM_SGP(namelist_defaults)

    namelist = deepcopy(namelist_defaults)
    namelist["meta"]["casename"] = "ARM_SGP"

    namelist["grid"]["nz"] = 88
    namelist["grid"]["dz"] = 50.0

    namelist["time_stepping"]["t_max"] = 3600.0 * 14.5
    namelist["time_stepping"]["dt_min"] = 2.0

    namelist["meta"]["simname"] = "ARM_SGP"
    namelist["meta"]["casename"] = "ARM_SGP"

    return namelist
end
function GATE_III(namelist_defaults)

    namelist = deepcopy(namelist_defaults)
    namelist["meta"]["casename"] = "GATE_III"

    namelist["grid"]["nz"] = 1700
    namelist["grid"]["dz"] = 10

    namelist["time_stepping"]["t_max"] = 3600.0 * 24.0
    namelist["time_stepping"]["dt_min"] = 2.0

    namelist["meta"]["simname"] = "GATE_III"
    namelist["meta"]["casename"] = "GATE_III"

    return namelist
end
function DYCOMS_RF01(namelist_defaults)

    namelist = deepcopy(namelist_defaults)

    namelist["meta"]["casename"] = "DYCOMS_RF01"

    namelist["grid"]["nz"] = 30
    namelist["grid"]["dz"] = 50

    namelist["time_stepping"]["t_max"] = 60 * 60 * 16.0
    namelist["time_stepping"]["dt_min"] = 6.0

    namelist["meta"]["simname"] = "DYCOMS_RF01"
    namelist["meta"]["casename"] = "DYCOMS_RF01"

    return namelist
end
function GABLS(namelist_defaults)

    namelist = deepcopy(namelist_defaults)

    namelist["meta"]["casename"] = "GABLS"

    namelist["grid"]["nz"] = 8
    namelist["grid"]["dz"] = 50.0

    namelist["time_stepping"]["t_max"] = 9 * 3600.0
    namelist["time_stepping"]["dt_min"] = 4.0
    namelist["time_stepping"]["dt_max"] = 8.0
    namelist["meta"]["simname"] = "GABLS"
    namelist["meta"]["casename"] = "GABLS"

    return namelist
end
# Not fully implemented yet - Ignacio
function SP(namelist_defaults)

    namelist = deepcopy(namelist_defaults)
    namelist["meta"]["casename"] = "SP"

    # this case is resolution dependent, we should check why
    namelist["grid"]["nz"] = 64
    namelist["grid"]["dz"] = 32

    namelist["time_stepping"]["t_max"] = 7200.0
    namelist["time_stepping"]["dt_min"] = 1.0
    namelist["meta"]["simname"] = "SP"
    namelist["meta"]["casename"] = "SP"

    return namelist
end
function DryBubble(namelist_defaults)
    namelist = deepcopy(namelist_defaults)
    namelist["meta"]["casename"] = "DryBubble"
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["surface_area"] = 0.0

    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["pressure_normalmode_buoy_coeff1"] = 0.12
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["pressure_normalmode_adv_coeff"] = 0.25
    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["pressure_normalmode_drag_coeff"] = 0.1

    namelist["grid"]["nz"] = 200
    namelist["grid"]["dz"] = 50.0

    namelist["stats_io"]["frequency"] = 10.0
    namelist["time_stepping"]["t_max"] = 1000.0
    namelist["time_stepping"]["dt_min"] = 0.5

    namelist["meta"]["simname"] = "DryBubble"
    namelist["meta"]["casename"] = "DryBubble"

    namelist["turbulence"]["EDMF_PrognosticTKE"]["entrainment_massflux_div_factor"] = 0.4

    return namelist
end

function LES_driven_SCM(namelist_defaults)
    namelist = deepcopy(namelist_defaults)
    namelist["grid"]["dz"] = 50.0
    namelist["grid"]["nz"] = 80

    namelist["stats_io"]["frequency"] = 10.0
    namelist["time_stepping"]["t_max"] = 3600.0 * 6
    namelist["time_stepping"]["dt_min"] = 1.0

    # use last 6 hours of LES simulation to drive LES
    namelist["t_interval_from_end_s"] = 3600.0 * 6
    # average in 1 hour interval around `t_interval_from_end_s`
    namelist["initial_condition_averaging_window_s"] = 3600.0

    namelist["meta"]["lesfile"] =
        joinpath(LESDrivenSCM_output_dataset_path, "Stats.cfsite23_HadGEM2-A_amip_2004-2008.07.nc")
    namelist["meta"]["simname"] = "LES_driven_SCM"
    namelist["meta"]["casename"] = "LES_driven_SCM"
    namelist["forcing"] = Dict()
    namelist["forcing"]["nudging_timescale"] = 6.0 * 3600.0

    return namelist
end

function write_file(namelist, root::String = ".")
    mkpath(root)

    @assert haskey(namelist, "meta")
    @assert haskey(namelist["meta"], "simname")

    casename = namelist["meta"]["casename"]
    open(joinpath(root, "namelist_$casename.in"), "w") do io
        JSON.print(io, namelist, 4)
    end

    return
end

if abspath(PROGRAM_FILE) == @__FILE__
    default_namelist(nothing)
end

end


text:
using Submodular

include("test_algorithms.jl")
include("test_combinatorial_sets.jl")
include("test_continuous_sets.jl")
include("test_operators.jl")
include("test_prox.jl")
include("test_submodular_functions.jl")


text:
"""
    A = JopPad(dom, pad..., [extend=false, accumulate=false])

where `dom::JetSpace` is the domain of `A`, and `pad::UnitRange...` determines the range of `A`.
If `extend=false`, then the padded region is set to zero.  If `extend=true`, then the padded
region is set from the boundary of the domain.  The `accumulate=true` option is specific to
the `Ginsu` operation in JetPackWave, and should not be used unless you really know what you
are doing.

# Examples:

## 1D
```julia
A = JopPad(JetSpace(Float64,2), -1:3)
m = [1.0, 2.0]
d = A*m # d = [0.0, 0.0, 1.0, 2.0, 0.0]
A = JopPad(JetSpace(Float64,2), -1:3, extend=true)
d = A*m # d = [1.0, 1.0, 1.0, 2.0, 2.0]
```

## 2D
```julia
A = JopPad(JetSpace(Float64,2,2), -1:3, 1:3)
m = [11. 12. ; 21. 22.]
d = A*m # d = [0. 0. 11. 12. 0. ; 0. 0. 21. 22. 0. ; 0. 0. 0. 0. 0.]
```

# Notes:
* This operator may also be used for truncation
* There may be overlap between the functionality of `JopPad` and `JopRestriction`, and it may be worth
thinking of how to consolidate them.
"""
function JopPad(dom::JetSpace{T}, pad::UnitRange...; extend=false, accumulate=false) where {T}
    JopLn(dom = dom, rng = JetSpace(T, map(length, pad)), df! = JopPad_df!, df′! = JopPad_df′!,
        s = (pad=pad, extend=extend, accumulate=accumulate))
end
export JopPad

function JopPad_df!(d::AbstractArray, m::AbstractArray; pad, extend, kwargs...)
    indices_rng = CartesianIndices(size(d))
    indices_dom = CartesianIndices(size(m))
    mpadₒ = modeloffset(pad)
    dpadₒ = dataoffset(pad)
    for didx in indices_rng
        midx = shiftidx(mpadₒ, didx)
        if midx ∈ indices_dom
            d[didx] = m[midx]
        elseif extend == false
            d[didx] = 0.0
        elseif extend == true
            d[didx] = m[modelidx_nn(mpadₒ, didx, size(m))]
        end
    end
    d
end

function JopPad_df′!(m::AbstractArray{T,N}, d::AbstractArray{T,N}; accumulate, pad, extend, kwargs...) where {T,N}
    if accumulate
        JopPad_df′!(Val(extend), m, d, (x,y)->x+y, pad)
    else
        m .= zero(T)
        JopPad_df′!(Val(extend), m, d, (x,y)->y, pad)
    end
    m
end

function JopPad_df′!(extend::Val{true}, m::AbstractArray{T,N}, d::AbstractArray{T,N}, f::Function, pad) where {T,N}
    indices_rng = CartesianIndices(size(d))
    indices_dom = CartesianIndices(size(m))
    mpadₒ = modeloffset(pad)
    dpadₒ = dataoffset(pad)
    for didx in indices_rng
        midx = shiftidx(mpadₒ, didx)
        if midx ∈ indices_dom
            m[midx] += d[didx]
        else
            m[modelidx_nn(mpadₒ, didx, size(m))] += d[didx]
        end
    end
    nothing
end

function JopPad_df′!(extend::Val{false}, m::AbstractArray{T,N}, d::AbstractArray{T,N}, f::Function, pad) where {T,N}
    dpadₒ = dataoffset(pad)
    for midx in intersection_domain(pad, size(m))
        didx = shiftidx(dpadₒ, midx)
        m[midx] = f(m[midx], d[didx])
    end
    nothing
end

function intersection_domain(pad::NTuple{N}, n::NTuple{N}) where {N}
    function _intersection_domain(i, pad, n)
        j = clamp(pad[i][1], 1, n[i])
        k = clamp(pad[i][end], 1, n[i])
        j:k
    end
    CartesianIndices(ntuple(i->_intersection_domain(i,pad,n), Val(N)))
end

modeloffset(pad::NTuple{N,UnitRange}) where {N} = ntuple(i->pad[i][1]-1, Val(N))
dataoffset(pad::NTuple{N,UnitRange}) where {N} = ntuple(i->1-pad[i][1], Val(N))
shiftidx(padₒ::NTuple{N,Int}, idx::CartesianIndex{N}) where {N} = CartesianIndex{N}(ntuple(dim->idx[dim]+padₒ[dim], N))
function modelidx_nn(padₒ::NTuple{N,Int}, idx::CartesianIndex{N}, n::NTuple{N,Int}) where N
    midx_try = shiftidx(padₒ, idx)
    CartesianIndex{N}(ntuple(dim->clamp(midx_try[dim], 1, n[dim]), Val(N)))
end


text:
"""
    SpectralVariable

Used to specify the type of the spectral variable which the refractive index is a function
of.  The available instances are `WAVELENGTH` for wavelength, `FREQUENCY` for frequency, and
`ENERGY` for photon energy.
"""
@enum SpectralVariable WAVELENGTH FREQUENCY ENERGY
for ins in instances(SpectralVariable); @eval export $(Symbol(ins)); end  # export all instances

"""
    read(path::String; <keyword arguments>)

Read the refractive index file at `path`, and return the wavelength vector `λ` and the
dielectric constant vector `ε` as a named tuple.

# Keywords
- `s_var::SpectralVariable=WAVELENGTH`: `WAVELENGTH`, `FREQUENCY`, or `ENERGY` if the refractive index file describes refractive index as a function of wavelength, frequency, or photon energy.
- `unit_prefix::Real=micro`: the prefix to multiply to the values of the spectral variable to convert their units to the SI units.  For example, `1e-6` for wavelengths in units of µm; `1.602176643` for photon energies in units of eV.
- `s_col::Integer=1`: the index of the spectral variable column in the refractive index file.
- `n_col::Integer=2`: the index of the refractive index column in the refractive index file.
- `kwargs...`: the keyword arguments to pass to `CSV.read()`
"""
function read(path::String;  # location of CSV file containing λ, n, and possibly k columns
              s_var::SpectralVariable=WAVELENGTH,
              unit_prefix::Real=micro,  # e.g., micro for μm; tera for THz; e⁺ for eV
              s_col::Integer=1,  # column index of spectral variable
              n_col::Integer=2,  # column index of refractive index
              kwargs...)  # keyword arguments for CSV.read; can be used to read .txt
    local df
    if splitpath(path)[1]=="https:"
        try
            print("First attempt to download from $path...  ")
            df = CSV.read(Downloads.download(path), DataFrame; kwargs...)
            println("Succeeds!")
        catch
            try
                println("Fails.")
                print("Second attempt to download from $path...  ")
                df = CSV.read(Downloads.download(path), DataFrame; kwargs...)
                println("Succeeds!")
            catch
                try
                    println("Fails.")
                    print("Final attempt to download from $path...  ")
                    df = CSV.read(Downloads.download(path), DataFrame; kwargs...)
                    println("Succeeds!")
                catch e
                    println("Fails.")
                    throw(e)
                end
            end
        end
    else
        df = CSV.read(path, DataFrame; kwargs...)
    end

    num_col = size(df, 2)
    num_col≥2 || @error "Data at $path contains $num_col columns, but should contain at least two columns corresponding to spectral variable λ and refractive index n."
    s = df[:,s_col] .* unit_prefix
    n = df[:,n_col]
    ε = n.^2

    if s_var==WAVELENGTH
        λ = s
    elseif s_var==FREQUENCY  # s = ν
        λ = c₀ ./ s  # λ = c₀ / ν
    else  # s_var==ENERGY; s = E
        λ = ℎc₀ ./ s  # λ = c₀ / ν = hc₀ / hν = ℎc₀ / E
    end

    # If sorted in descending order, reverse it.
    issorted(λ) || issorted(λ, rev=true) || @error "λ = $λ should be sorted in ascending or descending order."
    if issorted(λ, rev=true)  # sorted in descending order
        reverse!(λ)
        reverse!(ε)
    end

    return (; λ, ε)
end


text:
"""
    NoncentralT(ν, λ)
"""
struct NoncentralT{T<:Real} <: ContinuousUnivariateDistribution
    ν::T
    λ::T

    function NoncentralT{T}(ν::T, λ::T) where T
        @check_args(NoncentralT, ν > zero(ν))
        @check_args(NoncentralT, isreal(λ))
        new{T}(ν, λ)
    end
end

NoncentralT(ν::T, λ::T) where {T<:Real} = NoncentralT{T}(ν, λ)
NoncentralT(ν::Real, λ::Real) = NoncentralT(promote(ν, λ)...)
NoncentralT(ν::Integer, λ::Integer) = NoncentralT(Float64(ν), Float64(λ))

@distr_support NoncentralT -Inf Inf

### Conversions
convert(::Type{NoncentralT{T}}, ν::S, λ::S) where {T <: Real, S <: Real} = NoncentralT(T(ν), T(λ))
convert(::Type{NoncentralT{T}}, d::NoncentralT{S}) where {T <: Real, S <: Real} = NoncentralT(T(d.ν), T(d.λ))

### Parameters

params(d::NoncentralT) = (d.ν, d.λ)
@inline partype(d::NoncentralT{T}) where {T<:Real} = T

### Statistics

function mean(d::NoncentralT{T}) where T<:Real
    if d.ν > 1
        isinf(d.ν) ? d.λ :
        sqrt(d.ν/2) * d.λ * gamma((d.ν - 1)/2) / gamma(d.ν/2)
    else
        T(NaN)
    end
end

function var(d::NoncentralT{T}) where T<:Real
    d.ν > 2 ? d.ν*(1 + d.λ^2) / (d.ν - 2) - mean(d)^2 : T(NaN)
end

### Evaluation & Sampling

@_delegate_statsfuns NoncentralT ntdist ν λ

## sampling
function rand(rng::AbstractRNG, d::NoncentralT)
    ν = d.ν
    z = randn(rng)
    v = rand(rng, Chisq(ν))
    (z+d.λ)/sqrt(v/ν)
end