text:
# Autogenerated wrapper script for lapackpp_jll for i686-linux-musl-cxx03
export liblapackpp

using CompilerSupportLibraries_jll
using OpenBLAS32_jll
using blaspp_jll
JLLWrappers.@generate_wrapper_header("lapackpp")
JLLWrappers.@declare_library_product(liblapackpp, "liblapackpp.so")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, OpenBLAS32_jll, blaspp_jll)
    JLLWrappers.@init_library_product(
        liblapackpp,
        "lib/liblapackpp.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using CompScienceMeshes, BEAST
using LinearAlgebra
using Profile
using StaticArrays

function tau(x::SVector{U,T}) where {U,T}
    1.0-1.0/4.0
end

ntrc = X->ntrace(X,y)

T = tetmeshsphere(1.0,0.25)
X = nedelecd3d(T)
y = boundary(T)
@show numfunctions(X)

ϵ, μ, ω = 1.0, 1.0, 1.0; κ, η = ω * √(ϵ*μ), √(μ/ϵ)
ϵ_r =4.0
χ = tau
K, I, B = VIE.singlelayer(wavenumber=κ, tau=χ), Identity(), VIE.boundary(wavenumber=κ, tau=χ)
E = VIE.planewave(direction=ẑ, polarization=x̂, wavenumber=κ)
H = -1/(im*κ*η)*curl(E)

@hilbertspace j
@hilbertspace k
α = 1.0/ϵ_r
eq = @varform α*I[j,k]-K[j,k]-B[ntrc(j),k] == E[j]

dbvie = @discretise eq j∈X k∈X
u = solve(dbvie)

#postprocessing
Φ, Θ = [0.0], range(0,stop=2π,length=100)
pts = [point(cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)) for θ in Θ for ϕ in Φ]
ffd = potential(VIE.farfield(wavenumber=κ, tau=χ), pts, u, X)
ff = im*κ*ffd

using Plots

#Farfield
plot(xlabel="theta")
plot!(Θ, norm.(ff), label="far field", title="D-VIE")

#NearField
Z = range(-1,1,length=100)
Y = range(-1,1,length=100)
nfpoints = [point(0,y,z) for  y in Y, z in Z]

Enear = BEAST.grideval(nfpoints,α.*u,X)
Enear = reshape(Enear,100,100)

contour(real.(getindex.(Enear,1)))
heatmap(Z, Y,  real.(getindex.(Enear,1)), clim=(0.0,1.0))

ext: jl


text:
using Base64: Base64EncodePipe
using FileIO: save, @format_str, Stream

"""
    ImageContainer{format, S}

Stores raw image data as `format`.

supported format = [:jlc, :png, :svg, :jpg, :jpeg, :bmp, :gif, :mp4, :webm]
`:jlc` represents `Matrix{T<:Color}`.
"""
struct ImageContainer{format, S}
    content::S
end

function fileextension(file::AbstractString)
    _, ex = Base.Filesystem.splitext(file)
    if length(ex) <= 1
        ""
    else
        ex[2:end]
    end
end

"""
    storeimage(format::Symbol, data)
    storeimage(file::AbstractString)

Stores `data` as `format`.
Returns `ImageContainer{format, typeof(data)}`.

# Examples
```julia
c1 = storeimage(:png, read("image.png"))
c2 = storeimage("image.png")
```
"""
storeimage(format::Symbol, data) = ImageContainer{format, typeof(data)}(data)

function storeimage(file::AbstractString)
    ex = fileextension(file)
    if ex == ""
        error("Format identification failed: $file")
    else
        storeimage(Symbol(ex), read(file))
    end
end

const plainmimes = Dict(
    :png    => "image/png",
    :svg    => "image/svg+xml",
    :jpg    => "image/jpeg",
    :jpeg   => "image/jpeg",
    :bmp    => "image/bmp",     # for Juno
    :gif    => "image/gif"      # for Juno
)

# Without converting
for (fmt, mime) in plainmimes
    @eval function Base.show(io::IO, ::@MIME_str($mime),
                             c::ImageContainer{$(QuoteNode(fmt))})
        write(io, c.content)
    end
end

# With Base64 encoding
for fmt in (:gif, :bmp)
    @eval function Base.show(io::IO, ::MIME"text/html",
                             c::ImageContainer{$(QuoteNode(fmt))})
        write(io, "")
    end
end

# Videos
for fmt in (:mp4, :webm)
    @eval function Base.show(io::IO, ::MIME"text/html",
                             c::ImageContainer{$(QuoteNode(fmt))})
        write(io, "<video controls autoplay loop src=\"data:video/",
                  $(QuoteNode(fmt)), ";base64,")
        ioenc = Base64EncodePipe(io)
        write(ioenc, c.content)
        close(ioenc)
        write(io, "\" />")
    end
end

# Julia image = Matrix{<:Color}
# Needs converting to Bitmap
function Base.show(io::IO, ::MIME"text/html", c::ImageContainer{:jlc})
    write(io, "")
end

# Juno videos
for fmt in (:mp4, :webm)
    @eval function Base.show(io::IO,
                            ::MIME"application/prs.juno.plotpane+html",
                            c::ImageContainer{$(QuoteNode(fmt))})
        show(io, MIME("text/html"), c)
    end
end

# Juno Julia image
function Base.show(io::IO, ::MIME"image/bmp", c::ImageContainer{:jlc})
    save(Stream(format"BMP", io), c.content)
end


text:
# proj4 to sphere

# Pkg.add("NPZ")
# Pkg.add("FileIO")

using GLMakie, Proj4, NPZ, FileIO
using Colors
GLMakie.activate!()
let
    # download the data from here:
    # https://github.com/lazarusA/BeautifulMakie/tree/main/_assets/data
    x = npzread("aggregationOSMx.npy")
    y = npzread("aggregationOSMy.npy")
    aggOSM = npzread("aggregationOSM.npy")'
    img = load("osm_kbc.png") # osm_bmy.png
    # getting back to lon and lat coordinates
    source = Projection("+proj=webmerc +datum=WGS84")
    dest = Projection("+proj=longlat +datum=WGS84")
    lonlat = transform(source, dest, [x y])

    function toCartesian(lon, lat; r = 1, cxyz = (0, 0, 0))
        x = cxyz[1] + r * cosd(lat) * cosd(lon)
        y = cxyz[2] + r * cosd(lat) * sind(lon)
        z = cxyz[3] + r * sind(lat)
        return x, y, z
    end
    function lonlat3D(lon, lat, data; cxyz = (0.0, 0.0, 0.0))
        xyzw = zeros(size(data)..., 3)
        for (i, lon) in enumerate(lon), (j, lat) in enumerate(lat)
            x, y, z = toCartesian(lon, lat; cxyz = cxyz)
            xyzw[i, j, 1] = x
            xyzw[i, j, 2] = y
            xyzw[i, j, 3] = z
        end
        xyzw[:, :, 1], xyzw[:, :, 2], xyzw[:, :, 3]
    end
    x, y, z = lonlat3D(lonlat[:, 1], lonlat[:, 2], aggOSM)

    fig = Figure(resolution = (1200, 800), backgroundcolor = "#141414")
    ax = LScene(fig[1, 1], show_axis = false)
    GLMakie.surface!(ax, x, y, z; color = rotr90(img), shading = false,
        transparency = true)
    zoom!(ax.scene, cameracontrols(ax.scene), 0.65)
    GLMakie.rotate!(ax.scene, 12.0)
    display(fig)
end

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
"""
Using the algorithm proposed in the paper,
"Improving Walker's Algorithm to Run in Linear Time"
by Christoph Buchheim, Michael Junger, Sebastian Leipert
(http://dirk.jivas.de/papers/buchheim02improving.pdf)

Arguments
tree    Adjacency List that represents the given tree

Returns
positions     co-ordinates of the layout
"""
module Buchheim

using GeometryBasics

struct Tree{A<:AbstractVector,P<:AbstractVector,F}
    nodes::A
    mod::F
    thread::Vector{Int}
    ancestor::Vector{Int}
    prelim::F
    shift::F
    change::F
    positions::P
    nodesize::F
end

function Tree(tree::AbstractVector, nodesize)
    len = length(tree)
    mod = zeros(len)
    thread = zeros(Int, len)
    prelim = zeros(len)
    shift = zeros(len)
    change = zeros(len)
    ancestor = collect(1:len)
    nodes = copy(tree)
    positions = zeros(Point{2,Float64}, len)
    t = Tree(nodes, mod, thread, ancestor, prelim, shift, change, positions, nodesize)
    return t
end

function layout(t::AbstractVector; nodesize=ones(length(t)))
    return layout!(t, nodesize)
end

function layout!(t::AbstractVector, nodesize)
    tree = Tree(t, nodesize)
    first_walk(1, tree)
    second_walk(1, -tree.prelim[1], 0.0, tree)
    return tree.positions
end

function parent(v, t::Tree)
    tree = t.nodes
    for i in 1:length(tree)
        y = findall(x -> (x == v), tree[i])
        if length(y) != 0
            return i
        end
    end
    return nothing
end

function first_walk(v, t::Tree)
    prelim = t.prelim
    mod = t.mod
    tree = t.nodes
    nodesize = t.nodesize
    p = parent(v, t)
    if p != nothing
        index = findall(x -> (x == v), tree[p])[1]
    else
        index = 1
    end
    if length(tree[v]) == 0
        if v != tree[p][1]
            prelim[v] = prelim[tree[p][index - 1]] + (nodesize[tree[p][index - 1]])
        else
            prelim[v] = 0
        end
    else
        defaultAncestor = tree[v][1]
        for w in tree[v]
            first_walk(w, t)
            defaultAncestor = apportion(w, defaultAncestor, t)
        end
        execute_shifts(v, t)
        midpoint = (prelim[tree[v][1]] + prelim[tree[v][end]]) / 2
        if index > 1
            w = tree[p][index - 1]
            prelim[v] = prelim[w] + (nodesize[w] + 1.0)
            mod[v] = prelim[v] - midpoint
        else
            prelim[v] = midpoint
        end
    end
end

function apportion(v::T, defaultAncestor::T, t::Tree) where {T}
    tree = t.nodes
    ancestor = t.ancestor
    prelim = t.prelim
    mod = t.mod
    thread = t.thread
    p = parent(v, t)
    nodesize = t.nodesize
    if p != nothing
        index = findall(x -> (x == v), tree[p])[1]
    else
        index = 1
    end
    if index > 1
        w = tree[p][index - 1]
        v_in_right = v_out_right = v
        v_in_left = w
        v_out_left = tree[parent(v_in_right, t)][1]
        s_in_right = mod[v_in_right]
        s_out_right = mod[v_out_right]
        s_in_left = mod[v_in_left]
        s_out_left = mod[v_out_left]
        while next_right(v_in_left, t) != 0 && next_left(v_in_right, t) != 0
            v_in_left = next_right(v_in_left, t)
            v_in_right = next_left(v_in_right, t)
            v_out_left = next_left(v_out_left, t)
            v_out_right = next_right(v_out_right, t)
            ancestor[v_out_right] = v
            shift = (prelim[v_in_left] + s_in_left) - (prelim[v_in_right] + s_in_right) +
                    (nodesize[v_in_left])
            if shift > 0
                move_subtree(find_ancestor(v_in_left, v, defaultAncestor, t), v, shift, t)
                s_in_right += shift
                s_out_right += shift
            end
            s_in_left += mod[v_in_left]
            s_in_right += mod[v_in_right]
            s_out_left += mod[v_out_left]
            s_out_right += mod[v_out_right]
        end
        if next_right(v_in_left, t) != 0 && next_right(v_out_right, t) == 0
            thread[v_out_right] = next_right(v_in_left, t)
            mod[v_out_right] += s_in_left - s_out_right
        else
            if next_left(v_in_right, t) != 0 && next_left(v_out_left, t) == 0
                thread[v_out_left] = next_left(v_in_right, t)
                mod[v_out_left] += s_in_right - s_out_left
                defaultAncestor = v
            end
        end
    end
    return defaultAncestor
end

function number(v, t::Tree)
    p = parent(v, t)
    index = findall(x -> (x == v), t.nodes[p])[1]
    return index
end

function move_subtree(w_left::T, w_right::T, shift::Float64, t::Tree) where {T}
    change = t.change
    prelim = t.prelim
    tree = t.nodes
    shifttree = t.shift
    mod = t.mod
    n_wl = number(w_left, t)
    n_wr = number(w_right, t)
    subtrees = n_wr - n_wl
    change[w_right] -= shift / subtrees
    shifttree[w_right] += shift
    change[w_left] += shift / subtrees
    prelim[w_right] += shift
    return mod[w_right] += shift
end

function second_walk(v, m::Float64, depth::Float64, t::Tree)
    prelim = t.prelim
    mod = t.mod
    positions = t.positions
    nodesize = t.nodesize
    positions[v] = Point(prelim[v] + m, -depth)
    if length(t.nodes[v]) != 0
        maxdist = maximum([nodesize[i] for i in t.nodes[v]])
    else
        maxdist = 0
    end
    for w in t.nodes[v]
        second_walk(w, m + mod[v], Float64(depth + 1 + maxdist), t)
    end
end

function find_ancestor(w::T, v::T, defaultAncestor::T, tree::Tree) where {T}
    ancestor = tree.ancestor
    if ancestor[w] in tree.nodes[parent(v, tree)]
        return ancestor[w]
    else
        return defaultAncestor
    end
end

function execute_shifts(v, t::Tree)
    tree = t.nodes
    shift = t.shift
    change = t.change
    prelim = t.prelim
    mod = t.mod
    shiftnode = 0
    changenode = 0
    for w in reverse(tree[v])
        prelim[w] += shiftnode
        mod[w] += shiftnode
        changenode += change[w]
        shiftnode += shift[w] + changenode
    end
end

function next_left(v, t::Tree)
    tree = t.nodes
    thread = t.thread
    if length(tree[v]) != 0
        return tree[v][1]
    else
        return thread[v]
    end
end

function next_right(v, t::Tree)
    tree = t.nodes
    thread = t.thread
    if length(tree[v]) != 0
        return tree[v][end]
    else
        return thread[v]
    end
end

end #end of module


text:
"""
    ConstantGainApproximation()

Represents an approximation of the gain by a constant (in the Euclidean sense) vector field, given by the covariance of the observation function and x under the particle distribution.
"""
struct ConstantGainApproximation <: GainEstimationMethod end

function solve!(eq::PoissonEquation, method::ConstantGainApproximation)
    n, N, m = size(eq.gain)
    htilde = Htilde(eq)
    @inbounds for j in 1:m
        for i in 1:n
            g = 0.
            for l in 1:N
                g += eq.positions[i, l] * htilde[j, l]
            end
            for k in 1:N
                eq.gain[i, k, j] = g/N
            end
        end
    end
    return eq.gain
end

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
module SymbolicTracing

using SymbolicUtils, IRTools
using SymbolicUtils: Symbolic, Term, Sym
using IRTools: @dynamo, argument!, IR, isexpr

export @trace, trace

Base.@pure function sym_substitute(::Type{TSym}) where {T, TSym <: Symbolic{T}}
    Base.unwrap_unionall(T).abstract ? T : TSym
end
sym_substitute(::Type{T}) where {T} = T

@dynamo function sneakyinvoke(f, ::Type{T}, args...) where T<:Tuple
    ir = IR(f, T.parameters...)
    if ir === nothing
        sigs = " ::" * join(string.(args), " ::")
        error("Could not trace into $f($sigs). The method may not exist.\nRemember that trace will not recurse for Symbolic{T} for concretely typed T.")
    end
    argument!(ir, at = 2)
    return ir
end

@dynamo function trace(args...)
    ir = IR(args...)
    ir == nothing && return
    for (x, st) in ir
	isexpr(st.expr, :call) || continue
	ir[x] = Expr(:call, _trace, st.expr.args...)
    end
    return ir
end

_trace(f::Core.IntrinsicFunction, args...) = f(args...)
_trace(f) = trace(f)

@generated function _trace(f::F, args...) where {F}
    args′ = sym_substitute.(args)
    ex = if any(args .<: Symbolic)
        if fieldcount(F) == 0
            fi = F.instance
            if hasmethod(fi, Tuple{args...})
                :(f(args...))
            elseif fi isa Type
                :($Term(f, $Any[args...]))
            elseif hasmethod(fi, Tuple{args′...})
	        :($trace($sneakyinvoke, f, Tuple{$(args′...)}, args...))
            else
                sigs = " ::" * join(string.(args), ", ::")
                error("Could not trace into $fi($sigs). The method may not exist.\nRemember that trace will not recurse for Symbolic{T} for concretely typed T.")
            end
        else
            :($Term(f, $Any[args...]))
	end
    else
	:(f(args...))
    end
    "boo!"
    ci = expr_to_codeinfo(SymbolicTracing, [Symbol("#self#"), :f, :args], [:F], (F,), ex)
    ci.edges = []

    if fieldcount(F) == 0
        let T = Tuple{args...}, ftype = Tuple{f, args...}, ftype′ = Tuple{f, args′...}
            method_insts = Core.Compiler.method_instances(f.instance, T, typemax(UInt))
            covering_method_insts = [mi for mi in method_insts if ftype <: mi.def.sig || ftype′ <: mi.def.sig]
            method_doesnot_exist = isempty(covering_method_insts)
            # Now we add the edges so if a method is defined this recompiles
            if method_doesnot_exist
                # No method so attach to method table
                mt = f.name.mt
                typ = Base.rewrap_unionall(Tuple{f, Base.unwrap_unionall(T).parameters...}, T)
                push!(ci.edges, Core.Compiler.vect(mt, typ)...)
            else  # method exists, attach edges to all instances
                push!(ci.edges, covering_method_insts...)
            end
        end
    end
    return ci
end

"""
    expr_to_codeinfo(m::Module, argnames, spnames, sp, e::Expr)

Take an expr (usually a generated function generator) and convert it into a CodeInfo object
(Julia's internal, linear representation of code).

`m` is the module that the CodeInfo should be generated from (used for name resolution)

`argnames` must be an iterable container of symbols describing the CodeInfo's input arguments.
NOTE: the first argument should be given as `Symbol("#self#")`. So if the function is `f(x) = x + 1`,
then `argnames = [Symbol("#self#"), :x]`

`spnames` should be an iterable container of the names of the static parameters to the CodeInfo body
(e.g.) in `f(x::T) where {T <: Int} = ...`, `T` is a static parameter, so `spnames` should be `[:T]`

`sp` should be an iterable container of the static parameters to the CodeInfo body themselves (as
opposed to their names) (e.g.) in `f(x::T) where {T <: Int} = ...`, `T` is a static parameter,
so `sp` should be `[T]`

    `e` is the actual expression to lower to CodeInfo. This must be 'pure' in the same sense as generated
    function bodies.
    """
function expr_to_codeinfo(m::Module, argnames, spnames, sp, e::Expr)
    lam = Expr(:lambda, argnames,
               Expr(Symbol("scope-block"),
                    Expr(:block,
                         Expr(:return,
                              Expr(:block,
                                   e,
                                   )))))
    ex = if spnames === nothing || isempty(spnames)
        lam
    else
        Expr(Symbol("with-static-parameters"), lam, spnames...)
    end

    # Get the code-info for the generatorbody in order to use it for generating a dummy
    # code info object.
    ci = ccall(:jl_expand_and_resolve, Any, (Any, Any, Core.SimpleVector), ex, m, Core.svec(sp...))
    @assert ci isa Core.CodeInfo "Failed to create a CodeInfo from the given expression. This might mean it contains a closure or comprehension?\n Offending expression: $e"
    ci
end

macro trace(fcall)
    :($trace(() -> $fcall)) |> esc
end

end # module

max_stars_repo_name: MasonProtter/SymbolicTracing

max_issues_repo_name: MasonProtter/SymbolicTracing

max_forks_repo_name: MasonProtter/SymbolicTracing


text:
function init_s_grid(Ns)
    if any( (Ns .% 2) .== 0 )
        error("Ns must be odd numbers")
    end
    NN = round.(Int, (Ns-1)/2)
    Npoints = prod(Ns)
    s = zeros(3,Npoints)
    ip = 0
    for n = -NN[3]:NN[3]
        for m = -NN[2]:NN[2]
            for l = -NN[1]:NN[1]
                ip = ip + 1
                s[1,ip] = l/(2*NN[1] + 1)
                s[2,ip] = m/(2*NN[2] + 1)
                s[3,ip] = n/(2*NN[3] + 1)
            end
        end
    end
    return s
end

# r = hh * s
function s_to_r(s::Array{Float64,2}, hh)
    Npoints = size(s)[2]
    r = zeros(3,Npoints)
    for ip = 1:Npoints
        for α = 1:3
            r[α,ip] = 0.0
            for β = 1:3
                r[α,ip] = r[α,ip] + hh[α,β]*s[β,ip]
            end
        end
    end
    return r
end

# Function to evaluate LF
function eval_μ_lmn(Ns, hh, s_grid, ip, s)
    #
    NN = round.(Int, (Ns-1)/2)
    metric = 1/sqrt(det(hh))
    #
    s_lmn = s_grid[:,ip]
    #
    f = 0.0 + im*0.0
    # loop structure should be the same as the one used in init_s_grid
    for kz = -NN[3]:NN[3]
        for ky = -NN[2]:NN[2]
            for kx = -NN[1]:NN[1]
                ks = dot( [kx,ky,kz], (s - s_lmn) )
                f = f + cos(2π*ks)
            end
        end
    end
    return metric/sqrt(prod(Ns))*f
end

include("../../PW/common/gen_lattice.jl")

function test_main()
    # Sampling points
    Ns = [21,21,21]

    LL = gen_lattice_hexagonal(10)
    #LL = 10.0*eye(3,3)  # cubic
    hh = LL'

    s_grid = init_s_grid(Ns)

    ip = 2

    #r = [5.0, 5.0, 5.0]
    #s = inv(hh)*r
    s = s_grid[:,1]
    println( eval_μ_lmn(Ns, hh, s_grid, ip, s) )

    s = s_grid[:,2]
    println( eval_μ_lmn(Ns, hh, s_grid, ip, s) )
end

test_main()


text:
abstract type MeshAdaptiveDirectSearch end

mutable struct LTMADS <: MeshAdaptiveDirectSearch
  n::Int
  D::Array{Int,2}
  tau::Float64
  wminus::Int
  wplus::Int
  delta_m::Float64
  delta_p::Float64
  directionGenerator

  LTMADS(n, D = [Matrix{Float64}(I, n,n) -Matrix{Float64}(I, n,n)]) = begin
    # Default values taken from Audet2006:
    new(D, 4, -1, +1, 1, 1, LTMADSDirectionGenerator())
  end
end

# Eq 2.1 from Audet's paper:
function update_mesh_size(mads::MeshAdaptiveDirectSearch, meshSize)
  if improvedPointFound
    wk = rand(0:mads.wplus)
  else
    wk = rand(mads.wminus:-1)
  end
  meshSize * (mads.tau ^ wk)
end

# From page 203 of Audet2006:
function update_delta_m(m::LTMADS, isMinimalFrameCenter, improvedMeshPointFound)
  s = 4 # Audet2006 uses 4 but maybe we should adapt it instead?
  if isMinimalFrameCenter
    m.delta_m = m.delta_m / s     # Increases resolution
  elseif improvedMeshPointFound && (m.delta_m <= (1/s))
    m.delta_m = s * m.delta_m     # Decrease resolution
  else
    # do nothing
  end
end

# From box on page 203 in Audet2006:
mutable struct LTMADSDirectionGenerator
  cache::Dict{Int, Array{Int, 1}}  # Cache of previous directions used

  LTMADSDirectionGenerator() = begin
    new(Dict{Int, Array{Int, 1}}())
  end
end

function generate_direction(dg::LTMADSDirectionGenerator, l, n)
  if !haskey(dg.cache, l)
    # First create the general values
    twol = 2^l
    min = -twol + 1
    max = twol - 1
    bl = rand(min:max, n)

    # Then change one of the directions
    ihat = rand(1:n)
    bl[ihat] = (randbool() ? -1 : 1) * twol

    dg.cache[l] = bl
  end

  dg.cache[l]
end


text:
test_path = dirname(@__FILE__)
for test_file in readdir(test_path)
    test_file == "runtests.jl" && continue
    @show test_file
    include(joinpath(test_path, test_file))
end


text:
@testset "Transformations" begin
    function test_reverse(A, seq)
        revseq = reverse(LongSequence{A}(seq))
        @test convert(String, revseq) == reverse(seq)
        revseq = reverse!(LongSequence{A}(seq))
        @test convert(String, revseq) == reverse(seq)
    end

    function test_dna_complement(A, seq)
        comp = complement(LongSequence{A}(seq))
        @test convert(String, comp) == dna_complement(seq)
    end

    function test_rna_complement(A, seq)
        comp = complement(LongSequence{A}(seq))
        @test convert(String, comp) == rna_complement(seq)
    end

    function test_dna_revcomp(A, seq)
        revcomp = reverse_complement(LongSequence{A}(seq))
        @test convert(String, revcomp) == reverse(dna_complement(seq))
    end

    function test_rna_revcomp(A, seq)
        revcomp = reverse_complement(LongSequence{A}(seq))
        @test convert(String, revcomp) == reverse(rna_complement(seq))
    end

    @testset "Reverse" begin
        for len in [0, 1, 10, 32, 1000, 10000, 100000], _ in 1:10
            test_reverse(DNAAlphabet{4}, random_dna(len))
            test_reverse(RNAAlphabet{4}, random_rna(len))
            test_reverse(AminoAcidAlphabet, random_aa(len))
            test_reverse(CharAlphabet, random_aa(len))

            probs = [0.25, 0.25, 0.25, 0.25, 0.00]
            test_reverse(DNAAlphabet{2}, random_dna(len, probs))
            test_reverse(RNAAlphabet{2}, random_rna(len, probs))
        end
    end

    @testset "Complement" begin
        for len in [0, 1, 10, 32, 1000, 10000, 100000], _ in 1:10
            test_dna_complement(DNAAlphabet{4}, random_dna(len))
            test_rna_complement(RNAAlphabet{4}, random_rna(len))

            probs = [0.25, 0.25, 0.25, 0.25, 0.00]
            test_dna_complement(DNAAlphabet{2}, random_dna(len, probs))
            test_rna_complement(RNAAlphabet{2}, random_rna(len, probs))
        end
        seq_string = join(rand("-ACGTSWKMYRBDHVN", 1000))
        seq = complement(LongSequence{DNAAlphabet{4}}(seq_string))
        @test String(seq) == dna_complement(seq_string)

        seq_string = join(rand("-ACGUSWKMYRBDHVN", 1000))
        seq = complement(LongSequence{RNAAlphabet{4}}(seq_string))
        @test String(seq) == rna_complement(seq_string)
    end

    @testset "Reverse complement" begin
        for len in [0, 1, 10, 32, 1000, 10000, 100000], _ in 1:10
            test_dna_revcomp(DNAAlphabet{4}, random_dna(len))
            test_rna_revcomp(RNAAlphabet{4}, random_rna(len))

            probs = [0.25, 0.25, 0.25, 0.25, 0.00]
            test_dna_revcomp(DNAAlphabet{2}, random_dna(len, probs))
            test_rna_revcomp(RNAAlphabet{2}, random_rna(len, probs))
        end
        seq_string = join(rand("-ACGTSWKMYRBDHVN", 1000))
        seq = reverse_complement(LongSequence{DNAAlphabet{4}}(seq_string))
        @test String(seq) == reverse(dna_complement(seq_string))

        seq_string = join(rand("-ACGUSWKMYRBDHVN", 1000))
        seq = reverse_complement(LongSequence{RNAAlphabet{4}}(seq_string))
        @test String(seq) == reverse(rna_complement(seq_string))

        seq_string = join(rand("-ACGTSWKMYRBDHVN", 1000))
        seq = reverse_complement!(LongSequence{DNAAlphabet{4}}(seq_string))
        @test String(seq) == reverse(dna_complement(seq_string))

        seq_string = join(rand("-ACGTSWKMYRBDHVN", 1000))
        seq = LongSequence{DNAAlphabet{4}}(seq_string)
        seq2 = seq[100:200]
        reverse_complement!(seq2)
        @test String(seq) == seq_string
        @test String(seq2) == reverse(dna_complement(seq_string[100:200]))

	slice = randdnaseq(10)[2:9]
	@test reverse_complement(reverse_complement(slice)) == slice
    end

    @testset "Map" begin
        seq = dna""
        @test map(identity, seq) == dna""
        seq = dna"AAA"
        @test map(x -> DNA_C, seq) == dna"CCC"
        seq = dna"ACGTNACGTN"
        @test map(x -> x == DNA_N ? DNA_A : x, seq) == dna"ACGTAACGTA"
        @test seq == dna"ACGTNACGTN"
        @test map!(x -> x == DNA_N ? DNA_A : x, seq) === seq
        @test seq == dna"ACGTAACGTA"

        seq = rna""
        @test map(identity, seq) == rna""
        seq = rna"AAA"
        @test map(x -> RNA_C, seq) == rna"CCC"
        seq = rna"ACGUNACGUN"
        @test map(x -> x == RNA_N ? RNA_A : x, seq) == rna"ACGUAACGUA"
        @test seq == rna"ACGUNACGUN"
        @test map!(x -> x == RNA_N ? RNA_A : x, seq) === seq
        @test seq == rna"ACGUAACGUA"

        seq = aa""
        @test map(identity, seq) == aa""
        seq = aa"MMM"
        @test map(x -> AA_P, seq) == aa"PPP"
        seq = aa"XRNDCQXE"
        @test map(x -> x == AA_X ? AA_A : x, seq) == aa"ARNDCQAE"
        @test seq == aa"XRNDCQXE"
        @test map!(x -> x == AA_X ? AA_A : x, seq) === seq
        @test seq == aa"ARNDCQAE"
    end

    @testset "Filter" begin
        seq = dna""
        @test filter(x -> true, seq) == dna""
        @test filter(x -> false, seq) == dna""
        seq = dna"AAA"
        @test filter(x -> x == DNA_A, seq) == dna"AAA"
        @test filter(x -> x == DNA_C, seq) == dna""
        seq = dna"ACGTNACGTN"
        @test filter(x -> x == DNA_N, seq) == dna"NN"
        @test filter(x -> x != DNA_N, seq) == dna"ACGTACGT"
        @test seq == dna"ACGTNACGTN"
        @test filter!(x -> x != DNA_N, seq) == seq
        @test seq == dna"ACGTACGT"

        for len in 1:50, _ in 1:10
            str = random_dna(len)
            seq = LongDNASeq(str)
            @test filter(x -> x == DNA_N, seq) ==
                LongDNASeq(filter(x -> x == 'N', str))
            @test filter(x -> x != DNA_N, seq) ==
                LongDNASeq(filter(x -> x != 'N', str))
        end

        seq = rna""
        @test filter(x -> true, seq) == rna""
        @test filter(x -> false, seq) == rna""
        seq = rna"AAA"
        @test filter(x -> x == RNA_A, seq) == rna"AAA"
        @test filter(x -> x == RNA_C, seq) == rna""
        seq = rna"ACGUNACGUN"
        @test filter(x -> x == RNA_N, seq) == rna"NN"
        @test filter(x -> x != RNA_N, seq) == rna"ACGUACGU"
        @test seq == rna"ACGUNACGUN"
        @test filter!(x -> x != RNA_N, seq) == seq
        @test seq == rna"ACGUACGU"

        seq = aa""
        @test filter(x -> true, seq) == aa""
        @test filter(x -> false, seq) == aa""
        seq = aa"PPP"
        @test filter(x -> x == AA_P, seq) == aa"PPP"
        @test filter(x -> x != AA_P, seq) == aa""
        seq = aa"ARNDCQXGHILKMFPXTWYVOUX"
        @test filter(x -> x == AA_X, seq) == aa"XXX"
        @test filter(x -> x != AA_X, seq) == aa"ARNDCQGHILKMFPTWYVOU"
        @test seq == aa"ARNDCQXGHILKMFPXTWYVOUX"
        @test filter!(x -> x != AA_X, seq) == seq
        @test seq == aa"ARNDCQGHILKMFPTWYVOU"
    end
end


text:
function mw_marg_lf(x::Vector{Float64}, cov::Array{Float64,2})
# This is a scaled version of the density in King and Hillier (1985) Equation (2) @
# -- Note x is ncov x n ... compute n marginal likelihood values @
# Compute Marginal Likelihood value for x with covariance matrix cov
  x = x./sqrt(sum(x.^2, 1))'
  m = size(x, 1)
  d = det(cov)
  # (gauss code) x = chol(invpd(cov))
  c = chol(Hermitian(inv(x'x)))
  print(x)
  #x = c*x
  tmp = sum(x, 1)
  tmp = d*(tmp.^m)
  mlf = sqrt(1.0./tmp)[:]
  return mlf
end
function mw_marg_lf(x::Array{Float64,2}, cov::Array{Float64,2})
  nsim = size(x, 2)
  println(nsim)
  mlf_vec = Vector{Float64}(nsim)
  for isim ∈ 1:nsim
    mlf_vec[isim] = mw_marg_lf(x[:, isim], cov)[1]
  end
return mlf_vec
end


text:
# This file was generated, do not modify it. # hide
DTR = @load DecisionTreeRegressor pkg=DecisionTree

y, X = unpack(df, ==(:price))
train, test = partition(collect(eachindex(y)), 0.7, shuffle=true, rng=5)
tree = machine(DTR(), X, y)

fit!(tree, rows=train);

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using SlabWaveguides
using Test

@testset "SlabWaveguides.jl" begin
    # Write your tests here.
end


text:
# Julia wrapper for header: /usr/local/ffmpeg/2.4/include/libavutil/imgutils.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_image_fill_max_pixsteps,
    av_image_get_linesize,
    av_image_fill_linesizes,
    av_image_fill_pointers,
    av_image_alloc,
    av_image_copy_plane,
    av_image_copy,
    av_image_fill_arrays,
    av_image_get_buffer_size,
    av_image_copy_to_buffer,
    av_image_check_size,
    av_image_check_sar

function av_image_fill_max_pixsteps(max_pixsteps,max_pixstep_comps,pixdesc)
    ccall((:av_image_fill_max_pixsteps,libavutil),Void,(Ptr{Cint},Ptr{Cint},Ptr{AVPixFmtDescriptor}),max_pixsteps,max_pixstep_comps,pixdesc)
end

function av_image_get_linesize(pix_fmt::AVPixelFormat,width::Integer,plane::Integer)
    ccall((:av_image_get_linesize,libavutil),Cint,(AVPixelFormat,Cint,Cint),pix_fmt,width,plane)
end

function av_image_fill_linesizes(linesizes,pix_fmt::AVPixelFormat,width::Integer)
    ccall((:av_image_fill_linesizes,libavutil),Cint,(Ptr{Cint},AVPixelFormat,Cint),linesizes,pix_fmt,width)
end

function av_image_fill_pointers(data,pix_fmt::AVPixelFormat,height::Integer,ptr,linesizes)
    ccall((:av_image_fill_pointers,libavutil),Cint,(Ptr{Ptr{Uint8}},AVPixelFormat,Cint,Ptr{Uint8},Ptr{Cint}),data,pix_fmt,height,ptr,linesizes)
end

function av_image_alloc(pointers,linesizes,w::Integer,h::Integer,pix_fmt::AVPixelFormat,align::Integer)
    ccall((:av_image_alloc,libavutil),Cint,(Ptr{Ptr{Uint8}},Ptr{Cint},Cint,Cint,AVPixelFormat,Cint),pointers,linesizes,w,h,pix_fmt,align)
end

function av_image_copy_plane(dst,dst_linesize::Integer,src,src_linesize::Integer,bytewidth::Integer,height::Integer)
    ccall((:av_image_copy_plane,libavutil),Void,(Ptr{Uint8},Cint,Ptr{Uint8},Cint,Cint,Cint),dst,dst_linesize,src,src_linesize,bytewidth,height)
end

function av_image_copy(dst_data,dst_linesizes,src_data,src_linesizes,pix_fmt::AVPixelFormat,width::Integer,height::Integer)
    ccall((:av_image_copy,libavutil),Void,(Ptr{Ptr{Uint8}},Ptr{Cint},Ptr{Ptr{Uint8}},Ptr{Cint},AVPixelFormat,Cint,Cint),dst_data,dst_linesizes,src_data,src_linesizes,pix_fmt,width,height)
end

function av_image_fill_arrays(dst_data,dst_linesize,src,pix_fmt::AVPixelFormat,width::Integer,height::Integer,align::Integer)
    ccall((:av_image_fill_arrays,libavutil),Cint,(Ptr{Ptr{Uint8}},Ptr{Cint},Ptr{Uint8},AVPixelFormat,Cint,Cint,Cint),dst_data,dst_linesize,src,pix_fmt,width,height,align)
end

function av_image_get_buffer_size(pix_fmt::AVPixelFormat,width::Integer,height::Integer,align::Integer)
    ccall((:av_image_get_buffer_size,libavutil),Cint,(AVPixelFormat,Cint,Cint,Cint),pix_fmt,width,height,align)
end

function av_image_copy_to_buffer(dst,dst_size::Integer,src_data,src_linesize,pix_fmt::AVPixelFormat,width::Integer,height::Integer,align::Integer)
    ccall((:av_image_copy_to_buffer,libavutil),Cint,(Ptr{Uint8},Cint,Ptr{Ptr{Uint8}},Ptr{Cint},AVPixelFormat,Cint,Cint,Cint),dst,dst_size,src_data,src_linesize,pix_fmt,width,height,align)
end

function av_image_check_size(w::Integer,h::Integer,log_offset::Integer,log_ctx)
    ccall((:av_image_check_size,libavutil),Cint,(Uint32,Uint32,Cint,Ptr{Void}),w,h,log_offset,log_ctx)
end

function av_image_check_sar(w::Integer,h::Integer,sar::AVRational)
    ccall((:av_image_check_sar,libavutil),Cint,(Uint32,Uint32,AVRational),w,h,sar)
end


text:
import Flux: OneHotVector
import CuArrays: cudaconvert
import Base.Broadcast: BroadcastStyle, ArrayStyle

import .Basic: indices2onehot, onehot2indices, tofloat, _labelindex

Base.similar(o::OneHotArray{N, <:CuArray}) where N = similar(o, size(o))
Base.similar(o::OneHotArray{N, <:CuArray}, dims::NTuple{N, Int}) where N = similar(o, Bool, dims)
Base.similar(o::OneHotArray{N, <:CuArray}, dims::Int...) where N = similar(o, Bool, Tuple(dims))
Base.similar(o::OneHotArray{N, <:CuArray}, T::Type) where N = similar(o, T, size(o))
Base.similar(o::OneHotArray{N, <:CuArray}, T::Type, dims::Int...) where N = similar(o, T, Tuple(dims))
Base.similar(o::OneHotArray{N, <:CuArray}, T::Type, dims::NTuple{N, Int}) where N = CuArray{T}(undef, dims)

tofloat(::Type{F}, o::OneHotArray{N, <:CuArray}) where {F<:AbstractFloat,N} = CuArray{F}(o)

BroadcastStyle(::Type{<:OneHotArray{N, <:CuArray} where N}) = ArrayStyle{CuArray}()
cudaconvert(x::OneHotArray{<:CuArray}) = OneHotArray(x.dims, cudaconvert(x.data))

#gpu indices to onehot
indices2onehot(num::Int, xs::CuArray{Int}) = convert(CuArray{OneHotVector}, _labelindex(num, xs))

#gpu onehot to indices
onehot2indices(x::CuArray{OneHotVector}) = convert(CuArray{Int}, x)


text:
Base.show(io::IO,h::HilbertSpace) = write(io, "ℋ(", h.name, ")")
function Base.show(io::IO,h::ProductSpace)
    show(io, h.spaces[1])
    for i=2:length(h.spaces)
        write(io, " ⊗ ")
        show(io, h.spaces[i])
    end
end
function Base.show(io::IO,h::ClusterSpace)
    write(io, "$(h.N)x")
    show(io, h.original_space)
end

Base.show(io::IO,x::QSym) = write(io, x.name)
Base.show(io::IO,x::Create) = write(io, string(x.name, "′"))
Base.show(io::IO,x::Transition) = write(io, Symbol(x.name,x.i,x.j))

show_brackets = Ref(true)
function Base.show(io::IO,x::QTerm)
    show_brackets[] && write(io,"(")
    show(io, SymbolicUtils.arguments(x)[1])
    f = SymbolicUtils.operation(x)
    for i=2:length(SymbolicUtils.arguments(x))
        show(io, f)
        show(io, SymbolicUtils.arguments(x)[i])
    end
    show_brackets[] && write(io,")")
end

function Base.show(io::IO,x::QMul)
    if !SymbolicUtils._isone(x.arg_c)
        show(io, x.arg_c)
        show(io, *)
    end
    show_brackets[] && write(io, "(")
    show(io, x.args_nc[1])
    for i=2:length(x.args_nc)
        show(io, *)
        show(io, x.args_nc[i])
    end
    show_brackets[] && write(io, ")")
end

function SymbolicUtils.show_term(io::IO, t::SymbolicUtils.Term{<:AvgSym})
    write(io, "⟨")
    show_brackets[] = false
    show(io, SymbolicUtils.arguments(t)[1])
    show_brackets[] = true
    write(io, "⟩")
end

function Base.show(io::IO,de::AbstractMeanfieldEquations)
    for i=1:length(de.equations)
        write(io, "∂ₜ(")
        show(io, de.equations[i].lhs)
        write(io, ") = ")
        show(io, de.equations[i].rhs)
        write(io, "\n")
    end
end

function Base.show(io::IO, c::CorrelationFunction)
    show(io, average(c.op1*c.op2))
end
function Base.show(io::IO, S::Spectrum)
    write(io, "ℱ(")
    show(io, S.corr)
    write(io, ")(ω)")
end

Base.show(io::IO, c::CallableTransition) = write(io, c.name)

const T_LATEX = Union{<:QNumber,<:AbstractMeanfieldEquations, <:SymbolicUtils.Symbolic{<:CNumber},
        <:CorrelationFunction,<:Spectrum}
Base.show(io::IO, ::MIME"text/latex", x::T_LATEX) = write(io, latexify(x))


text:
#####################
# general parser infrastructure for circuits
#####################

# The `ParserCombinator` library works correctly but is orders of magnitude too slow.
# Instead here we hardcode some simpler parsers to speed things up

"""
Load a logical circuit from file.
Support file formats:
 * ".sdd" for SDD files
 * ".psdd" for PSDD files
 * ".circuit" for Logistic Circuit files
"""
function load_logical_circuit(file::String)::UnstLogicalΔ
    compile_logical(parse_circuit_file(file))
end

"""
Load a smooth logical circuit from file.
Support file formats:
 * ".psdd" for PSDD files
 * ".circuit" for Logistic Circuit files
"""
function load_smooth_logical_circuit(file::String)::UnstLogicalΔ
    compile_smooth_logical(parse_circuit_file(file))
end

"""
Load a smooth structured logical circuit from file.
Support circuit file formats:
 * ".psdd" for PSDD files
 * ".circuit" for Logistic Circuit files
Supported vtree file formats:
 * ".vtree" for VTree files
"""
function load_struct_smooth_logical_circuit(circuit_file::String, vtree_file::String)::Tuple{StructLogicalΔ{PlainVtreeNode},PlainVtree}
    circuit_lines = parse_circuit_file(circuit_file)
    vtree_lines = parse_vtree_file(vtree_file)
    compile_smooth_struct_logical(circuit_lines, vtree_lines)
end

#####################
# parse based on file extension
#####################

function parse_circuit_file(file::String)::CircuitFormatLines
    if endswith(file,".circuit")
        parse_lc_file(file)
    elseif endswith(file,".psdd")
        parse_psdd_file(file)
    elseif endswith(file,".sdd")
        parse_sdd_file(file)
    else
        throw("Cannot parse this file type as a circuit: $file")
    end
end

#####################
# parser of logistic circuit file format
#####################

const parens = r"\(([^\)]+)\)"

function parse_lc_decision_line(ln::String)::DecisionLine{LCElement}
    @assert startswith(ln, "D")
    head::SubString, tail::SubString = split(ln,'(',limit=2)
    head_tokens = split(head)
    head_ints::Vector{UInt32} = map(x->parse(UInt32,x),head_tokens[2:4])
    elems_str::String = "("*tail
    elems = Vector{LCElement}()
    for x in eachmatch(parens::Regex, elems_str)
        tokens = split(x[1], limit=3)
        weights::Vector{Float64} = map(x->parse(Float64,x), split(tokens[3]))
        elem = LCElement(parse(UInt32,tokens[1]), parse(UInt32,tokens[2]), weights)
        push!(elems,elem)
    end
    DecisionLine(head_ints[1],head_ints[2],head_ints[3],elems)
end

function parse_lc_literal_line(ln::String)::WeightedLiteralLine
    @assert startswith(ln, "T") || startswith(ln, "F")
    tokens = split(ln)
    head_ints = map(x->parse(UInt32,x),tokens[2:4])
    weights = map(x->parse(Float64,x), tokens[5:end])
    lit = var2lit(head_ints[3])
    if startswith(ln, "F")
       lit = -lit # negative literal
    end
    WeightedLiteralLine(head_ints[1],head_ints[2],lit,true,weights)
end

function parse_comment_line(ln::String)
    @assert startswith(ln, "c")
    CircuitCommentLine(lstrip(chop(ln, head = 1, tail = 0)))
end

function parse_lc_header_line(ln::String)
    @assert (ln == "Logistic Circuit") || (ln == "Logisitic Circuit")
    LcHeaderLine()
end

function parse_bias_line(ln::String)::BiasLine
    @assert startswith(ln, "B")
    tokens = split(ln)
    weights = map(x->parse(Float64,x), tokens[2:end])
    BiasLine(weights)
end

parse_lc_file(file::String)::CircuitFormatLines = open(parse_lc_file, file)

function parse_lc_file(file::Core.IO)::CircuitFormatLines
    q = Vector{CircuitFormatLine}()
    for ln in eachline(file)
        @assert !isempty(ln)
        if ln[1] == 'D'
            push!(q, parse_lc_decision_line(ln))
        elseif ln[1] == 'T' || ln[1] == 'F'
            push!(q, parse_lc_literal_line(ln))
        elseif ln[1] == 'c'
            push!(q, parse_comment_line(ln))
        elseif ln[1] == 'L'
            push!(q, parse_lc_header_line(ln))
        elseif ln[1] == 'B'
            push!(q, parse_bias_line(ln))
        else
            error("Cannot parse logistic circuit file format line '$ln'")
        end
    end
    q
end

#####################
# parser for PSDD circuit file format
#####################

function parse_psdd_decision_line(ln::String)::DecisionLine{PSDDElement}
    @assert startswith(ln, "D")
    tokens = split(ln)
    head_ints::Vector{UInt32} = map(x->parse(UInt32,x),tokens[2:4])
    elems = Vector{PSDDElement}()
    for (p,s,w) in Iterators.partition(tokens[5:end],3)
        prime = parse(UInt32,p)
        sub = parse(UInt32,s)
        weight = parse(Float64,w)
        elem = PSDDElement(prime, sub, weight)
        push!(elems,elem)
    end
    DecisionLine(head_ints[1],head_ints[2],head_ints[3],elems)
end

function parse_psdd_true_leaf_line(ln::String)::WeightedNamedConstantLine
    @assert startswith(ln, "T")
    tokens = split(ln)
    @assert length(tokens)==5
    head_ints = map(x->parse(UInt32,x),tokens[2:4])
    weight = parse(Float64,tokens[5])
    WeightedNamedConstantLine(head_ints[1],head_ints[2],head_ints[3],weight)
end

function parse_literal_line(ln::String, normalized::Bool)::UnweightedLiteralLine
    @assert startswith(ln, "L")
    tokens = split(ln)
    @assert length(tokens)==4 "line has too many tokens: $ln"
    head_ints = map(x->parse(UInt32,x),tokens[2:3])
    lit = parse(Int32,tokens[4])
    UnweightedLiteralLine(head_ints[1],head_ints[2],lit,normalized)
end

function parse_psdd_header_line(ln::String)
    @assert startswith(ln, "psdd")
    PsddHeaderLine(parse(Int,split(ln)[2]))
end

parse_psdd_file(file::String)::CircuitFormatLines = open(parse_psdd_file, file)

function parse_psdd_file(file::Core.IO)::CircuitFormatLines
    q = Vector{CircuitFormatLine}()
    for ln in eachline(file)
        @assert !isempty(ln)
        if ln[1] == 'D'
            push!(q, parse_psdd_decision_line(ln))
        elseif ln[1] == 'T'
            push!(q, parse_psdd_true_leaf_line(ln))
        elseif ln[1] == 'L'
            push!(q, parse_literal_line(ln, true))
        elseif ln[1] == 'c'
            push!(q, parse_comment_line(ln))
        elseif startswith(ln,"psdd")
            push!(q, parse_psdd_header_line(ln))
        else
            error("Cannot parse PSDD file format line '$ln'")
        end
    end
    q
end

#####################
# parser for SDD circuit file format
#####################

function parse_sdd_decision_line(ln::String)::DecisionLine{SDDElement}
    @assert startswith(ln, "D")
    tokens = split(ln)
    head_ints::Vector{UInt32} = map(x->parse(UInt32,x),tokens[2:4])
    elems = Vector{SDDElement}()
    for (p,s) in Iterators.partition(tokens[5:end],2)
        prime = parse(UInt32,p)
        sub = parse(UInt32,s)
        elem = SDDElement(prime, sub)
        push!(elems,elem)
    end
    DecisionLine(head_ints[1],head_ints[2],head_ints[3],elems)
end

function parse_sdd_constant_leaf_line(ln::String)::AnonymousConstantLine
    @assert startswith(ln, "T") || startswith(ln, "F")
    tokens = split(ln)
    @assert length(tokens)==2
    AnonymousConstantLine(parse(UInt32,tokens[2]), startswith(ln, "T"), false)
end

function parse_sdd_header_line(ln::String)
    @assert startswith(ln, "sdd")
    SddHeaderLine(parse(Int,split(ln)[2]))
end

parse_sdd_file(file::String)::CircuitFormatLines = open(parse_sdd_file, file)

function parse_sdd_file(file::Core.IO)::CircuitFormatLines
    q = Vector{CircuitFormatLine}()
    for ln in eachline(file)
        @assert !isempty(ln)
        if ln[1] == 'D'
            push!(q, parse_sdd_decision_line(ln))
        elseif ln[1] == 'T' || ln[1] == 'F'
            push!(q, parse_sdd_constant_leaf_line(ln))
        elseif ln[1] == 'L'
            push!(q, parse_literal_line(ln, false))
        elseif ln[1] == 'c'
            push!(q, parse_comment_line(ln))
        elseif startswith(ln,"sdd")
            push!(q, parse_sdd_header_line(ln))
        else
            error("Cannot parse SDD file format line '$ln'")
        end
    end
    q
end

#####################
# loader for CNF file format
#####################

"""
Load a CNF as a logical circuit from file.
Supppor file formats:
* ".cnf" for CNF files
"""
function load_cnf(file::String)::UnstLogicalΔ
    @assert endswith(file, ".cnf")

    # linearized circuit nodes
    circuit = Vector{UnstLogicalΔNode}()

    # linearized clauses (disjunctions)
    clauses = Vector{⋁Node}()

    # literal cache is responsible for making leaf literals nodes unique and adding them to `circuit`
    lit_cache = Dict{Lit,LogicalLeafNode}()
    literal_node(l::Lit) = get!(lit_cache, l) do
        leaf = LiteralNode(l)
        push!(circuit,leaf) # also add new leaf to linearized circuit before caller
        leaf
    end

    # record the current clause
    clause = ⋁Node([])

    count_clauses = 0

    open(file) do file

        for ln in eachline(file)
            @assert !isempty(ln)
            if ln[1] == 'c' || startswith(ln, "p cnf")
                # skip comment and header lines
                continue
            else
                tokens = split(ln)
                for token in tokens
                    if !occursin(r"^\s*[-]?[0-9]+\s*$", token)
                        error("Cannot parse CNF file format line '$ln'")
                    end
                    literal = parse(Lit, token)
                    if literal == 0
                        push!(clauses, clause)
                        push!(circuit, clause)
                        clause = ⋁Node([])
                        count_clauses += 1
                    else
                        push!(clause.children, literal_node(literal))
                    end
                end
            end
        end

        # handle the last clause
        if length(clause.children) > 0
            push!(clauses, clause)
            push!(circuit, clause)
            count_clauses + 1
        end

        # create the root conjunction node
        if length(clauses) > 0
            push!(circuit, ⋀Node(clauses))
        end
    end

    circuit
end

#####################
# loader for DNF file format
#####################

"""
Load a CNF as a logical circuit from file.
Supppor file formats:
* ".cnf" for CNF files
"""
function load_dnf(file::String)::UnstLogicalΔ
    @assert endswith(file, ".dnf")

    # linearized circuit nodes
    circuit = Vector{UnstLogicalΔNode}()

    # linearized clauses (conjunctions)
    clauses = Vector{⋀Node}()

    # literal cache is responsible for making leaf literals nodes unique and adding them to `circuit`
    lit_cache = Dict{Lit,LogicalLeafNode}()
    literal_node(l::Lit) = get!(lit_cache, l) do
        leaf = LiteralNode(l)
        push!(circuit,leaf) # also add new leaf to linearized circuit before caller
        leaf
    end

    # record the current clause
    clause = ⋀Node([])

    open(file) do file

        for ln in eachline(file)
            @assert !isempty(ln)
            if ln[1] == 'c' || startswith(ln, "p dnf")
                # skip comment and header lines
                continue
            else
                tokens = split(ln)
                for token in tokens
                    if !occursin(r"^\s*[-]?[0-9]+\s*$", token)
                        error("Cannot parse CNF file format line '$ln'")
                    end
                    literal = parse(Lit, token)
                    if literal == 0
                        push!(clauses, clause)
                        push!(circuit, clause)
                        clause = ⋀Node([])
                    else
                        push!(clause.children, literal_node(literal))
                    end
                end
            end
        end

        # handle the last clause
        if length(clause.children) > 0
            push!(clauses, clause)
            push!(circuit, clause)
        end

        # create the root disjunction node
        if length(clauses) > 0
            push!(circuit, ⋁Node(clauses))
        end
    end

    circuit
end


text:
###############################################################################
#
# Vary the traits in the SPACSimple struct
#
###############################################################################
"""
    vary_spac!(node::SPACSimple{FT},
               weat_years::DataFrame,
               factor_to_vary::String,
               ratio_to_vary::FT
    ) where {FT<:AbstractFloat}

Generalized function to vary the SoilPlantAirContinuum to run sensitivity
    analysis, given
- `node` SPACSimple type struct
- `weat_years` Weather data
- `factor_to_vary` Which parameter to vary
- `ratio_to_vary` How much to vary
"""
function vary_spac!(
            node::SPACSimple{FT},
            weat_years::DataFrame,
            factor_to_vary::String,
            ratio_to_vary::FT
) where {FT<:AbstractFloat}
    # choose which parameter(s) to vary
    if factor_to_vary=="wb"
        (node.hs.root.vc).b *= ratio_to_vary;
        (node.hs.stem.vc).b *= ratio_to_vary;
        (node.hs.leaf.vc).b *= ratio_to_vary;
    elseif factor_to_vary=="wc"
        (node.hs.root.vc).c *= ratio_to_vary;
        (node.hs.stem.vc).c *= ratio_to_vary;
        (node.hs.leaf.vc).c *= ratio_to_vary;
    elseif factor_to_vary=="wk"
        node.hs.root.k_max *= ratio_to_vary;
        node.hs.stem.k_max *= ratio_to_vary;
        node.hs.leaf.k_sla *= ratio_to_vary;
        node.hs.root.k_element .*= ratio_to_vary;
        node.hs.stem.k_element .*= ratio_to_vary;
        node.hs.leaf.k_element .*= ratio_to_vary;
    elseif factor_to_vary=="kw"
        node.hs.root.k_max *= ratio_to_vary;
        node.hs.stem.k_max *= ratio_to_vary;
        node.hs.root.k_element .*= ratio_to_vary;
        node.hs.stem.k_element .*= ratio_to_vary;
    elseif factor_to_vary=="kl"
        node.hs.leaf.k_sla *= ratio_to_vary;
        node.hs.leaf.k_element .*= ratio_to_vary;
    elseif factor_to_vary=="cv"
        node.c_vmax *= ratio_to_vary;
    elseif factor_to_vary=="cc"
        node.c_cons *= ratio_to_vary;
    elseif factor_to_vary=="gm"
        node.g_max *= ratio_to_vary;
    elseif factor_to_vary=="ga"
        node.gaba *= ratio_to_vary;
        node.lai  /= ratio_to_vary;
    elseif factor_to_vary=="sd"
        node.h_soil *= ratio_to_vary;
        node.hs.root.p_gravity .*= ratio_to_vary;
    elseif factor_to_vary=="ca"
        weat_years[:,"CO2"] *= ratio_to_vary;
    elseif factor_to_vary=="rh"
        weat_years[:, "D"] *= ratio_to_vary;
    elseif factor_to_vary=="ta"
        _multip1 = saturation_vapor_pressure.((weat_years).Tair .+ T_0(FT));
        weat_years[:, "Tair"] .+= ratio_to_vary;
        _multip2 = saturation_vapor_pressure.((weat_years).Tair .+ T_0(FT));
        weat_years[:, "D"] .*= _multip2 ./ _multip1;
    else
        @warn "Invalid parameter provided, nothing has been changed";
    end

    return nothing
end


text:
type linear_solver_JULIA <: abstract_linear_system_solver
	_SparseMatrix::SparseMatrixCSC{Float64,Int64}
	_factor # TO DO, give type

	# options
	sym::Symbol

  function linear_solver_JULIA()
      return new();
  end
end

function ls_factor!(solver::linear_solver_JULIA, n::Int64, m::Int64)
			if solver.sym == :unsymmetric
				 solver._factor = lufact(solver._SparseMatrix);
			elseif solver.sym == :definite
				solver._factor = cholfact(solver._SparseMatrix);
			else
				error("this.options.sym = " * string(solver.sym) * " not supported")
			end

			return 1; # inertia
end

function ls_solve!(solver::linear_solver_JULIA, my_rhs::Array{Float64,1}, my_sol::Array{Float64,1})
  #GLOBAL_timer::class_algorithm_timer
  #GLOBAL_timer.start("Solve/julia")
  my_sol[1:length(my_sol)] = solver._factor \ my_rhs; #::UmfpackLU{Float64,Int64}
  #GLOBAL_timer.stop("Solve/julia")
end

function ls_solve(solver::linear_solver_JULIA, my_rhs::AbstractArray)
  #GLOBAL_timer::class_algorithm_timer
  #GLOBAL_timer.start("Solve/julia")
  return solver._factor \ my_rhs; #::UmfpackLU{Float64,Int64}
  #GLOBAL_timer.stop("Solve/julia")
end

max_forks_repo_name: ohinder/Everest


text:
using ..SymEngine
using ..SymEngine: @vars, Basic
export @vars, Basic, subs

include("register.jl")
include("instruct.jl")
include("blocks.jl")
include("patch.jl")


text:
##############
# Bands
##############
#
# TODO: add extrapolation too (second half of make_bands)

export map_onto_bands, make_1Dglacier, map_back_to_2D, map_back_to_2D!

# used to round (up or down) the binsize to the next decimal place
_binround(binsize::Number) = -floor(Int, log10(abs(binsize)))
_binround(binsize) = 0

"""
    make_1Dglacier(dem::Gridded, binsize_or_bins, glaciermask=BitArray([]);
                        binround=_binround(binsize_or_bins),
                        window_dem_smooth=0.0,
                        window_width_smooth=0.0,
                        alpha_min=deg2rad(0.4),
                        alpha_max=deg2rad(60.0),
                        FILL=-9999999.0)

Makes a 1D glacier from a 2D DEM by using Huss' elevation band trick.

Returns:
- bands -- elevation bands.  The i-th band is (bands[i], bands[i]+step(bands))
- bandi -- linear index into dem which assigns each cell to a elevation band
- alphas, areas, lengths, widths, x, xmid -- elevation band slope, area, etc
- dem -- the used DEM, a possibly smoothed version of the input DEM
- alpha2d -- slopes at each point of `dem`
"""
function make_1Dglacier(dem::Gridded, binsize_or_bins, glaciermask=trues(size(dem.v));
                        binround=_binround(binsize_or_bins),
                        min_bin_number_ends=0,
                        min_bands=4,
                        window_dem_smooth=0.0,
                        window_width_smooth=0.0,
                        alpha_min=deg2rad(0.4),
                        alpha_max=deg2rad(60.0),
                        FILL=-9999999.0,
                        verbose=true)
    dx = step(dem.x)
    # Smooth dem to get smooth alpha, smoother bands.  This is in
    # particular important when there is cross-flow bumpiness, such as
    # on Uaar.  However, it can also be bad.  YMMV, check!
    if window_dem_smooth>0
        dem = deepcopy(dem)
        nofillmask = dem.v.!=FILL
        mask = nofillmask .& glaciermask
        dem.v[:] = boxcar(dem.v, round(Int,window_dem_smooth/dx), mask, (!).(mask))
        dem.v[(!).(nofillmask)] = FILL
    end
    # no FILL inside glaciermask
    @assert !any(dem.v[glaciermask].==FILL)
    @assert !any(isnan.(dem.v[glaciermask]))

    # 2D slopes
    ret_nans = false
    alpha2d = absslope(dem, glaciermask, ret_nans)

    bands, bandi = bin_grid(dem, binsize_or_bins, glaciermask,
                            binround=binround, min_bin_number_ends=min_bin_number_ends,
                            min_bands=min_bands)
    @assert length(bands)>=min_bands "Need at least $min_bands elevation bins, only got $(length(bands))"

    nb = length(bands)
    cellsz = step(dem.x)^2

    totalarea = sum(glaciermask)*cellsz

    malphas, widths, lengths, areas = (zeros(nb) for i=1:4)
    dzs = zeros(nb)
    for i=1:nb
        ind = bandi[i]
        if i!=nb
            dzs[i] = abs(bands[i+1]-bands[i])
        else
            # TODO: a hack. Fix
            dzs[i] = abs(bands[i]-bands[i-1])
        end
        # this is the critical step:
        malphas[i] = band_slope!(alpha2d[ind], i, alpha_min, alpha_max)
        areas[i] = length(ind)*cellsz
        lengths[i] = dzs[i]/tan(malphas[i])
        widths[i] = areas[i]/lengths[i]
    end
    # update missing bands
    for i=1:nb
        if malphas[i]==-9999
            ma1 = malphas[max(1,i-1)]
            ma2 = malphas[min(nb,i+1)]
            if ma1==-9999 && ma2==-9999
                error("Too many consecutive bands for which no slope could be calculated: $(max(1,i-1):min(nb,i+1))")
            elseif ma1==-9999
                malphas[i] = ma2
            elseif ma2==-9999
                malphas[i] = ma1
            else
                malphas[i] = 1/2*(ma1+ma2)
            end
            lengths[i] = dzs[i]/tan(malphas[i])
            widths[i] = areas[i]/lengths[i]
        end
    end

    # Smooth the width:
    # Note, this can make the malphas noisy!
    if window_width_smooth>0
        widths = boxcar(widths, round(Int, window_width_smooth/mean(dzs) ))
        for i=1:nb
            # make sure length is zero when area is zero
            lengths[i] = areas[i]==widths[i]==0 ? zero(lengths[i]) : areas[i]/widths[i]
            malphas[i] = atan(dzs[i]/lengths[i])
            @assert malphas[i]>=0
        end
    end

    x = vcat(0,cumsum(lengths))
    xmid = x[1:end-1] + diff(x)/2

    # tests
    # if abs(totalarea-sum(areas))>1.0
    #     error("Something's amiss, sum of area of bands $(sum(areas)) not equal total area $totalarea,")
    # end
    # check band length against diagonal
    tmp = sum(glaciermask,2)
    xextent = (findlast(tmp.>0)-findfirst(tmp.>0))*dx
    tmp = sum(glaciermask,1)
    yextent = (findlast(tmp.>0)-findfirst(tmp.>0))*dx
    box_diag = sqrt(xextent^2 + yextent^2)
    if verbose && abs(sum(lengths)-box_diag)/box_diag>0.4
        warn("Glacier length from might be wrong. Band-length: $(sum(lengths)/1e3)km, bounding-box diag: $(box_diag/1e3)km")
    end

    return bands, bandi, malphas, areas, lengths, widths, x, xmid, dem, alpha2d
end

"""
    band_slope!(alphas, bandnr, alpha_min=deg2rad(0.4), alpha_max=deg2rad(60.0))

Need to calculate a meaningful mean of the slopes in a elevation bin.
*This is tricky but critical!*

One check can be that all the bin-lengths should add up to the total
glacier length.

Input:
- alphas -- slope angles in one band (these are sorted in place, thus the ! in the function name)
- bandnr -- which band those alphas belong to (only used for error message)
- alpha_max, alpha_min -- maximal and minimal allowed slope
"""
function band_slope!(alphas, bandnr, alpha_min, alpha_max)
    # parameters
    ratio_fac = 2
    f_q5 = 0.05
    f_q20 = 0.2
    f_q80 = 0.8
    q_band = (0.55, 0.95)

    n = length(alphas)
    if n==0
        return -9999*one(alpha_min)
        # error("Band $bandnr has no elements!")
        # return deg2rad(45)
    end
    # magic slope calculation
    sort!(alphas)
    # calculate indices of 5th, 20th and 80th percentiles:
    iq5  = max(1,round(Int,n*f_q5))
    iq20 = max(1,round(Int,n*f_q20))
    iq80 = min(n,round(Int,n*f_q80))
    # angle of those quantiles:
    q5, q20, q80 = [max(rad2deg.(i),eps(alpha_min)) for i in (alphas[iq5], alphas[iq20], alphas[iq80])]
    # Now some of Matthias' magic:
    a = (q20/q80)*ratio_fac # 2x ratio of angles
    a = min(a, q_band[2])
    a = max(a, q_band[1]) # scaled to lie [0.55, 0.95]
    iq_magic = round(Int,n*a)  # set a "new" magic quantile to that value
    q_magic = rad2deg(alphas[iq_magic])
    # only use indices within those quantiles
    ind = q5 .<= rad2deg.(alphas) .< q_magic
    out = sum(ind)>1 ? mean(alphas[ind]) : mean(alphas)
    # limit alphas
    out = max(alpha_min, out)
    out = min(alpha_max, out)
    return out
end

"""
Bins a gird into bands.  Often used to bin a DEM into elevation bands.

- g -- to be binned ::Gridded or ::Matrix
- binsize_or_bins -- bin size or the bins (a Range)
- mask -- specify if not all locations of a gird should be binned.

KW:
- binround -- floor the bin-start using this many digits (see help of floor)
- min_bin_number_ends -- minimum number of elements in the uppermost and lowermost
                         band.  If below, merge those cells into the first viable band.
- min_bands -- minimum number of bands produced (3).  Not integrated with min_bin_number_ends!

Return:
- bands -- a range of the bands, e.g. 0.0:10.0:100.0
- bandi -- a Vector{Vector{Int}} of length(bands) with each element
           containing the indices of cells in the band
"""
function bin_grid(g::Gridded, binsize_or_bins, mask=BitArray([]);
                  binround=_binround(binsize_or_bins), min_bin_number_ends=0, min_bands=4)
    bin_grid(g.v, binsize_or_bins, mask;
             binround=binround, min_bin_number_ends=min_bin_number_ends,
             min_bands=min_bands)
end
function bin_grid(v::Matrix, binsize::Number, mask=BitArray([]);
                  binround=_binround(binsize), min_bin_number_ends=0, min_bands=4)/
    if isempty(mask)
        v = v
        ginds = 1:length(v)
    else
        @assert size(mask)==size(v)
        v = v[mask]
        ginds = find(mask[:])
    end
    nv = length(v)
    bins = 1.0:-1
    while length(bins)<=min_bands
        mi, ma = minimum(v), maximum(v)
        if binsize>=0
            binstart = floor(mi, binround)
            binend = floor(ma, binround) # better: `ceil(ma, binround) - binsize` ?
        else
            binstart = ceil(ma, binround)
            binend = ceil(mi, binround) # better: `ceil(ma, binround) - binsize` ?
        end
        @assert !isnan(binstart) && !isnan(binend)
        bins = binstart:binsize:binend # these are the start of the bins

        # decreas binsize for next round
        binsize = step(bins)/2
    end
    return _bin_grid_kernel(bins, nv, v, ginds, min_bin_number_ends)
end
function bin_grid(v::Matrix, bins::AbstractRange, mask=BitArray([]); min_bin_number_ends=0, kw...)
    if isempty(mask)
        v = v
        ginds = 1:length(v)
    else
        @assert size(mask)==size(v)
        v = v[mask]
        ginds = find(mask[:])
    end
    nv = length(v)
    _bin_grid_kernel(bins, nv, v, ginds, min_bin_number_ends)
end

@inbounds function _bin_grid_kernel(bins, nv, v, ginds, min_bin_number_ends)
    # initialize output
    bandi = Vector{Int}[]
    for b in bins
        ind = Int[]
        push!(bandi, ind)
    end
    # fill it
    for j=1:nv
        if step(bins)>0
            i = searchsortedlast(bins, v[j])
            i = i==0 ? 1 : i # if smaller then add to lowest bin
        else
            # https://github.com/JuliaLang/julia/issues/18653
            i = searchsortedlast(collect(bins), v[j], rev=true)
            i = i==0 ? 1 : i # if smaller then add to highest bin (i.e. bins[1])
        end
        push!(bandi[i], ginds[j])
    end
    # remove top and bottom bins if too few elements
    inds2pop = [1,length(bandi)] # remove up to and including these bandi
    num = [0,0]
    for i = 1:length(bandi)
        inds2pop[1] = i
        if length(bandi[i])+num[1]>=min_bin_number_ends
            break
        end
        num[1] +=length(bandi[i])
    end
    for i = length(bandi):-1:1
        inds2pop[2] = i
        if length(bandi[i])+num[2]>=min_bin_number_ends
            break
        end
        num[2] +=length(bandi[i])
    end
    # add the dropped cells to the next/previous band
    append!(bandi[inds2pop[1]], vcat(bandi[1:inds2pop[1]-1]...))
    append!(bandi[inds2pop[2]], vcat(bandi[inds2pop[2]+1:end]...))
    return bins[inds2pop[1]:inds2pop[2]], bandi[inds2pop[1]:inds2pop[2]]
end

import Interpolations
"""
Bins a trajectory using a grid

- tr -- trajectory to be binned
- g -- g.v of grid used for binning
- binsize_or_bins -- bin size or the bins (a Range)
- mask -- specify if not all locations of a gird should be binned.]

KW:
- binround -- floor the bin-start using this many digits (see help of floor)
"""
function bin_traj(tr::Traj, g::Gridded, binsize_or_bins, mask=trues(size(g.v)); binround=_binround(binsize_or_bins))
    @assert size(mask)==size(g.v)

    demi  = Interpolations.interpolate((g.x, g.y), g.v, Interpolations.Gridded(Interpolations.Linear()) )
    maski  = Interpolations.interpolate((g.x, g.y), mask, Interpolations.Gridded(Interpolations.Constant()) ) # returns Int!

    v = [demi[x,y] for (x,y) in zip(tr.x,tr.y)]
    vm = Bool[maski[x,y] for (x,y) in zip(tr.x,tr.y)]
    v = v[vm]
    ginds = find(vm)
    nv = length(v)

    if isa(binsize_or_bins, Number)
        mi, ma = minimum(v), maximum(v)
        binstart = floor(mi, binround)
        binend = floor(ma, binround)
        bins = binstart:binsize_or_bins:binend # these are the start of the bins
    else
        bins = binsize_or_bins
    end
    _bin_grid_kernel(bins, nv, v, ginds, 0)
end

"""
   map_onto_bands(bandi, field, fn=mean)

Map a field onto the (elevation) bands.  The field needs to have the
same size as the original binned-grid.

Input:
- bandi -- as returned by bin_grid
- field -- the field, either a Matrix or a Gridded
- fn -- the function to do the reduction with.  Default==mean
- fill -- fill value, if set, ignore those points

Output
- the value of the field on the bands.  If no values were found in a band,
  then return NaN.
"""
function map_onto_bands(bandi, field::Matrix, fn=mean, fill=nothing)
    resT = typeof(fn(field[bandi[1]])) # to get result type
    out = zeros(resT, length(bandi))
    for i in 1:length(bandi)
        count = 0
        for j=1:length(bandi[i])
            val = field[bandi[i][j]]
            if val!=fill
                # only use the ones which have no fill
                out[i] += val
                count+=1
            end
        end
        out[i] /=count
    end
    return out
end
map_onto_bands(bandi, field::Gridded, fn=mean, fill=NaN) = map_onto_bands(bandi, field.v, mean, fill)

"""
    map_back_to_2D(dims2d, bandi, field1d)

Maps 1D field back onto 2D.  More or less inverse of map_onto_bands.
"""
function map_back_to_2D(dims2d, bandi, field1d)
    out = zeros(eltype(field1d), dims2d)
    map_back_to_2D!(out, bandi, field1d)
    out
end
"""
    map_back_to_2D!(out2d, bandi, field1d)

Maps 1D field back onto 2D.  More or less inverse of map_onto_bands.
"""
function map_back_to_2D!(out, bandi, field1d)
    for (i,is) in enumerate(bandi)
        out[is] = field1d[i]
    end
    nothing
end

"""
    bins2matrix(g::Union{Gridded,Matrix}, bands, bandi) -> binmat

Return a matrix which gives the bin-number of each its (i,j) locations.
Locations not binned (i.e. masked) are ==0.
"""
bins2matrix(g::Gridded, bands, bandi) = bins2matrix(g.v, bands, bandi)
function bins2matrix(g::Matrix, bands, bandi)
    out = zeros(Int, size(g))
    for (n,b) in enumerate(bandi)
        for i in b
            out[i] = n
        end
    end
    return out
end

"""
    bandi_for_other_grid(bands, bandi, g::Gridded,
                         othergrid::Gridded, othermask=trues(size(othergrid.v))
    bandi_for_other_grid(bands, bandi, g::Gridded,
                         othergrid::Gridded, othermask=trues(size(othergrid.v)))

Returns vector of indices (bandi) to map a different grid (othergird) onto the
bands encoded in `binmat` (or `bands, bandi`) and grid `g`.  It only
maps points onto bands which are within the mask applied to generate
the bands.  Additionally & optionally, a mask for the othergrid can
also be given.  Return:

    bandi
"""
function bandi_for_other_grid(bands, bandi::Vector{Vector{Int}}, g::Gridded,
                              othergrid::Gridded, othermask=trues(size(othergrid.v)))
    binmat=bins2matrix(g, bands, bandi)
    bandi_for_other_grid(bands, bandi, binmat, g, othergrid, othermask)
end
function bandi_for_other_grid(bands, bandi, binmat::Matrix{Int}, g::Gridded,
                              othergrid::Gridded, othermask=trues(size(othergrid.v)) )
    og = othergrid
    @assert size(othergrid)==size(othermask)
    if g.x!=og.x || g.y!=og.y
        bandi_ = [Int[] for i=1:length(bands)]
        s2i = LinearIndices(og.v)
        itpm = Interpolations.interpolate((g.x, g.y), binmat,
                            Interpolations.Gridded(Interpolations.Constant()) );
        itpm = Interpolations.extrapolate(itpm, 0);
        for j=1:size(og.v,2)
            for i=1:size(og.v,1)
                if othermask[i,j]
                    ind = convert(Int, itpm(og.x[i], og.y[j]))
                    if ind>0
                        push!(bandi_[ind], s2i[i, j])
                    end
                end
            end
        end
    else
        bandi_ = deepcopy(bandi)
    end
    return bandi_
end

#####################
# Extrapolation of IV
#####################

"""
To specify which edge of a cell is meant.  `_noedge` can be used if none
is used.
"""
@enum Loc _noedge=0 left=1 right=2 lower=3 upper=4

"Orientation of a line"
@enum Orientation nohand=0 lefthand=1 righthand=2

"""
One cell-edge of a regular gird of cells.  Specified as a cell
and which edge of the four edges of a cell.

TODO:
Arguably not the best datastructure for what is done below.
"""
struct Edge
    i::Int # cell ind
    j::Int # cell ind
    loc::Loc
end
const noedge = Edge(-1,-1,_noedge)
Base.show(io::IO,e::Edge) = println(io, "($(e.i),$(e.j),$(e.loc))")

"""
    get_nodes(e::Edge)

Returns the start and end nodes of the edge on the staggered grid.

The nodes are returned such that the cell is on the right of the edge.
"""
function get_nodes(e::Edge)
    i,j,loc = e.i, e.j, e.loc
    if loc==left
        return (i,j), (i,j+1)
    elseif loc==right
        return (i+1,j+1), (i+1,j)
    elseif loc==lower
        return (i+1,j), (i,j)
    elseif loc==upper
        return (i,j+1), (i+1,j+1)
    else
        error("!")
    end
end

"""
    orientation(e1::Edge) -> always left
    orientation(e1::Edge, e2::Edge)

Return whether the cells are on left, right or both of the two cells.
If not connected throws an error or return `_noedge`, depending on
`throwerror` flag.

"""
function orientation(e1::Edge, e2::Edge, throwerror=true)::Orientation
    n1,n2 = get_nodes(e1)
    m1,m2 = get_nodes(e2)
    if n2==m1
        return righthand
    elseif n1==m2
        return lefthand
    elseif n1==m1 || n2==m2
        return nohand
    elseif throwerror
        error("The two edges:\n $e1 $e2 are not connected")
    else
        return nohand
    end
end
orientation(e1::Edge)::Orientation = lefthand

"""
A line made up of a continuous (and sorted) collection of edges.  Note
that the edges themselves also have a cell associated with them.  So,
a Line also a collection of (possibly repeated) cells.

The line has an orientation: when traversing the Line from 1 to n, all
cells are on the right of their edge.
"""
struct Line
    edges::Vector{Edge}
end
#Base.show(io::IO,l::Line) = show(io, l.edges)

Base.getindex(l::Line,i) = l.edges[i]
Base.length(l::Line) = length(l.edges)

"""
    orientation(l::Line)

Returns the orientation of a line: `left` means cell is on left of line.
"""
orientation(l::Line)::Orientation = orientation(l.edges)
function orientation(l::Vector)::Orientation
    if length(l)==0
        error("Line has length 0")
    end
    if length(l)==1
        # println("Line has length $(length(l)).  Too short to establish orientation.")
        return orientation(l[1]) # returns `left`
    end
    # establish orientation of first two segments
    ori = orientation(l[1],l[2])
    for i=3:length(l)
        ori2 = orientation(l[i-1],l[i])
        if ori2!=ori
            println("""
                    Line changes orientation between edge # $(i-2) and $i from $ori to $ori2.
                    e0: $(l[i-2])
                    e1: $(l[i-1])
                    e2: $(l[i])
                    """)
            return _noedge
        end
    end
    return ori
end

"""
    next_edge!(edges::Set{Edge}, edge::Edge, testori::Orientation)

Return an adjacent edge of input edge such that the orientation of the
two edges is equal `testori`.  Pops returned edge off `edges`.  If no
edge is found, return `Edge(-1,-1,nodege)`.

"""
function next_edge!(edges::Set{Edge}, edge::Edge, testori::Orientation)
    edge.loc==_noedge && error("Not a proper edge: $edge")
    locs = (left,right,lower,upper)
    for i=[0,-1,1],j=[0,-1,1]
        for loc in locs
            test_edge = Edge(edge.i+i,edge.j+j,loc)
            if orientation(edge, test_edge,false)==testori && (test_edge in edges)
                return pop!(edges, test_edge)
            end
        end
    end
    # nothing found
    return noedge
end

"""
    get_ux_uy(dem, mask)

The direction and magnitude of steepest descent on each point within the mask.
"""
function get_ux_uy(dem, mask)
    ux,uy = (-).(gradient3by3(dem, mask))
    return ux, uy
end

"""
    get_edges_on_boundary(bands, bandi, binmat, landmask=nothing) -> es

Return a list of sets of cell-edges which are at the boundary.

    typeof(es) == Dict{Tuple{Int,Int},Set{Edge}}()

which maps (bandnr,otherbandnr) => Set of edges

"""
function get_cells_on_boundary(bands, bandi, binmat, landmask=nothing)
    error("Not updated to Julia 1.0 yet")
    dims = size(binmat)
    i2s = CartesianIndices(binmat)

    if landmask!=nothing
        # encode sea cells in binmat
        binmat = copy(binmat)
        binmat[(landmask.==0) .& (binmat.==0)] = -1
    end

    # The boundaries of all bands:
    # (bandnr,otherbandnr) => Set of edges
    edges_at_boundaries = Dict{Tuple{Int,Int},Set{Edge}}()

    # Loop to identify all boundary cells of each band and determine
    # which cell-edges are on the boundary:
    for (ib,bup) in enumerate(bands)
        for I in bandi[ib]
            i,j = i2s[I].I
            loc = 1
            # left-right
            for ii=-1:2:1
                if 1<=i+ii<=dims[1] # skip if outside of binmat
                    iother = binmat[i+ii,j]
                    if iother!=ib
                        push!(get!(edges_at_boundaries, (ib, iother), Set{Edge}()),
                              Edge(i,j,loc))
                    end
                end
                loc+=1
            end
            # lower-upper
            for jj=-1:2:1
                if 1<=j+jj<=dims[2] # skip if outside of binmat
                    iother = binmat[i,j+jj]
                    if iother!=ib
                        push!(get!(edges_at_boundaries, (ib, iother), Set{Edge}()),
                              Edge(i,j,loc))
                    end
                end
                loc+=1
            end
        end
    end
    return edges_at_boundaries
end

"""
    calc_boundaries(bands, bandi, binmat, landmask=nothing)

Calculate interface between bands and also to outside (bandnr==0) and
sea-cells (bandnr==-1).

Returns `boundaries` which is `bandnr => otherbandnr => Vector{Line}`
i.e. the collection of all boundary lines for each
`(bandnr,otherbandnr)`.  Each line has the elevation band on its
*left*.  For cells outside the glacier `otherbandnr==0` except if
'landmask` is given, then sea-cells have `otherbandnf==-1`.

Of type `Vector{Dict{Int,Vector{Line}}}()`.
"""
function calc_boundaries(bands, bandi, binmat, landmask=nothing)
    dims = size(binmat)

    if landmask!=nothing
        # encode sea cells in binmat
        binmat = copy(binmat)
        binmat[(landmask.==0) .& (binmat.==0)] = -1
    end

    # The boundaries of all bands:
    # (bandnr,otherbandnr) => Set of edges
    edges_at_boundaries = get_cells_on_boundary(bands, bandi, binmat, landmask)

    # bandnr => otherbandnr => Vector{Line}
    boundaries = [Dict{Int,Vector{Line}}() for i=1:length(bands)]

    for ((band,otherband), edges) in edges_at_boundaries

        # # delete:
        # band, otherband = 2,0
        # edges = boundaries[(band,otherband)]

        output = Line[]

        # Take one starting point and go in both directions until no
        # more neighbors found.
        while length(edges)>0
            firstedge = first(edges)
            delete!(edges, firstedge)
            out = [Edge[], Edge[]]
            for ii=1:2 # first lefthand Orientation(1), then righthand Orientation(2)
                curedge = firstedge
                while curedge!=noedge
                    push!(out[ii], curedge)
                    curedge = next_edge!(edges, curedge, Orientation(ii))
                end
            end
            prepend!(out[2], reverse(out[1][2:end]))
            # make sure the elevation band is on the line's left:
            ori = orientation(out[2])
            if ori == righthand
                push!(output, Line(reverse(out[2])))
            elseif ori==lefthand
                push!(output, Line(out[2]))
            elseif ori==nohand
                error("""
                      Constructed line has not a consistent orientation.
                      Band: $band
                      Otherband: $otherband
                      """
                      )
            end
        end
        boundaries[band][otherband] = output
    end

    return boundaries
end

"""
    calc_fluxdir(l::Line, ux, uy, window, dx, bands) -> flux,ii,jj,fluxdirx,fluxdiry

Calculated the flux-direction on a Line by taking a running average of (ux,uy)
over some part of it of length window.  Also returns the flux to use and the upstream
cell `(ii,jj)`.  The flux across the k-th edge is then given by:

    flux[k] * h[ii[k],jj[k]] * u[ii[k],jj[k]]

where h is the thickness and u the depth averaged velocity.

The flux across the line is

    sum(flux .* h_line .* u_line)

Input:
- l::Line
- ux,uy: unsmoothed velocity field (calculated as fall-line of surface DEM with get_ux_uy)
- window : length of smoothing window (along band) in number of edges
- dx : grid spacing
- bands::AbstractRange : the elevation bands

Output:
- flux : the flux
- ii, jj : indices of upstream cell
- fluxdirx,fluxdiry : flux direction at upstream cell ii, jj

Notes:
- this is not 100% correct but probably close enough.
- now the averaging is over all edges.  Arguably it could be done over cells?
"""
function calc_fluxdir(l::Line, ux, uy, window::Int, dx, bands)
    edges = l.edges
    ne = length(edges)
    # flux in x and y-direction on the cell associated with a Line edge:
    fluxdirx = zeros(ne)
    fluxdiry = zeros(ne)
    # the flux which goes with the edge:
    flux = zeros(ne)
    # indices of upstream cell
    ii = zeros(Int,ne)
    jj = zeros(Int,ne)

    R = CartesianIndices(size(edges))
    I1, Iend = first(R), last(R)
    for I in R
        fx, fy = 0.0, 0.0
        for J in CartesianIndices(max(I1, I-I1*window), min(Iend, I+I1*window))
            i,j,loc = edges[J].i, edges[J].j, edges[J].loc
            # average ux,uy in the two cells
            if loc==left
                fx += (ux[i-1,j] + ux[i,j])
                fy += (uy[i-1,j] + uy[i,j])
            elseif loc==right
                fx += (ux[i+1,j] + ux[i,j])
                fy += (uy[i+1,j] + uy[i,j])
            elseif loc==lower
                fx += (ux[i,j-1] + ux[i,j])
                fy += (uy[i,j-1] + uy[i,j])
            elseif loc==upper
                fx += (ux[i,j+1] + ux[i,j])
                fy += (uy[i,j+1] + uy[i,j])
            else
                error()
            end
        end
        norm = sqrt(fx^2+fy^2)
        fx, fy = fx/norm, fy/norm
        fluxdirx[I] = fx
        fluxdiry[I] = fy
        # figure out upstream cell
        i,j,loc = edges[I].i, edges[I].j, edges[I].loc
        ii[I],jj[I],flux[I] =
            if loc==left
                fx<0 ? (i,j,fx*dx) : (i-1,j,fx*dx)
            elseif loc==right
                fx<0 ? (i+1,j,-fx*dx) : (i,j,-fx*dx)
            elseif loc==lower
                fy<0 ? (i,j,fy*dx) : (i,j-1,fy*dx)
            else #if loc==upper
                fy<0 ? (i,j+1,-fy*dx) : (i,j,-fy*dx)
            end
    end
    # for reverse ordered bands, need a sign flip
    sig = sign(bands[2]-bands[1])

    return sig*flux,ii,jj,fluxdirx,fluxdiry
end

"""
    calc_u(q1d, boundaries, u_trial, thick,
                ux, uy, dx, mask, bands, bandi, lengths,
                flux_dir_window, # in [m]
                boxcarM::AbstractMatrix or window_frac,
                scale_u=ones(q1d);
                plotyes=false,
                x=nothing, y=nothing)

Calculates a 2D field of depth averaged ice flow speed `ubar` which has the
same flux across elevation band boundaries as the supplied 1D flux
`q1d`.  Only `q1d` is a flow-band variable.  The flux across elevation
bands is calculated with `calc_fluxdir`.

Input:
- q1d: 1d flux (m/a) on elevation bands
- boundaries: from calc_boundaries
- u_trial: 2D velocity field which has the right shape (say cross-valley) but not the right magnitude
- thick: ice thickness
- ux,uy: unsmoothed velocity field (calculated as fall-line of surface DEM with get_ux_uy)
- dx: grid size
- mask: true where glacier is
- bands : elevation bands
- bandi
- lengths: length of each elevation band
- flux_dir_window: as used in calc_fluxdir
- boxcarM or window_frac:
  - pre-calculated boxcar operator
  - window_frac: window over which to smooth as fraction of the maximum(length)
- scale_u (ones(q1d)): if provided scale output by this much -> so make surface speeds

Output:
- ubar2d: the smoothed IV
- ubar: the IV at the elevation-band boundaries
- mask_ubar: true where ubar contains a value
- facs: scaling factor in 1D

Note:
- An assumption about the distribution of the u is needed, which will
  be scaled to conform to mass conservation.  At the moment a function
  `thick[i,j].^u_exp` is used.
"""
function calc_u(q1d, boundaries, u_trial, thick,
                ux, uy, dx, mask, bands, bandi, lengths,
                flux_dir_window, # in [m]
                filter, scale_u=ones(q1d);
                plotyes=false,
                x=nothing, y=nothing)
    u, mask_u, facs = _calc_u(q1d, boundaries, u_trial, thick,
                                    ux, uy, dx, mask, bands,
                                    flux_dir_window, # in [m]
                                    plotyes, x, y, scale_u)
    u_ = copy(u)
    u_[isnan.(u_)] = 0  # remove NaNs to allow filtering

    # type assertion is needed for type-stability
    out::Array{eltype(q1d),2} = if filter isa AbstractMatrix
        VAWTools.apply_boxcar_matrix(filter, u_)
    else
        boxcar(u_, Int((filter*maximum(lengths))÷dx)+1, mask_u, (!).(mask) )
    end

    # Make depth-averaged flow speed into a surface flow speed.
    scale_u2d = map_back_to_2D(size(out), bandi, scale_u)
    out = out .* scale_u2d

    return out, u, mask_u, facs
end
# this helper function is needed for type stability.
function _calc_u(q1d, boundaries, u_trial, thick,
                 ux, uy, dx, mask, bands,
                 flux_dir_window,
                 plotyes,
                 x, y, scale_u) # these are only needed for plotting
    #plotyes && figure()
    dims = size(mask)

    # calculate the u at all elevation band boundaries:
    u2d = zeros(dims)*NaN
    facs = Float64[] # scaling factor in 1D
    for (ib,bnd) in enumerate(boundaries)
        if ib<length(bands)
            # Find the receiving band, can be a number further than 1.
            ibb = 0
            for ib_=ib+1:length(bands)
                if haskey(bnd, ib_)
                    ibb = ib_
                    break
                end
            end
            ibb==0 && error("No band below band $ib, but $ib is not bottom band!  This means that the domain is likely disjoint.")
            bb=bnd[ibb]
        else # outflow at terminus
            # This probably needs updating where several elevation bands contribute (tide-water)
            # -> no, only ever the last band does outflow in our 1D model
            #
            # This errors at times, for example on RGI60-14.11814 (which is not sea-terminating)
            bb = get(bnd, -1, bnd[0]) # if sea-terminating use those edges.
        end
        # loop over segments
        ffs = Float64[]
        is = Int[]
        js = Int[]
        for l in bb
            #@assert orientation(l)
            ff, fi, fj, fx, fy = calc_fluxdir(l, ux, uy, Int(flux_dir_window÷dx), dx, bands)
            # if plotyes
            #     quiver(x[fi],y[fj],fx,fy)
            # end
            append!(is, fi)
            append!(js, fj)
            append!(ffs, ff)
        end
        u_t = [u_trial[i,j] for (i,j) in zip(is,js)]
        h_line = [thick[i,j] for (i,j) in zip(is,js)]
        q_trial = sum(ffs.*u_t.*h_line)
        fac = max(q1d[ib],0)/q_trial
        push!(facs, fac)
        u = u_t*fac
        # m,ij = findmin(u)
        #@assert all(u.>=0) "$i, $ij, $(u_trial[ij]), $(ff[ij]), $(q_trial), $(q1d[i]), $(h_line[ij])"
        for (n,(i,j)) in enumerate(zip(is,js))
            if u[n]<0;
                println("Flow speed<0: $(u[n]) at location ($i,$j).  This should not happen!  Setting to zero.")
                u[n]=0
            end
            # also scale u:
            u2d[i,j] = u[n] * scale_u[ib]
        end
    end
    mask_u = mask .& ((!).(isnan.(u2d))) # location of all cells for which `u` was calculated
    # if plotyes
    #     # imshow(binmat',origin="lower", extent=(x[1],x[end],y[1],y[end]), cmap="flag"); colorbar();
    #     imshow(u',origin="lower", extent=(x[1],x[end],y[1],y[end]),); colorbar(); clim(0,50)
    # end

    return u2d, mask_u, facs
end

"""
    get_iv_boxcar_M(F, dem, mask, bands, bandi, lengths, iv_window_frac)

Precalculate the boxcar operator for the IV calculation (this is the
most expensive part).
"""
function get_iv_boxcar_M(F, dem, mask, bands, bandi, lengths, iv_window_frac)
    ux,uy = get_ux_uy(dem, mask)
    binmat = bins2matrix(dem, bands, bandi)
    boundaries = calc_boundaries(bands,bandi,binmat)
    q1d = ones(bands)
    u_trial = ones(dem.v)
    thick = u_trial
    dx = step(dem.x)
    flux_dir_window = 2
    mask_u = _calc_u(q1d, boundaries, u_trial, thick,
                     ux, uy, dx, mask,
                     bands,
                     flux_dir_window,
                     false,nothing,nothing,ones(q1d))[2]

    return boxcar_matrix(F, Int((iv_window_frac*maximum(lengths))÷dx)+1, mask_u, (!).(mask)),
           boundaries, ux, uy
end

#################
# Plotting
################
"""
    plot_bands(dem, bands, bandi; bands2plot=1:length(bands))

Plots the bands in 2D.  For length(bands2plot)<=16 the default
colorscale will show one color per band.

Needs Plots.jl imported in Main (REPL)
"""
function plot_bands(dem, bands, bandi; bands2plot=1:length(bands))
    binmat = Float64.(VAWTools.bins2matrix(dem, bands, bandi))
    for i in eachindex(binmat)
        if !(binmat[i] in bands2plot)
            binmat[i] = NaN
        end
    end
    Main.PyPlot.contourf(dem.x,dem.y,binmat',aspect_ratio=:equal)
end


text:
using JsonGrinder, JSON, Test, SparseArrays, Flux
using JsonGrinder: ExtractScalar, ExtractCategorical, ExtractArray, ExtractBranch
using Mill: catobs
using LinearAlgebra

@testset "Testing scalar conversion" begin
	sc = ExtractScalar(Float64,2,3)
	@test all(sc("5").data .== [9])
	@test all(sc(5).data .== [9])
	@test all(sc(nothing).data .== [0])
end

@testset "Testing array conversion" begin
	sc = ExtractArray(ExtractCategorical(2:4))
	@test all(sc([2,3,4]).data.data .== Matrix(1.0I, 4, 3))
	@test ismissing(sc(nothing).data)
	@test all(sc(nothing).bags.bags .== [0:-1])
	sc = ExtractArray(ExtractScalar(Float64))
	@test all(sc([2,3,4]).data.data .== [2 3 4])
	@test ismissing(sc(nothing).data)
	@test all(sc(nothing).bags.bags .== [0:-1])
end

@testset "Testing ExtractBranch" begin
	vector = Dict("a" => ExtractScalar(Float64,2,3),"b" => ExtractScalar(Float64));
	other = Dict("c" => ExtractArray(ExtractScalar(Float64,2,3)));
	br = ExtractBranch(vector,other)
	a1 = br(Dict("a" => 5, "b" => 7, "c" => [1,2,3,4]))
	a2 = br(Dict("a" => 5, "b" => 7))
	a3 = br(Dict("a" => 5, "c" => [1,2,3,4]))
	@test all(catobs(a1,a1).data[1].data .==[7 7; 9 9])
	@test all(catobs(a1,a1).data[2].data.data .== [-3 0 3 6 -3 0 3 6])
	@test all(catobs(a1,a1).data[2].bags .== [1:4,5:8])

	@test all(catobs(a1,a2).data[1].data .==[7 7; 9 9])
	@test all(catobs(a1,a2).data[2].data.data .== [-3 0 3 6])
	@test all(catobs(a1,a2).data[2].bags .== [1:4,0:-1])

	@test all(catobs(a2,a3).data[1].data .==[7 0; 9 9])
	@test all(catobs(a2,a3).data[2].data.data .== [-3 0 3 6])
	@test all(catobs(a2,a3).data[2].bags .== [0:-1,1:4])

	@test all(catobs(a1,a3).data[1].data .==[7 0; 9 9])
	@test all(catobs(a1,a3).data[2].data.data .== [-3 0 3 6 -3 0 3 6])
	@test all(catobs(a1,a3).data[2].bags .== [1:4,5:8])

	br = ExtractBranch(vector,nothing)
	a1 = br(Dict("a" => 5, "b" => 7, "c" => [1,2,3,4]))
	a2 = br(Dict("a" => 5, "b" => 7))
	a3 = br(Dict("a" => 5, "c" => [1,2,3,4]))
	@test all(a1.data .==[7; 9])
	@test all(a2.data .==[7; 9])
	@test all(a3.data .==[0; 9])

	br = ExtractBranch(nothing,other)
	a1 = br(Dict("a" => 5, "b" => 7, "c" => [1,2,3,4]))
	a2 = br(Dict("a" => 5, "b" => 7))
	a3 = br(Dict("a" => 5, "c" => [1,2,3,4]))

	@test all(a1.data.data .== [-3 0 3 6])
	@test all(a1.bags .== [1:4])
	@test all(catobs(a1,a1).data.data .== [-3 0 3 6 -3 0 3 6])
	@test all(catobs(a1,a1).bags .== [1:4,5:8])

	@test all(catobs(a1,a2).data.data .== [-3 0 3 6])
	@test all(catobs(a1,a2).bags .== [1:4,0:-1])

	@test all(a3.data.data .== [-3 0 3 6])
	@test all(a3.bags .== [1:4])
	@test all(catobs(a3,a3).data.data .== [-3 0 3 6 -3 0 3 6])
	@test all(catobs(a3,a3).bags .== [1:4,5:8])
end

@testset "Testing Nested Missing Arrays" begin
	other = Dict("a" => ExtractArray(ExtractScalar(Float64,2,3)),"b" => ExtractArray(ExtractScalar(Float64,2,3)));
	br = ExtractBranch(nothing,other)
	a1 = br(Dict("a" => [1,2,3], "b" => [1,2,3,4]))
	a2 = br(Dict("b" => [2,3,4]))
	a3 = br(Dict("a" => [2,3,4]))
	a4 = br(Dict{String,Any}())

	@test all(catobs(a1,a2).data[1].data.data .== [-3.0  0.0  3.0  6.0  0.0  3.0  6.0])
	@test all(catobs(a1,a2).data[1].bags .== [1:4, 5:7])
	@test all(catobs(a1,a2).data[2].data.data .== [-3.0  0.0  3.0])
	@test all(catobs(a1,a2).data[2].bags .== [1:3, 0:-1])

	@test all(catobs(a2,a3).data[1].data.data .== [0.0  3.0  6.0])
	@test all(catobs(a2,a3).data[1].bags .== [1:3, 0:-1])
	@test all(catobs(a2,a3).data[2].data.data .== [0 3 6])
	@test all(catobs(a2,a3).data[2].bags .== [0:-1, 1:3])

	@test all(catobs(a1,a4).data[1].data.data .== [-3.0  0.0  3.0  6.0])
	@test all(catobs(a1,a4).data[1].bags .== [1:4, 0:-1])
	@test all(catobs(a1,a4).data[2].data.data .== [-3.0  0.0  3.0])
	@test all(catobs(a1,a4).data[2].bags .== [1:3, 0:-1])
end

@testset "ExtractOneHot" begin
	samples = ["{\"name\": \"a\", \"count\" : 1}",
		"{\"name\": \"b\", \"count\" : 2}",]
	vs = JSON.parse.(samples)

	e = ExtractOneHot(["a","b"], "name", "count")
	@test e(vs).data[:] ≈ [1, 2, 0]
	@test e(nothing).data[:] ≈ [0, 0, 0]
	@test typeof(e(vs).data) == SparseMatrixCSC{Float32,Int64}
	@test typeof(e(nothing).data) == SparseMatrixCSC{Float32,Int64}

	e = ExtractOneHot(["a","b"], "name", nothing)
	@test e(vs).data[:] ≈ [1, 1, 0]
	@test e(nothing).data[:] ≈ [0, 0, 0]
	@test typeof(e(vs).data) == SparseMatrixCSC{Float32,Int64}
	@test typeof(e(nothing).data) == SparseMatrixCSC{Float32,Int64}
	vs = JSON.parse.(["{\"name\": \"c\", \"count\" : 1}"])
	@test e(vs).data[:] ≈ [0, 0, 1]
	@test typeof(e(vs).data) == SparseMatrixCSC{Float32,Int64}
end

@testset "ExtractCategorical" begin
	e = ExtractCategorical(["a","b"])
	@test e("a").data[:] ≈ [1, 0, 0]
	@test e("b").data[:] ≈ [0, 1, 0]
	@test e("z").data[:] ≈ [0, 0, 1]
	@test e(nothing).data[:] ≈ [0, 0, 1]
	@test typeof(e("a").data) == Flux.OneHotMatrix{Array{Flux.OneHotVector,1}}
	@test typeof(e(nothing).data) == Flux.OneHotMatrix{Array{Flux.OneHotVector,1}}

	@test e(["a", "b"]).data ≈ [1 0; 0 1; 0 0]
	@test typeof(e(["a", "b"]).data) == Flux.OneHotMatrix{Array{Flux.OneHotVector,1}}
end

@testset "show" begin
	e = ExtractCategorical(["a","b"])
	buf = IOBuffer()
	Base.show(buf, e)
	str_repr = String(take!(buf))
	@test str_repr == """Categorical d = 3\n"""

	e = ExtractOneHot(["a","b"], "name", nothing)
	buf = IOBuffer()
	Base.show(buf, e)
	str_repr = String(take!(buf))
	@test str_repr == """OneHot d = 3\n"""

	other = Dict("a" => ExtractArray(ExtractScalar(Float64,2,3)),"b" => ExtractArray(ExtractScalar(Float64,2,3)));
	br = ExtractBranch(nothing,other)
	buf = IOBuffer()
	Base.show(buf, br)
	str_repr = String(take!(buf))
	@test str_repr ==
"""
: struct
  Empty vec
  Other:
  ├── a: Array of
  │    └── Float64
  └── b: Array of
      └── Float64
"""

	vector = Dict("a" => ExtractScalar(Float64,2,3),"b" => ExtractScalar(Float64))
	other = Dict("c" => ExtractArray(ExtractScalar(Float64,2,3)))
	br = ExtractBranch(vector,other)
	buf = IOBuffer()
	Base.show(buf, br)
	str_repr = String(take!(buf))
	@test str_repr ==
"""
: struct
  Vec:
  ├── a: Float64
  └── b: Float64
  Other:
  └── c: Array of
      └── Float64
"""

	other1 = Dict("a" => ExtractArray(ExtractScalar(Float64,2,3)),"b" => ExtractArray(ExtractScalar(Float64,2,3)))
	br1 = ExtractBranch(nothing,other1)
	other = Dict("a" => ExtractArray(br1), "b" => ExtractScalar(Float64,2,3))
	br = ExtractBranch(nothing,other)
	buf = IOBuffer()
	Base.show(buf, br)
	str_repr = String(take!(buf))
	@test str_repr ==
"""
: struct
  Empty vec
  Other:
  ├── a: Array of
  │    └── : struct
  │          Empty vec
  │          Other:
  │          ├── a: Array of
  │          │    └── Float64
  │          └── b: Array of
  │              └── Float64
  └── b: Float64
"""
end

@testset "equals and hash test" begin
	other1 = Dict(
		"a" => ExtractArray(ExtractScalar(Float64,2,3)),
		"b" => ExtractArray(ExtractScalar(Float64,2,3)),
		"c" => ExtractCategorical(["a","b"]),
	)
	br1 = ExtractBranch(nothing,other1)
	other11 = Dict(
		"a" => ExtractArray(br1),
		"b" => ExtractScalar(Float64,2,3),
	)
	br11 = ExtractBranch(nothing,other11)

	other2 = Dict(
		"a" => ExtractArray(ExtractScalar(Float64,2,3)),
		"b" => ExtractArray(ExtractScalar(Float64,2,3)),
		"c" => ExtractCategorical(["a","b"]),
	)
	br2 = ExtractBranch(nothing,other2)
	other22 = Dict("a" => ExtractArray(br2), "b" => ExtractScalar(Float64,2,3))
	br22 = ExtractBranch(nothing,other22)

	@test hash(br11) === hash(br22)
	@test hash(br11) !== hash(br1)
	@test hash(br11) !== hash(br2)
	@test hash(br1) === hash(br2)

	@test br11 == br22
	@test br11 != br1
	@test br11 != br2
	@test br1 == br2
end


text:
"""
  This interface is designed to provide information about the number and types of Phases supported by the component that implements it. It defines all the Phases that a component such as a Physical Property Calculator can handle. It does not provide information about the Phases that are actually present in a Material Object. This function is provided by the GetPresentPhases method of the ICapeThermoMaterial interface.
"""
module ICapeThermoPhases

    export getnumphases!,getphaseinfo!,getphaselist!
    """
      Returns the number of Phases.
    """
    function getnumphases!(
        #= [retval][out] =# num::Int32)
		end

    """
      Returns information on an attribute associated with a Phase for the purpose of understanding what lies behind a Phase label.
    """
    function getphaseinfo!(
        #= [in] =# phaseLabel::Vector{AbstractString},
        #= [in] =# phaseAttribute::Vector{AbstractString},
        #= [retval][out] =# value::Any)
		end

    """
      Returns Phase labels and other important descriptive information for all the Phases supported.
    """
    function getphaselist!(
        #= [out][in] =# phaseLabels!::Vector{AbstractString},
        #= [out][in] =# stateOfAggregation!::Vector{AbstractString},
        #= [out][in] =# keyCompoundId!::Vector{AbstractString})
		end

end


text:
module TestGraphAlgorithms
using Test

using Catlab.Graphs

# Connectivity
##############

g = Graph(6)
add_edges!(g, [2,3], [4,5])
@test connected_components(g) == [[1], [2,4], [3,5], [6]]

# DAGs
######

# Topological sorting
#--------------------

# Discrete graph.
g = Graph(5)
@test topological_sort(g) == 1:5

# Path graph.
g = Graph(5)
add_edges!(g, 5:-1:2, 4:-1:1)
@test topological_sort(g) == 5:-1:1

# Diamond.
g = Graph(4)
add_edges!(g, [1,1,2,3], [2,3,4,4])
@test topological_sort(g) == 1:4

# Multiple edges.
g = Graph(3)
add_edges!(g, [1,1,1,2], [2,2,3,3])
@test topological_sort(g) == 1:3

# Error handling: cyclic graph.
g = Graph(3)
add_edges!(g, [1,2,3], [2,3,1])
@test_throws Exception topological_sort(g)

# Transitivity
#-------------

g = Graph(3)
add_edges!(g, [1,2], [2,3])
g₀ = copy(g)
@test transitive_reduction!(g) == g₀
add_edge!(g, 1, 3)
@test transitive_reduction!(g) == g₀

g = Graph(5)
add_edges!(g, [1,1,1,1,2,3,3,4], [2,3,4,5,4,4,5,5])
transitive_reduction!(g)
@test sort!(collect(zip(src(g), tgt(g)))) == [(1,2),(1,3),(2,4),(3,4),(4,5)]

end


text:
using WebAPI
using HTTP
using JSON3
using Sockets
using Test

# for server tests
function bhaskara(a, b, c)
    Δ = b^2 - 4*a*c

    Δ < 0 && (Δ = complex(Δ))

    x₁ = (-b + √Δ) / 2a
    x₂ = (-b - √Δ) / 2a
    return x₁, x₂
end

@testset "WebAPI.jl" begin
    include("params_query.jl")
    include("req_body.jl")
    include("res_constructors.jl")
    include("add_route.jl")
    include("utility_functions.jl")
    include("prints.jl")
    include("server.jl")
end