text:
# Taken and adapted from LRUCache.jl
mutable struct LRUNode{K, V}
    k::K
    v::V
    next::LRUNode{K, V}
    prev::LRUNode{K, V}

    # All new created nodes are self referential only
    function LRUNode{K, V}(k::K, v::V) where {K, V}
        x = new{K, V}(k, v)
        x.next = x
        x.prev = x
        return x
    end
end

mutable struct LRUList{K, V}
    first::Union{LRUNode{K, V}, Nothing}
    size::Int64

    LRUList{K, V}() where {K, V} = new{K, V}(nothing, 0)
end

Base.first(l::LRUList) = !isempty(l) ? l.first : error("LRUList is empty")
Base.last(l::LRUList) = !isempty(l) ? l.first.prev : error("LRUList is empty")

Base.length(l::LRUList) = l.size
Base.isempty(l::LRUList) = length(l) == 0

function Base.show(io::IO, l::LRUNode{K, V}) where {K, V}
    print(io, "LRUNode{", K, ", ", V, "}(")
    show(io, l.k)
    print(io, ", ")
    show(io, l.v)
    print(io, ")")
end

function Base.show(io::IO, l::LRUList{K, V}) where {K, V}
    print(io, "LRUList{", K, ", ", V, "}(")
    if length(l) != 0
        f = first(l)
        show(io, f.k)
        print(io, "=>")
        show(io, f.v)
        n = f.next
        while n !== f
            print(io, ", ")
            show(io, n.k)
            print(io, "=>")
            show(io, n.v)
            n = n.next
        end
    end
    print(io, ")")
end

function Base.push!(list::LRUNode{K, V}, new::LRUNode{K, V}) where {K, V}
    new.next = list
    new.prev = list.prev
    list.prev.next = new
    list.prev = new
    return list
end

function Base.push!(l::LRUList{K, V}, el::LRUNode{K, V}) where {K, V}
    if isempty(l)
        l.first = el
    else
        push!(l.first, el)
    end
    l.size += 1
    return l
end

function Base.pop!(l::LRUList{K, V}, n::LRUNode{K, V}=last(l)) where {K, V}
    if n.next === n
        l.first = nothing
    else
        if n === first(l)
            l.first = n.next
        end
        n.next.prev = n.prev
        n.prev.next = n.next
    end
    l.size -= 1
    return n
end

function Base.pushfirst!(l::LRUList{K, V}, el::LRUNode{K, V}) where {K, V}
    push!(l, el)
    rotate!(l)
end

# Rotate one step forward, so last element is now first
function rotate!(l::LRUList)
    if length(l) > 1
        l.first = first(l).prev
    end
    return l
end

# Move the node n to the front of the list
function move_to_front!(l::LRUList{T}, n::LRUNode{T}) where {T}
    if first(l) !== n
        pop!(l, n)
        pushfirst!(l, n)
    end
    return l
end

function Base.delete!(l::LRUList{K, V}, n::LRUNode{K, V}) where {K, V}
    pop!(l, n)
    return l
end

function Base.empty!(l::LRUList{K, V}) where {K, V}
    l.first = nothing
    l.size = 0
    return l
end

# Default cache size
const __maxfraction__ = 0.5

mutable struct LRU{K,V} <: AbstractDict{K,V}
    ht::Dict{K, LRUNode{K, V}}
    hts::Dict{K, Int64}
    q::LRUList{K, V}
    currentsize::Int64
    maxsize::Int64

    LRU{K, V}(m::Integer) where {K, V} = new{K, V}(Dict{K, V}(), Dict{K, Int64}(), LRUList{K, V}(), 0, m)
end
function LRU{K, V}(; maxsize::Integer = 0, maxrelsize::Real = 0) where {K, V}
    if maxrelsize == 0 && maxsize == 0
        m = floor(Int64, __maxfraction__ * Sys.total_memory())
    else
        m = max(maxsize, floor(Int64, maxrelsize*Sys.total_memory()))
    end
    LRU{K,V}(m)
end
LRU(; kwargs...) = LRU{Any, Any}(; kwargs...)

Base.show(io::IO, lru::LRU{K, V}) where {K, V} = print(io,"LRU{$K, $V}($(lru.maxsize))")

function Base.iterate(lru::LRU, state...)
    next = iterate(lru.ht, state...)
    if next === nothing
        return nothing
    else
        (k, node), state = next
        return k=>node.v, state
    end
end

Base.length(lru::LRU) = length(lru.q)
Base.isempty(lru::LRU) = isempty(lru.q)
function Base.sizehint!(lru::LRU, n::Integer)
    sizehint!(lru.ht, n)
    sizehint!(lru.hts, n)
    return lru
end

Base.haskey(lru::LRU, key) = haskey(lru.ht, key)
Base.get(lru::LRU, key, default) = haskey(lru, key) ? lru[key] : default

function Base.get!(default::Base.Callable, lru::LRU{K, V}, key::K) where {K,V}
    if haskey(lru, key)
        return lru[key]
    else
        value = default()
        lru[key] = value
        return value
    end
end

function Base.get!(lru::LRU{K,V}, key::K, default::V) where {K,V}
    if haskey(lru, key)
        return lru[key]
    else
        lru[key] = default
        return default
    end
end

function Base.getindex(lru::LRU, key)
    node = lru.ht[key]
    move_to_front!(lru.q, node)
    return node.v
end

function Base.setindex!(lru::LRU{K, V}, v, key) where {K, V}
    if haskey(lru, key)
        item = lru.ht[key]
        lru.currentsize -= lru.hts[key]
        item.v = v
        s = Base.summarysize(v)
        lru.currentsize += s
        lru.hts[key] = s
        move_to_front!(lru.q, item)
    else
        item = LRUNode{K, V}(key, v)
        pushfirst!(lru.q, item)
        lru.ht[key] = item
        s = Base.summarysize(v)
        lru.currentsize += s
        lru.hts[key] = s
    end
    while lru.currentsize > lru.maxsize
        rm = pop!(lru.q)
        lru.currentsize -= lru.hts[rm.k]
        delete!(lru.ht, rm.k)
        delete!(lru.hts, rm.k)
    end
    return lru
end

function setsize!(lru::LRU; maxsize::Integer = 0, maxrelsize::Real = 0)
    @assert 0 <= maxsize
    @assert 0 <= maxrelsize < 1
    if maxrelsize == 0 && maxsize == 0
        m = floor(Int64, __maxfraction__ * Sys.total_memory())
    else
        m = max(maxsize, floor(Int64, maxrelsize*Sys.total_memory()))
    end
    lru.maxsize = m
    while lru.currentsize > lru.maxsize
        rm = pop!(lru.q)
        lru.currentsize -= lru.hts[rm.k]
        delete!(lru.ht, rm.k)
        delete!(lru.hts, rm.k)
    end
end

function Base.delete!(lru::LRU, key)
    item = lru.ht[key]
    lru.currentsize -= lru.hts[key]
    delete!(lru.q, item)
    delete!(lru.ht, key)
    delete!(lru.hts, key)
    return lru
end

function Base.empty!(lru::LRU)
    lru.currentsize = 0
    empty!(lru.q)
    empty!(lru.ht)
    empty!(lru.hts)
end


text:
@testset "precondition" begin
    A = [2..4 -2..1; -1..2 2..4]
    b = [-2..2, -2..2]

    np = NoPrecondition()
    idp = InverseDiagonalMidpoint()
    imp = InverseMidpoint()

    A1, b1 = np(A, b)
    @test A1 == A && b1 == b

    A2, b2 = idp(A, b)
    Acorrect = [2/3..4/3 -2/3..1/3; -1/3..2/3 2/3..4/3]
    bcorrect = [-2/3..2/3, -2/3..2/3]
    @test all(interval_isapprox.(A2, Acorrect)) && all(interval_isapprox.(b2, bcorrect))

    A3, b3 = imp(A, b)
    Acorrect = [22/37..52/37 -20/37..20/37;-20/37..20/37 22/37..52/37]
    bcorrect = [-28/37..28/37, -28/37..28/37]
    @test all(interval_isapprox.(A3, Acorrect)) && all(interval_isapprox.(b3, bcorrect))
end

ext: jl


text:
include("downloads.jl")
using GLFW
tutorials = readdir(dirname(@__FILE__))
filter!(tutorials) do file
    endswith(file, ".jl") &&
    basename(file) != "downloads.jl" &&
    basename(file) != "run_all.jl"
end
for tutorial in tutorials
    GLFW.Init()
    println("Now showing: ", tutorial)
    include(tutorial)
    GLFW.Terminate() # clean up gl context
end


text:
# Browser Window (Borrowed from Blink.jl)
@static if is_apple()
    launch(x) = run(`open $x`)
elseif is_linux()
    launch(x) = run(`xdg-open $x`)
elseif is_windows()
    launch(x) = run(`cmd /C start $x`)
end

# const connections = Dict{Int,WebSocket}() # WebSocket.id => WebSocket
const conditions = Dict{String,Condition}()
conditions["connected"] = Condition()
conditions["unloaded"] = Condition()

Endpoint("/pages.js") do request::HTTP.Request
    readstring(joinpath(dirname(@__FILE__),"pages.js"))
end

# ws = WebSocketHandler() do request::Request, client::WebSocket
#     while true
#         msg = JSON.parse(String(read(client)))
#         route = msg["route"]
#         if !haskey(pages[route].sessions,client.id)
#             pages[route].sessions[client.id] = client
#         end
#         haskey(msg,"args") ? callbacks[msg["name"]].callback(client,msg["args"]) : callbacks[msg["name"]].callback(client)
#     end
# end

# For external websocket connections not from a page served locally, e.g. IJulia
const external = Dict{String,WebSocket}()

function start(p = 8000)
    global port = p
    HTTP.listen(ip"127.0.0.1",p) do http
        if HTTP.WebSockets.is_upgrade(http.message)
            HTTP.WebSockets.upgrade(http) do client
                while !eof(client);
                    data = String(readavailable(client))
                    msg = JSON.parse(data)
                    name = pop!(msg,"name"); route = pop!(msg,"route"); id = pop!(msg,"id")
                    if haskey(pages,route)
                        if !haskey(pages[route].sessions,id)
                            pages[route].sessions[id] = client
                        end
                    else
                        if !haskey(external,id)
                            external[id] = client
                        end
                    end
                    if haskey(callbacks,name)
                        callbacks[name].callback(client,route,id,msg)
                    end
                end
            end
        else
            route = HTTP.URI(http.message.target).path
            if haskey(pages,route)
                HTTP.Servers.handle_request(pages[route].handler,http)
            else
                HTTP.Servers.handle_request((req) -> HTTP.Response(404),http)
            end
        end
    end
end


text:
module gmi_types

using apf_types
import apf_types.ModelEntity  # == gmi_ent*

export ModelEntity, NullModelEntity, Model, Gmi_set, _Gmi_iter, Gmi_iter

const NullModelEntity = ModelEntity(C_NULL)

"""
  gmi_model*
"""
struct Model
  p::Ptr{Void}
end

"""
  Gmi_set*.  Has an extra field `v` that contains the ModelEntities themselves.

  Users can call `finalize` on this object to free it if needed.
"""
mutable struct Gmi_set
  p::Ptr{Void}
  v::Vector{ModelEntity}
end
#TODO implement AbstractAray interface?

"""
  gmi_iter*
"""
struct _Gmi_iter
  p::Ptr{Void}
end

"""
  A Julia iterator that behaves like `_Gmi_iter` does in C

  There is no need to free these object when done, they are automatically
  freed when iteration is complete.
"""
struct Gmi_iter
  it::_Gmi_iter
  g::Model
end

end  # end module


text:
include("../src/mathFoo.jl")
include("../src/textFoo.jl")

d = summa(1,2,3)

calcQuadRoots(1,5,2)


text:
using AdventOfCodeSolutions
using Test

function input(puzzle::Puzzle{2020, 11, n}) where n
    io = openInput(puzzle)
    return split(read(io, String), '\n', keepempty=false)
end

function parseInput(input)
    return permutedims(hcat(map(collect, input)...))
end

function countTakenSeats(A, i, j)
    (n, m) = size(A)
    count(
        idx -> A[idx[1], idx[2]] == '#',
        filter(
            !=([i,j]),
            [[l,k] for l = (max(1, i - 1) : min(i + 1, n)), k=max(1, j - 1):min(j + 1, m)]
        )
    )
end

function countTakenSeatsInSight(A, i, j)
    (n, m) = size(A)
    directions = filter(!=([0 0]), [[i j] for i=-1:1, j=-1:1])
    count(directions) do (direction)
        p = [i j]

        p += direction
        while (p[1] ∈ 1:n && p[2] ∈ 1:m)
            if (A[p[1], p[2]] == '#')
                return true
            end
            if (A[p[1], p[2]] == 'L')
                return false
            end

            p += direction
        end
        return false
    end
end

"""
If a seat is empty (L) and there are no occupied seats adjacent to it, the seat becomes occupied.
If a seat is occupied (#) and four or more seats adjacent to it are also occupied, the seat becomes empty.
Otherwise, the seat's state does not change.
"""
function iterateInnerNeighbours(A, result, i, j)
    if A[i, j] == '.'
        return
    end

    takenSeats = countTakenSeats(A, i, j)
    if A[i, j] == 'L' && takenSeats == 0
        result[i, j] = '#'
    end
    if A[i, j] == '#' && takenSeats >= 4
        result[i, j] = 'L'
    end
end

"""
If a seat is empty (L) and there are no occupied seats adjacent to it, the seat becomes occupied.
If a seat is occupied (#) and five or more seats in sight are also occupied, the seat becomes empty.
Otherwise, the seat's state does not change.
"""
function iterateInnerSight(A, result, i, j)
    if A[i, j] == '.'
        return
    end

    takenSeats = countTakenSeatsInSight(A, i, j)
    if A[i, j] == 'L' && takenSeats == 0
        result[i, j] = '#'
    end
    if A[i, j] == '#' && takenSeats >= 5
        result[i, j] = 'L'
    end
end

function iterate(A, iterateFn = iterateInnerNeighbours)
    result = copy(A)
    n, m = size(A)

    for i = 1:n, j =1:m
        iterateFn(A, result, i, j)
    end

    return result
end
iterateNeighbours(A) = iterate(A, iterateInnerNeighbours)
iterateSight(A) = iterate(A, iterateInnerSight)

function solveProblem(A, iterateFn)
    converged = false
    iterations = 0
    while converged != true
        iterations += 1
        previousA = copy(A)
        A = iterateFn(previousA)

        if (A == previousA)
            break
        end

        if iterations > 10000
            @show "HOI"
            break
        end
    end
    count(v -> v == '#', A)
end

function solve(::Puzzle{2020, 11, 1}, input)
    A = parseInput(input)
    return solveProblem(A, iterateNeighbours)
end

function solve(::Puzzle{2020, 11, 2}, input)
    A = parseInput(input)
    return solveProblem(A, iterateSight)
end

@testset "Day 11" begin
    @testset "Examples" begin

        @testset "Neighbouring taken seats" begin
            @test iterateNeighbours(['.' '.' '.'; '.' 'L' '.'; '.' '.' '.';]) == ['.' '.' '.'; '.' '#' '.'; '.' '.' '.';]
            @test iterateNeighbours(['.' '.' '.'; '.' '#' '.'; '.' '.' '.';]) == ['.' '.' '.'; '.' '#' '.'; '.' '.' '.';]
            @test iterateNeighbours(['#' '.' '#'; '#' '#' '#'; '#' '.' '#';]) == ['#' '.' '#'; '#' 'L' '#'; '#' '.' '#';]
        end

        @testset "Takenseats in sight" begin
            INPUT = """
.......#.
...#.....
.#.......
.........
..#L....#
....#....
.........
#........
...#.....
"""

            A = split(INPUT, '\n', keepempty=false) |> parseInput
            @test countTakenSeatsInSight(A, 5, 4) == 8

            INPUT = """
.............
.L.L.#.#.#.#.
.............
"""
            A = split(INPUT, '\n', keepempty=false) |> parseInput
            @test countTakenSeatsInSight(A, 2, 2) == 0

            INPUT = """
.##.##.
#.#.#.#
##...##
...L...
##...##
#.#.#.#
.##.##.
"""
            A = split(INPUT, '\n', keepempty=false) |> parseInput
            @test countTakenSeatsInSight(A, 4, 4) == 0

            INPUT = """
L.LL.LL.LL
LLLLLLL.LL
L.L.L..L..
LLLL.LL.LL
L.LL.LL.LL
L.LLLLL.LL
..L.L.....
LLLLLLLLLL
L.LLLLLL.L
L.LLLLL.LL
"""
            A = split(INPUT, '\n', keepempty=false) |> parseInput

            INPUT = """
#.##.##.##
#######.##
#.#.#..#..
####.##.##
#.##.##.##
#.#####.##
..#.#.....
##########
#.######.#
#.#####.##"""
            expectedA = split(INPUT, '\n', keepempty=false) |> parseInput

            @test iterateSight(A) == expectedA

            INPUT="""
#.LL.LL.L#
#LLLLLL.LL
L.L.L..L..
LLLL.LL.LL
L.LL.LL.LL
L.LLLLL.LL
..L.L.....
LLLLLLLLL#
#.LLLLLL.L
#.LLLLL.L#
"""
            expectedA2 = split(INPUT, '\n', keepempty=false) |> parseInput
            @test iterateSight(expectedA) == expectedA2

            INPUT="""
#.L#.##.L#
#L#####.LL
L.#.#..#..
##L#.##.##
#.##.#L.##
#.#####.#L
..#.#.....
LLL####LL#
#.L#####.L
#.L####.L#
"""
            expectedA3 = split(INPUT, '\n', keepempty=false) |> parseInput
            @test iterateSight(expectedA2) == expectedA3
        end

    end
    TEST_INPUT_1 = """
L.LL.LL.LL
LLLLLLL.LL
L.L.L..L..
LLLL.LL.LL
L.LL.LL.LL
L.LLLLL.LL
..L.L.....
LLLLLLLLLL
L.LLLLLL.L
L.LLLLL.LL
"""

    @test(solve(Puzzle(11, 1), split(TEST_INPUT_1, '\n', keepempty=false)) == 37)
    @test(solve(Puzzle(11, 1), input(Puzzle(11, 1))) == 2489)

    @test(solve(Puzzle(11, 2), split(TEST_INPUT_1, '\n', keepempty=false)) == 26)
    @test(solve(Puzzle(11, 2), input(Puzzle(11, 2))) == 2180)
end


text:
module HeckeProfile
using Hecke

function Profile.callersf(matchfunc::Function, bt::Vector, lidict::Profile.LineInfoFlatDict)
  counts = Dict{Profile.StackFrame, Int}()
  lastmatched = false
  for id in bt
      if id == 0
	  lastmatched = false
	  continue
      end
      li = lidict[id]
      if lastmatched
	  if haskey(counts, li)
	      counts[li] += 1
	  else
	      counts[li] = 1
	  end
      end
      lastmatched = matchfunc(li)
  end
  k = collect(keys(counts))
  v = collect(values(counts))
  p = sortperm(v, rev=true)
  return [(v[i], k[i]) for i in p]
end

function Profile.callers(funcname::Symbol, bt::Vector, lidict::Profile.LineInfoDict; filename = nothing, linerange = nothing)
   bt, lid = Profile.flatten(bt, lidict)

    if filename === nothing && linerange === nothing
        return Profile.callersf(li -> li.func == funcname, bt, lid)
    end
    filename === nothing && throw(ArgumentError("if supplying linerange, you must also supply the filename"))
    if linerange === nothing
        return Profile.callersf(li -> li.func == funcname && li.file == filename, bt, lid)
    else
        return Profile.callersf(li -> li.func == funcname && li.file == filename &&
            in(li.line, linerange), bt, lid)
    end
end

Profile.callers(funcname::Symbol; kwargs...) = Profile.callers(funcname, Profile.retrieve()...; kwargs...)
Profile.callers(func::Function, bt::Vector, lidict::Profile.LineInfoDict; kwargs...) =
    Profile.callers(Symbol(func), bt, lidict; kwargs...)
Profile.callers(func::Function; kwargs...) = Profile.callers(Symbol(func), Profile.retrieve()...; kwargs...)

function Base.process_backtrace(process_func::Function, t::Vector, limit::Int=typemax(Int); skipC = !true)
    n = 0
    last_frame = StackTraces.UNKNOWN
    count = 0
    for i = eachindex(t)
        lkups = StackTraces.lookup(t[i])
        for lkup in lkups
            if lkup === StackTraces.UNKNOWN
                continue
            end

            if lkup.from_c && skipC; continue; end
            if i == 1 && lkup.func == :error; continue; end
            count += 1
            if count > limit; break; end

            if lkup.file != last_frame.file || lkup.line != last_frame.line || lkup.func != last_frame.func
                if n > 0
                    process_func(last_frame, n)
                end
                n = 1
                last_frame = lkup
            else
                n += 1
            end
        end
    end
    if n > 0
        process_func(last_frame, n)
    end
end

#=
  from base/profile.jl
  bt is a vector if UInt
    each complete backtrace corresponds to a block in bt between to 0
    iz = find(bt .== 0)
    then bt[1:iz[1]-1] is a complete trace
         bt[iz[1]+1:iz[2]-1] the next, ...
=#

function callers_frame(fname::Symbol, bt::Vector, lidict::Profile.LineInfoDict)
  counts = []
  st = typeof(bt)()
  skip = false
  found = false
  for id in bt
    skip = skip && id != 0
    skip && continue

    if id == 0
      st = typeof(bt)()
      continue
    end
    li = lidict[id]
    push!(st, id)
    for x in li
#      if contains(String(x.func), "eval_user_input")
      if x.func == :eval_user_input
        if found
          push!(st, 0)
          push!(counts, st)
          found = false
        end
        skip = true
        break
      end
      if x.func == fname
        found = true
        break
      end
    end
  end
  return counts
end

function leafs(bt::Vector, lidict::Profile.LineInfoDict, skipC::Bool = false)
  M = Hecke.MSet{Symbol}()
  skip = false
  for f = bt
    if !skip
      if f == 0
        continue
      end
      for fr = lidict[f]
        if skipC && fr.from_c
          continue
        end
        push!(M, fr.func)
        skip = true
        break
      end
    else
      skip = f!=0
    end
  end
  return M
end

function counts(bt::Vector, lidict::Profile.LineInfoDict, skipC::Bool = false)
  M = Hecke.MSet{Symbol}()
  for f = bt
    if f == 0
      continue
    end
    for fr = lidict[f]
      if skipC
        if !fr.from_c
          push!(M, fr.func)
        end
      else
        push!(M, fr.func)
      end
    end
  end
  return M
end

#all frames starting with fname
function prune(fname::Symbol, bt::Vector, lidict::Profile.LineInfoDict)
  rt = typeof(bt)()
  last = typeof(bt)()
  bingo = false
  for i=length(bt):-1:1
    g = bt[i]
    if g==0
      if bingo
        append!(rt, reverse(last))
        push!(rt, g)
        last = typeof(bt)()
      end
      @assert length(last) == 0
      bingo = false
      continue
    end
    for f = g
      if lidict[f][1].func == fname
        bingo = true
      end
      if bingo
        push!(last, g)
        break
      end
    end
  end
  if bingo
    append!(rt, reverse(last))
    push!(rt, UInt(0))
  end
  return rt
end

mutable struct Graph{T <: Any}
  v::Dict{T, Int}
  e::Dict{Tuple{T, T}, Int}
  function Graph{T}() where T
    r = new()
    r.v = Dict{T, Int}()
    r.e = Dict{Tuple{T, T}, Int}()
    return r
  end
end

function symb(a::StackFrame)
  if a.func==:Type
    if isdefined(a, :linfo) && isdefined(a.linfo, :value) && isdefined(a.linfo.value, :specTypes)
      return Symbol(a.linfo.value.specTypes)
    else
      return a.func
    end
  else
    return a.func
  end
end

function graph(bt::Vector, lidict::Profile.LineInfoDict, skipC::Bool = true, skipMacro::Bool=true)
  g = Graph{Symbol}()
  v = g.v
  e = g.e

  last = :rt
  Mac = Symbol("macro expansion")
  for f = bt
    if f==0
      last = :rt
      continue
    end
    for y = lidict[f]
      if skipC && y.from_c
        continue
      end
      fn = symb(y)
      if skipMacro && fn == Mac
        continue
      end
      if haskey(v, fn)
        v[fn] += 1
      else
        v[fn] = 1
      end
      t = (fn, last)
      if haskey(e, t)
        e[t] += 1
      else
        e[t] = 1
      end
      last = fn
    end
  end
  return g
end

graph(skipC::Bool = true, skipMacro::Bool=true) = graph(Profile.retrieve()..., skipC, skipMacro)

function parents{T <:Any}(g::Graph{T}, c::T)
  return [a for (a,b) = keys(g.e) if b==c]
end

function parents_with_counts{T <:Any}(g::Graph{T}, c::T)
  return [(a, d) for ((a,b), d) = g.e if b==c]
end

function children{T <: Any}(g::Graph{T}, c::T)
  return [b for (a,b) = keys(g.e) if a==c]
end

function children_with_count{T <: Any}(g::Graph{T}, c::T)
  return [(b, d) for ((a,b), d) = g.e if a==c]
end

end


text:
"""
Simple Tucker (HOSVD) type
# Data
- `core::Array{T, N}`
- `factors::NTuple{N, Matrix{T}}`

Tucker factors are stored as tall matrices
"""
struct Tucker{T, N}
    core::Array{T, N}
    factors::NTuple{N, Matrix{T}}
    #props::Dict{Symbol, Any}
end

function Tucker_tot(A::Array{T,N}; thresh=-1, verbose=0) where {T,N}
    core,factors = tucker_decompose_tot(A, thresh=thresh, verbose=verbose)
    return Tucker{T,N}(core, NTuple{N}(factors))
end
function Tucker(A::Array{T,N}; thresh=-1, max_number=nothing, verbose=0, type="magnitude") where {T,N}
    core,factors = tucker_decompose(A, thresh=thresh, max_number=max_number, verbose=verbose, type=type)
    return Tucker{T,N}(core, NTuple{N}(factors))
end
recompose(t::Tucker{T,N}) where {T<:Number, N} = tucker_recompose(t.core, t.factors)
dims_large(t::Tucker{T,N}) where {T<:Number, N} = return [size(f,1) for f in t.factors]
dims_small(t::Tucker{T,N}) where {T<:Number, N} = return [size(f,2) for f in t.factors]
Base.length(t::Tucker) = prod(dims_small(t))
Base.size(t::Tucker) = size(t.core)
function Base.permutedims(t::Tucker{T,N}, perm) where {T,N}
    #t.core .= permutedims(t.core, perm)
    return Tucker{T,N}(permutedims(t.core, perm), t.factors[perm])
end

"""
    add(tucks::Vector{Tucker{T,N}}; thresh=1e-10, max_number=nothing) where {T,N}

Add together multiple Tucker instances. Assumed non-orthogonal.

# Arguments
- `tucks::Vector{Tucker{T,N}}`: Vector of Tucker objects
"""
function add(tucks::Vector{Tucker{T,N}}; thresh=1e-10, max_number=nothing, type="magnitude") where {T,N}

    # sort the Tucker objects to add. This puts them in a well-defined order for reproducibility.
    norms = norm.(tucks)
    perm = sortperm(norms,rev=true)
    tucks = tucks[perm]
    #display(norm.(tucks))
    length(tucks) > 0 ||  error("not enough Tuckers to add", length(tucks))
    length(tucks) > 1 ||  return tucks[1]

    # first, ensure the Tuckers all correpsond to appropriately dimensioned uncompressed tensors
    dims = dims_large(first(tucks))
    for tuck in tucks
        all(dims_large(tuck) .== dims) || error("dimension error.", dims_large.(tucks))
    end

    # first, get orthogonal basis spanning all Tucker objects for each index
    new_factors = Vector{Matrix{T}}()
    for i in 1:N
        Ui = Vector{Matrix{T}}()
        for tuck in tucks
            push!(Ui, tuck.factors[i])
        end
        Ui = hcat(Ui...)

        F = svd(Ui)

        nkeep = 0
        if type == "magnitude"
            for si in F.S
                if si > thresh
                    nkeep += 1
                end
            end
        elseif type == "sum"
            target = sum(F.S )
            curr = 0.0
            for si in F.S
                if abs(curr-target) > thresh
                    nkeep += 1
                    curr += si*si
                end
            end
        else
            error("wrong type")
        end
        if max_number != nothing
            nkeep = min(nkeep, max_number)
        end
        Ui = F.U[:,1:nkeep]
        push!(new_factors, Ui)
    end

    new_core = zeros([size(new_factors[i],2) for i in 1:N]...)

    for tuck in tucks
        transforms = Vector{Matrix{Float64}}()
        for i in 1:N
            push!(transforms, new_factors[i]'*tuck.factors[i])
        end

        new_core += recompose(Tucker{T,N}(tuck.core, Tuple(transforms)))
    end
    return Tucker{T,N}(new_core, Tuple(new_factors))

end
Base.:(+)(t1::Tucker, t2::Tucker) = add([t1, t2])
nonorth_add(tl::Vector{Tucker{T,N}}; thresh=1e-7) where {T,N}= add(tl, thresh=thresh)
nonorth_add(t1::Tucker, t2::Tucker; thresh=1e-7) = add([t1,t2], thresh=thresh)
nonorth_dot(t1, t2) = dot(t1,t2)
function scale(t1::Tucker, a)
    t2core = t1.core .* a
    return Tucker(t2core, t1.factors)
end
function scale!(t1::Tucker, a)
    t1.core .*= a
    return
end

"""
    compress(t::Tucker{T,N}; thresh=1e-7, max_number=nothing) where {T,N}

Try to compress further
"""
function compress(t::Tucker{T,N}; thresh=1e-7, max_number=nothing, type="magnitude") where {T,N}

    length(t) > 0 || return t
    tt = Tucker(t.core, thresh=thresh, max_number=max_number, type=type)

    new_factors = [zeros(1,1) for i in 1:N]

    for i in 1:N
        new_factors[i]  = t.factors[i] * tt.factors[i]
    end

    return Tucker(tt.core, NTuple{N}(new_factors))
#    return Tucker(recompose(t), thresh=thresh, max_number=max_number)
end

"""
    dot(t1::Tucker{T,N}, t2::Tucker{T,N}) where {T,N}

Note: This doesn't assume `t1` and `t2` have the same compression vectors
"""
function dot(t1::Tucker{T,N}, t2::Tucker{T,N}) where {T,N}
    #overlaps = []
    overlaps = Dict{Int,Matrix{T}}()
    all(dims_large(t1) .== dims_large(t2)) || error(" t1 and t2 don't have same dimensions")
    for f in 1:N
        #push!(overlaps, t1.factors[f]' * t2.factors[f])
        overlaps[f] = t1.factors[f]' * t2.factors[f]
    end
    return sum(transform_basis(t1.core, overlaps) .* t2.core)
    #return sum(tucker_recompose(t1.core, overlaps) .* t2.core)
end

"""
    function tucker_decompose(A::Array{T,N}; thresh=1e-7, max_number=nothing, verbose=1, type="magnitude") where {T,N}

Tucker Decomposition of dense tensor:
    A ~ X *(1) U1 *(2) U2 ....
where cluster states are discarded based on the corresponding SVD
#Arguments
- `A`: matrix to decompose
- `thresh`: threshold for discarding tucker factors
- `max_number`: limit number of tucker factors to this value
- `type`: type of trunctation. "magnitude" discards values smaller than this number.
    "sum" discards values such that the sum of discarded values is smaller than `thresh`.
"""
function tucker_decompose(A::Array{T,N}; thresh=1e-7, max_number=nothing, verbose=1, type="magnitude") where {T,N}
    factors = Vector{Matrix{T}}()
    if verbose > 0
        println(" Tucker Decompose:", size(A))
    end
    length(A) > 0 || error(" can't decompose array with zero data")
    for i in 1:ndims(A)
        idx = collect(1:ndims(A))
        idx[i] = -1
        perm = sortperm(idx)
        U,Σ, = svd(reshape(permutedims(A,perm), size(A,i), length(A)÷size(A,i)))

#        idx_l = collect(1:ndims(A))
#        idx_r = collect(1:ndims(A))
#        idx_l[i] = -1
#        idx_r[i] = -2
#        G = tensorcontract(A,idx_l,A,idx_r)
#        #G = @ncon([A, A], [idx_l, idx_r])
#        F = eigen((G .+ G') .* .5) # should be symmetric, but sometimes values get very small and numerical error builds up
#        perm = sortperm(real(F.values), rev=true)
#        l = F.values[perm]
#        v = F.vectors[:,perm]

        nkeep = 0
        if verbose > 0
            @printf(" index dimension: %6i\n", size(A)[i])
        end
        nkeep = 0
        if type == "magnitude"
            for (idx,Σi) in enumerate(Σ)
                if abs(Σi) > thresh
                    nkeep += 1
                    if verbose > 0
                        @printf("   Singular Value %4i = %12.8f\n", idx, Σi)
                    end
                end
            end
        elseif type == "sum"
            target = sum(l)
            curr = 0.0
            for (idx,Σi) in enumerate(Σ)
                if abs(curr-target) > thresh
                    nkeep += 1
                    curr += Σi
                    if verbose > 0
                        @printf("   Singular Value %4i = %12.8f\n", idx, Σi)
                    end
                end
            end
        else
            error("wrong type")
        end
        if max_number != nothing
            nkeep = min(nkeep, max_number)
        end

        push!(factors, U[:,1:nkeep])
    end
    return transform_basis(A,factors), factors
end

"""
    function tucker_decompose_tot(A::Array{T,N}; thresh=1e-7, verbose=1) where {T,N}

Tucker Decomposition of dense tensor:
    A ~ X *(1) U1 *(2) U2 ....
where the cluster states are discarded to ensure ||V-v||_F < thresh, with V and v being the full and
approximated tensors. This isn't quite ready for use, as it doens't fall back to the optimal case for
SVD of a matrix as it doesn't consider the possibility of being diagonal in svd basis.
"""
function tucker_decompose_tot(A::Array{T,N}; thresh=1e-7, verbose=1) where {T,N}
    factors = Vector{Matrix{T}}()
    values = []
    inds = []
    if verbose > 0
        println(" Tucker Decompose:", size(A))
    end
    for i in 1:ndims(A)
        idx_l = collect(1:ndims(A))
        idx_r = collect(1:ndims(A))
        idx_l[i] = -1
        idx_r[i] = -2
        G = tensorcontract(A,idx_l,A,idx_r)
        #G = @ncon([A, A], [idx_l, idx_r])
        F = eigen((G .+ G') .* .5) # should be symmetric, but sometimes values get very small and numerical error builds up
        perm = sortperm(real(F.values), rev=true)
        l = F.values[perm]
        v = F.vectors[:,perm]

        println(i)
        for li in l
            if verbose > 0
                @printf(" Eigenvalue = %12.8f\n", li)
            end
        end
        #if max_number != nothing
        #    nkeep = min(nkeep, max_number)
        #end

        push!(factors, v)
        append!(values, sqrt.(abs.(l)))
        append!(inds, [i for j in 1:length(l)])
    end

    perm = sortperm(real(values), rev=true)
    values = values[perm]
    inds = inds[perm]

    dims = zeros(Int,N)
    target = sum(values)
    curr = 0.0
    nkeep = 0
    for (idx,i) in enumerate(values)
        if abs(curr-target) > thresh
            nkeep += 1
            curr += i
        end
    end
    for i in 1:nkeep
        dims[inds[i]] += 1
    end
    for i in 1:N
        factors[i] = factors[i][:,1:dims[i]]
    end
    display(target)
    display(size.(factors,2))
    return transform_basis(A,factors), factors
end

"""
    tucker_recompose(core, factors)

Recompose Tucker Decomposition
"""
tucker_recompose(core, factors) = transform_basis(core, factors, trans=true)

"""
"""
function transform_basis(v::Array{T,N}, transform_list::Dict{Int,Matrix{T}}; trans=false) where {T,N}
  #={{{=#
    #
    #   e.g.,
    #   v(i,j,k,l) U(iI) U(jJ) U(lL) = V(I,J,k,L)
    #
    #   vv(i,jkl) = reshape(vv(i,j,k,l))
    #
    #   vv(jkl,I) = v(i,jkl)' * U(i,I)
    #   vv(j,klI) = reshape(vv(jkl,I))
    #
    #   vv(klI,J) = v(j,klI)' * U(j,J)
    #   vv(k,lIJ) = reshape(vv(klI,J))
    #
    #   vv(lIJ,k) = v(k,lIJ)'
    #   vv(l,IJK) = reshape(vv(lIJ,k))
    #
    #   vv(IJk,L) = v(l,IJk)' * U(l,L)
    #   vv(I,JKL) = reshape(vv(IJK,L))

    length(transform_list) > 0 || return v
    #display(("v:",size(v)))
    #display(("t:",[(i,size(j)) for (i,j) in transform_list]))
    vv = deepcopy(v)
    dims = [size(vv)...]

    vv = reshape(vv,dims[1],prod(dims[2:end]))

    for i in 1:N
        if haskey(transform_list, i)

            if trans
                vv = vv' * transform_list[i]'
                dims[1] = size(transform_list[i])[1]
            else
                vv = vv' * transform_list[i]
                dims[1] = size(transform_list[i])[2]
            end

            dims = circshift(dims, -1)
            vv = reshape(vv,dims[1],prod(dims[2:end]))
        else
            vv = vv'
            dims = circshift(dims, -1)
            vv = reshape(vv,dims[1],prod(dims[2:end]))
        end
    end

    return reshape(vv,dims...)
end
#=}}}=#

"""
"""
function transform_basis(v::Array{T,N}, transforms::NTuple{N,Matrix{T}}; trans=false) where {T,N}
  #={{{=#
    vv = deepcopy(v)
    dims = [size(vv)...]

    vv = reshape(vv,dims[1],prod(dims[2:end]))

    for i in 1:N

        if trans
            vv = vv' * transforms[i]'
            dims[1] = size(transforms[i])[1]
        else
            vv = vv' * transforms[i]
            dims[1] = size(transforms[i])[2]
        end

        dims = circshift(dims, -1)
        vv = reshape(vv,dims[1],prod(dims[2:end]))
    end

    return reshape(vv,dims...)
end
#=}}}=#

transform_basis(v::Array{T,N}, transforms; trans=false) where {T,N} = transform_basis(v, NTuple{N, Matrix{T}}(transforms), trans=trans)

function unfold(A::AbstractArray{T,N}, i::Integer) where {T,N}
    # https://github.com/Jutho/TensorOperations.jl/issues/13
    B = permutedims(A,vcat(i,setdiff(1:N,i)))
    d = size(B,1)
    return reshape(B,(d,div(length(B),d)))
end

LinearAlgebra.norm(A::Tucker{T,N}) where {T,N} = norm(A.core)

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
using Printf

mutable struct TimeText; end

function draw_item(timetext::TimeText, canvas::Gtk.GtkCanvas, scene::Scene)
    ctx = Gtk.getgc(canvas)
    h = height(canvas)
    w = width(canvas)

    set_source_rgb(ctx, 1, 1, 1)
    rectangle(ctx, 0, 0, w, h)
    fill(ctx)

    set_source_rgb(ctx, 0.5, 0.0, 0.7)
    timesincelastframe = time()-scene.lastTimeCalled
    txt = @sprintf("%07.3fs since last frame", timesincelastframe)
    text(ctx,0,18,txt)
    scene.lastTimeCalled = time()
    return nothing
end


text:
# Script to obtain sectoral-country trade elasticities

# Using trade elasticities from Imbs and Mejean (2017) (http://www.isabellemejean.com/data.html)
# http://www.isabellemejean.com/data.html in the folder for the paper "Trade elasticities", "/data_trade_elasticities/sigma_v9195.txt"

using DataFrames, CSV, XLSX, Statistics

# ---------- Match elasticities with WIOD ISO3 names -------------------------------------------------------------------------------------------------------

# elasticities
raw_sigma = CSV.read("data/sigma_v9195.txt", DataFrame) # data for elasticities
rename!(raw_sigma, :j => :m49) # change column name to join on identifier "m49"

# M49 vis-a-vis ISO3 country names
raw_names = CSV.read("data/country_m49_iso3.csv", DataFrame)

# join tables
c_table_ctry = leftjoin(raw_sigma, raw_names, on = :m49) # join table
c_table_ctry = c_table_ctry[completecases(c_table_ctry), :] # delete rows with missing values
transform!(c_table_ctry, [:isic3, :model, :method, :country, :iso3] .=> ByRow(string) .=> [:isic3, :model, :method, :country, :iso3], renamecols=false) # change type to string to use XLSX.writetable()

XLSX.writetable("clean/table_elasticities.xlsx", c_table_ctry, overwrite=true)

# ---------- Correspondence table for ISIC to NACE -------------------------------------------------------------------------------------------------------

# Cannot match directly between ISIC Rev.3 and NACE Rev.2
# Solution: ISIC Rev.3 to ISIC Rev.3.1 to ISIC Rev.4 to NACE Rev.2 to classification used in paper (i.e. small aggregations)

# ISIC Rev.3 vis-a-vis ISIC Rev.3.1
raw_ISIC3_ISIC31 = CSV.read("data/ISIC3_ISIC31.txt", DataFrame)[:, Not(["partial3", "partial31", "Activity"])]
rename!(raw_ISIC3_ISIC31, [:ISIC3, :ISIC31])
raw_ISIC3_ISIC31 = string.(raw_ISIC3_ISIC31)
raw_ISIC3_ISIC31 = lpad.(raw_ISIC3_ISIC31, 4, '0') # pad with leading zeros so all have 4 digits

# ISIC Rev.3.1 vis-a-vis ISIC Rev.4
raw_ISIC31_ISIC4 = CSV.read("data/ISIC31_ISIC4.txt", DataFrame)[:, Not(["partialISIC31", "partialISIC4", "Detail"])]
rename!(raw_ISIC31_ISIC4, [:ISIC31, :ISIC4])
raw_ISIC31_ISIC4 = string.(raw_ISIC31_ISIC4)
raw_ISIC31_ISIC4 = lpad.(raw_ISIC31_ISIC4, 4, '0') # pad with leading zeros so all have 4 digits

# ISIC Rev.4 vis-a-vis NACE Rev.2
raw_ISIC4_NACE2 = CSV.read("data/ISIC4_NACE2.csv", DataFrame)[2:end, 1:2]
rename!(raw_ISIC4_NACE2, [:ISIC4, :CPA2008]) # rename to CPA2008 as only numerical NACE codes (i.e. no A01 but just 01 which correspond to CPA2008)
raw_ISIC4_NACE2 = string.(raw_ISIC4_NACE2) # remove all one digit codes (otherwise bounds error below since we select first two characters)
filter!(row -> length(row.CPA2008) > 1, raw_ISIC4_NACE2)
f(x) = x[1:2] # select first two elements of cell
raw_ISIC4_NACE2[:, :CPA2008] = f.(raw_ISIC4_NACE2[:, :CPA2008]) # aggregates CPA2008 to two digit level (i.e. match all codes to just two digit level)

# CPA2008 (2 digits) vis-a-vis NACE Rev.2 (and the WIOD classification of sectors which aggregates some NACE)
raw_CPA_NACE = CSV.read("data/CPA2008_NACE2.csv", DataFrame)
transform!(raw_CPA_NACE, names(raw_CPA_NACE) .=> ByRow(string) .=> names(raw_CPA_NACE), renamecols=false) # change type to string
raw_CPA_NACE[:, :CPA2digits] = lpad.(raw_CPA_NACE[:, :CPA2digits], 2, '0') # pad with leading zeros so all have 2 digits

# join table to create a single correspondence table
c_table_ISIC_NACE = leftjoin(raw_ISIC3_ISIC31, raw_ISIC31_ISIC4, on=:ISIC31)
c_table_ISIC_NACE = leftjoin(c_table_ISIC_NACE, raw_ISIC4_NACE2, on=:ISIC4)
rename!(c_table_ISIC_NACE, :CPA2008 => :CPA2digits)
c_table_ISIC_NACE = leftjoin(c_table_ISIC_NACE, raw_CPA_NACE, on=:CPA2digits)

XLSX.writetable("clean/table_correspondence_ISIC_NACE.xlsx", c_table_ISIC_NACE, overwrite=true)

# ---------- Match elasticities with NACE Rev.2 -------------------------------------------------------------------------------------------------------

rename!(c_table_ctry, :isic3 => :ISIC3)

# lose too many sectors when matching to 4 digits with adding a zero, therefore, change ISIC into 2 digits and match
#c_table_ctry[:, :ISIC3] = rpad.(c_table_ctry[:, :ISIC3], 4, '1') # add a zero since since we match with 4 digit codes (lose too many sectors)
c_table_ctry[:, :ISIC3] = f.(c_table_ctry[:, :ISIC3])
c_table_ISIC_NACE[:, :ISIC3] = f.(c_table_ISIC_NACE[:, :ISIC3])

elas_NACE = leftjoin(c_table_ctry, c_table_ISIC_NACE, on=:ISIC3)

dropmissing!(elas_NACE) # only two sectors which are not matched (38 and 39 - actually could not find ISIC sectors with these codes!)
begin
    EU27 = ["AUT", "BEL", "BGR", "HRV", "CYP", "CZE", "DNK", "EST", "FIN", "FRA", "DEU", "GRC", "HUN", "IRL", "ITA", "LVA",
    "LTU", "LUX", "MLT", "NLD", "POL", "PRT", "ROU", "SVK", "SVN", "ESP", "SWE"]
end
filter!(row -> row.iso3 in [EU27; "GBR"], elas_NACE) # only keep EU27 and GBR

# summary statistics per country
gdf = groupby(elas_NACE, [:iso3])
d_elas_NACE = combine(gdf, :sigma => mean => :avg, :sigma => maximum => :max, :sigma => minimum => :min) # elasticities only available for 11 countries

# summary statistics per country-sector pairs (some sectors have multiple elasticities - need to average over them)
gdf = groupby(elas_NACE, [:iso3, :WIOD])
elas_NACE = combine(gdf, :sigma => mean => :sigma)

# ---------- Final table containinng all elasticities, countries, sectors -------------------------------------------------------------------------------------------

# use the country-sector elasticities and assume an elasticity of -4 per sector for all other country-sector pairs

WIOD = raw_CPA_NACE.WIOD
unique!(WIOD)

N = size(EU27, 1) + 1 # EU27 + GBR
S = size(WIOD, 1)

country = repeat([EU27; "GBR"], inner=S)
WIOD = repeat(WIOD, outer=N)
d_ctry_WIOD = DataFrame(iso3 = country, WIOD = WIOD)

all_elasticities = leftjoin(d_ctry_WIOD, elas_NACE, on=[:iso3, :WIOD]) # have elasticities for about 15 sectors per available country
all_elasticities.sigma = all_elasticities.sigma .- 1 # elasticities are reported with +1 (investigate this further)
all_elasticities = coalesce.(all_elasticities, 4.0) # assume elasticity of -4 for all the missing country-sector pairs
sort!(all_elasticities, [:iso3, :WIOD])

XLSX.writetable("clean/all_elasticities.xlsx", all_elasticities, overwrite=true)


text:
"""
The ControlUnit.jl module defines the ControlUnit type and methods

------------------------------------------------------------------------------
COPYRIGHT/LICENSE. This file is part of the XYZ package. It is subject to
the license terms in the LICENSE file found in the top-level directory of
this distribution. No part of the XYZ package, including this file, may be
copied, modified, propagated, or distributed except according to the terms
contained in the LICENSE file.
------------------------------------------------------------------------------
"""
# --- Exports

# ------ Types

export ControlUnit
export BuiltInControlSignal
export ControlUnitResponse

# ------ Functions

export run!, process_control_signal!
export is_running, is_terminated
export wait_for_input_ready

# --- Type definitions

using ZMQ

"""
    enum ControlUnitState

Built-in control unit states.
"""

@enum ControlUnitState begin
    NOT_RUNNING
    RUNNING
    TERMINATED
end

"""
    enum ControlUnitResponse

Built-in control unit responses.
"""

@enum ControlUnitResponse begin
    SUCCESS
    FAILED
    UNKNOWN_SIGNAL
end

"""
    enum BuiltInControlSignal

Built-in control unit signals.

TODO: document reserved control signal values, including "nothing"
"""

@enum BuiltInControlSignal begin
    START = -1
    STOP = -2
    TERMINATE = -3
end

BUILT_IN_CONTROL_SIGNALS = instances(BuiltInControlSignal)

"""
    struct ControlUnit

Component that listens for and processes control signals. The `logic_core`
field (1) provides methods to handle control signals and (2) stores the
current control state (if any).
"""
struct ControlUnit
    #=
      Fields
      ------
      * `state`: control unit state

      * `logic_core`: control logic core

      * `url`: URL that ControlUnit listens for control signals on

      * `socket`: ZMQ Socket for receiving control signals and
                          returning status
    =#
    state::Dict
    logic_core::AbstractControlLogicCore
    url::String
    socket::Socket
end

"""
    ControlUnit(logic_core::AbstractControlLogicCore,
                url::String;
                copy_logic_core=true,
                use_bind=true)

Construct a control unit ... TODO
"""
function ControlUnit(logic_core::AbstractControlLogicCore,
                     url::String;
                     copy_logic_core=true,
                     use_bind=true)

    # Initialize state
    state = Dict("mode"=>NOT_RUNNING)

    # Create Socket to listen for control signals
    socket = Socket(REP)

    # Connect socket to URL
    if use_bind
        bind(socket, url)
    else
        connect(socket, url)
    end

    # Return new ControlUnit
    if copy_logic_core
        ControlUnit(state, deepcopy(logic_core), url, socket)
    else
        ControlUnit(state, logic_core, url, socket)
    end
end

# --- Method definitions

"""
    run!(control_unit::ControlUnit)

Start control signal processing loop for `control_unit`.
"""
function run!(control_unit::ControlUnit)
    # Set state to RUNNING
    control_unit.state["mode"] = RUNNING

    # Start control signal processing loop within a Task
    @async begin
        while !is_terminated(control_unit)
            process_control_signal!(control_unit)
        end
    end
end

"""
    process_control_signal!(control_unit::ControlUnit)

Receive and process next control signal.
"""
function process_control_signal!(control_unit::ControlUnit)
    # --- Wait for control signal

    message = recv(control_unit.socket, Vector{UInt8})

    # --- Decode message

    signal = nothing

    # Attempt to interpret as a built-in control signal
    try
        signal = _decode_control_signal(ControlUnit, message)
    catch ArgumentError
        signal = nothing
    end

    # Attempt to interpret as a custom control signal
    if isnothing(signal)
        try
            signal = decode_control_signal(typeof(control_unit.logic_core),
                                           message)
        catch ArgumentError
            signal = nothing
        end
    end

    # --- Process signal

    # Handle unknown and invalid signals
    if isnothing(signal)
        response = UNKNOWN_SIGNAL
        send(control_unit.socket, response)
        return
    end

    # Process valid signal
    if signal in BUILT_IN_CONTROL_SIGNALS
        try
            response = _process_in_control_signal!(control_unit, signal)
        catch
            response = FAILED

            # TODO: add log message
        end
    else
        # Process custom control signals
        try
            response = process_control_signal!(control_unit.logic_core, signal)
        catch
            response = FAILED

            # TODO: add log message
        end
    end

    # --- Send response

    if isnothing(response)
        response = ""
    end
    send(control_unit.socket, response)
end

"""
    is_running(control_unit::ControlUnit)

Return true if the `control_unit` is RUNNING state; return false otherwise.
"""
is_running(control_unit::ControlUnit) = control_unit.state["mode"] == RUNNING

"""
    is_terminated(control_unit::ControlUnit)

Return true if the `control_unit` is TERMINATED state; return false otherwise.
"""
is_terminated(control_unit::ControlUnit) =
    control_unit.state["mode"] == TERMINATED

"""
    wait_for_input_ready(control_unit::ControlUnit,
                         input_channels::Vector{InputChannel})

Wait for input channels to be ready for processing.
"""
function wait_for_input_ready(control_unit::ControlUnit,
                              input_channels::Vector{InputChannel})

    wait_for_input_ready(control_unit.logic_core, input_channels)
end

# --- Private helper functions

"""
    _decode_control_signal(::control_unit::ControlUnit, bytes::Vector{UInt8})

Convert `bytes` to a signal that is understood by `_process_control_signal!()`.
Throw an ArgumentError if `bytes` does not convert to a valid signal.
"""
function _decode_control_signal(::Type{ControlUnit}, bytes::Vector{UInt8})

    signal = reinterpret(BuiltInControlSignal, bytes)[1]

    if !(signal in BUILT_IN_CONTROL_SIGNALS)
        throw(ArgumentError)
    end

    return signal
end

"""
    _process_in_control_signal!(control_unit::ControlUnit, signal)

Process built-in control signals. Return response
"""
function _process_in_control_signal!(control_unit::ControlUnit, signal)

    if signal == START
        control_unit.state["mode"] = RUNNING
    elseif signal == STOP
        control_unit.state["mode"] = NOT_RUNNING
    elseif signal == TERMINATE
        control_unit.state["mode"] = TERMINATED
    else
        return UNKNOWN_SIGNAL
    end

    return SUCCESS
end


text:
using PumasTutorials
PumasTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
function [Grid]=GridToKiteGridNeu(Grid,OrientFace)
KiteGrid.nBar=[ 0  1   0   -1
                1  0  -1    0];
KiteGrid.Dim=3;
KiteGrid.Type='Quad';

for iF=1:Grid.NumFaces
  Grid.Faces(iF).Kite.N=zeros(1,2*size(Grid.Faces(iF).N,2)+1);
  Grid.Faces(iF).Kite.N(1,1:size(Grid.Faces(iF).N,2))=...
    Grid.Faces(iF).N;
  Grid.Faces(iF).Kite.NumNodes=size(Grid.Faces(iF).N,2);
  Grid.Faces(iF).Kite.E=zeros(1,3*size(Grid.Faces(iF).N,2));
  Grid.Faces(iF).Kite.NumEdges=0;
  Grid.Faces(iF).Kite.F=zeros(1,size(Grid.Faces(iF).N,2));
  Grid.Faces(iF).Kite.NumFaces=0;
end

NumNodes=Grid.NumNodes+Grid.NumEdges+Grid.NumFaces;
Nodes(1:NumNodes)=Node([0 0 0],0);
NodeNumber=1;
for iN=1:Grid.NumNodes
  Nodes(NodeNumber)=Node(Grid.Nodes(iN).P,NodeNumber);
  NodeNumber=NodeNumber+1;
end
for iE=1:Grid.NumEdges
  P=0.5*(Grid.Nodes(Grid.Edges(iE).N(1)).P...
    +Grid.Nodes(Grid.Edges(iE).N(2)).P);
  Nodes(NodeNumber)=Node(P,NodeNumber);
  for iFLoc=1:size(Grid.Edges(iE).F,2)
    iF=Grid.Edges(iE).F(iFLoc);
    Grid.Faces(iF).Kite.NumNodes=Grid.Faces(iF).Kite.NumNodes+1;
    Grid.Faces(iF).Kite.N(1,Grid.Faces(iF).Kite.NumNodes)=NodeNumber;
  end
  NodeNumber=NodeNumber+1;
end
for iF=1:Grid.NumFaces
  P=[0 0 0];
  for iN=1:size(Grid.Faces(iF).N,2)
    P=P+Grid.Nodes(Grid.Faces(iF).N(iN)).P;
  end
  P=(1/size(Grid.Faces(iF).N,2))*P;
  Nodes(NodeNumber)=Node(P,NodeNumber);
  Grid.Faces(iF).Kite.NumNodes=Grid.Faces(iF).Kite.NumNodes+1;
  Grid.Faces(iF).Kite.N(1,Grid.Faces(iF).Kite.NumNodes)=NodeNumber;
  NodeNumber=NodeNumber+1;
end
KiteGrid.Nodes=Nodes;
KiteGrid.NumNodes=NumNodes;

NumEdges=2*Grid.NumEdges;
for iF=1:Grid.NumFaces
  NumEdges=NumEdges+size(Grid.Faces(iF).N,2);
end

Edges(1:NumEdges)=Edge([1 2],KiteGrid,0,0);
EdgeNumber=1;
for iE=1:Grid.NumEdges
  N1=Grid.Edges(iE).N(1);
  N2=iE+Grid.NumNodes;
  N3=Grid.Edges(iE).N(2);
  Edges(EdgeNumber)=Edge([N1 N2],KiteGrid,EdgeNumber,EdgeNumber);
  EdgeNumber=EdgeNumber+1;
  Edges(EdgeNumber)=Edge([N2 N3],KiteGrid,EdgeNumber,EdgeNumber);
%   for iFLoc=1:size(Grid.Edges(iE).F,2)
%     iF=Grid.Edges(iE).F(iFLoc);
%     Grid.Faces(iF).Kite.NumEdges=Grid.Faces(iF).Kite.NumEdges+1;
%     Grid.Faces(iF).Kite.E(1,Grid.Faces(iF).Kite.NumEdges)=EdgeNumber-1;
%     Grid.Faces(iF).Kite.NumEdges=Grid.Faces(iF).Kite.NumEdges+1;
%     Grid.Faces(iF).Kite.E(1,Grid.Faces(iF).Kite.NumEdges)=EdgeNumber;
%   end
  EdgeNumber=EdgeNumber+1;
end
for iF=1:Grid.NumFaces
  for i=1:size(Grid.Faces(iF).E,2)
    iE=Grid.Faces(iF).E(i);
    if Grid.Edges(iE).N(1)==Grid.Faces(iF).N(i);
      Grid.Faces(iF).Kite.E(2*i-1)=2*iE-1;
      Grid.Faces(iF).Kite.E(2*i)=2*iE;
    else
      Grid.Faces(iF).Kite.E(2*i-1)=2*iE;
      Grid.Faces(iF).Kite.E(2*i)=2*iE-1;
    end
  end
  Grid.Faces(iF).Kite.NumEdges=2*size(Grid.Faces(iF).E,2);
end
NumInnerNode=Grid.NumNodes+Grid.NumEdges+1;
for iF=1:Grid.NumFaces
  N2=NumInnerNode;
  for iE=1:size(Grid.Faces(iF).E,2)
    N1=Grid.NumNodes+Grid.Edges(Grid.Faces(iF).E(iE)).E;
    Edges(EdgeNumber)=Edge([N1 N2],KiteGrid,EdgeNumber,EdgeNumber);
    Grid.Faces(iF).Kite.NumEdges=Grid.Faces(iF).Kite.NumEdges+1;
    Grid.Faces(iF).Kite.E(1,Grid.Faces(iF).Kite.NumEdges)=EdgeNumber;
    EdgeNumber=EdgeNumber+1;
  end
  NumInnerNode=NumInnerNode+1;
  Grid.Faces(iF).Kite.NumEdgesI=Grid.Faces(iF).Kite.NumEdges;
  Grid.Faces(iF).Kite.NumEdgesB=0;
end

KiteGrid.Edges=Edges;
KiteGrid.NumEdges=NumEdges;
KiteGrid.NumEdgesI=NumEdges;
KiteGrid.NumEdgesB=0;
NumFaces=0;
for iF=1:Grid.NumFaces
  NumFaces=NumFaces+size(Grid.Faces(iF).N,2);
end
Faces(1:NumFaces)=Face([0 0 0 0],KiteGrid,0,'x');
FaceNumber=1;
NumInnerEdge=2*Grid.NumEdges+1;
Type='o';
for iF=1:Grid.NumFaces
  for iE=1:size(Grid.Faces(iF).E,2)-1
    if Grid.Edges(Grid.Faces(iF).E(iE)).N(1)...
        ==Grid.Faces(iF).N(iE+1)
      E1=2*Grid.Edges(Grid.Faces(iF).E(iE)).E-1;
    else
      E1=2*Grid.Edges(Grid.Faces(iF).E(iE)).E;
    end
    if Grid.Edges(Grid.Faces(iF).E(iE+1)).N(1)...
        ==Grid.Faces(iF).N(iE+1)
      E2=2*Grid.Edges(Grid.Faces(iF).E(iE+1)).E-1;
    else
      E2=2*Grid.Edges(Grid.Faces(iF).E(iE+1)).E;
    end
    E3=NumInnerEdge+1;
    E4=NumInnerEdge;
    [Faces(FaceNumber),KiteGrid]=Face([E1 E2 E3 E4],KiteGrid,...
      FaceNumber,Type,OrientFace);
    Grid.Faces(iF).Kite.NumFaces=Grid.Faces(iF).Kite.NumFaces+1;
    Grid.Faces(iF).Kite.F(1,Grid.Faces(iF).Kite.NumFaces)=FaceNumber;
    FaceNumber=FaceNumber+1;
    NumInnerEdge=NumInnerEdge+1;
  end
  iE=size(Grid.Faces(iF).E,2);
  if Grid.Edges(Grid.Faces(iF).E(iE)).N(1)...
      ==Grid.Faces(iF).N(1)
    E1=2*Grid.Edges(Grid.Faces(iF).E(iE)).E-1;
  else
    E1=2*Grid.Edges(Grid.Faces(iF).E(iE)).E;
  end
  if Grid.Edges(Grid.Faces(iF).E(1)).N(1)...
      ==Grid.Faces(iF).N(1)
    E2=2*Grid.Edges(Grid.Faces(iF).E(1)).E-1;
  else
    E2=2*Grid.Edges(Grid.Faces(iF).E(1)).E;
  end
  E3=NumInnerEdge+1-size(Grid.Faces(iF).E,2);
  E4=NumInnerEdge;
  [Faces(FaceNumber),KiteGrid]=Face([E1 E2 E3 E4],KiteGrid,...
    FaceNumber,Type,OrientFace);
  Grid.Faces(iF).Kite.NumFaces=Grid.Faces(iF).Kite.NumFaces+1;
  Grid.Faces(iF).Kite.F(1,Grid.Faces(iF).Kite.NumFaces)=FaceNumber;
  FaceNumber=FaceNumber+1;
  NumInnerEdge=NumInnerEdge+1;
end
KiteGrid.Faces=Faces;
KiteGrid.NumFaces=NumFaces;

KiteGrid=Orientation(KiteGrid);
KiteGrid=Renumbering(KiteGrid);
Grid.KiteGrid=KiteGrid;
end


text:
using Test
using LogicCircuits
using LogicCircuits.LoadSave: zoo_lc_file, zoo_vtree_file

include("../helper/plain_logic_circuits.jl")

@testset "Circuit saver tests" begin

  mktempdir() do tmp
    circuit, vtree = load_struct_smooth_logic_circuit(zoo_lc_file("mnist-medium.circuit"), zoo_vtree_file("balanced.vtree"))
    # the circuit above does not conform the the SDD requirements
    @test_throws Exception save_circuit("$tmp/temp.sdd", circuit, vtree)
  end

  mktempdir() do tmp
    mgr = SddMgr(7, :balanced)
    v = Dict([(i => compile(mgr, Lit(i))) for i=1:7])
    c = (v[1] | !v[2] | v[3]) &
        (v[2] | !v[7] | v[6]) &
        (v[3] | !v[4] | v[5]) &
        (v[1] | !v[4] | v[6])
    @test num_edges(c) == 147

    @test_nowarn save_circuit("$tmp/temp.sdd", c, mgr)
    l = load_logic_circuit("$tmp/temp.sdd")
    @test num_edges(l) == 147
    @test num_variables(l) == 7

    c2 = PlainStructLogicCircuit(mgr,c)

    @test_nowarn save_circuit("$tmp/temp2.sdd", c2, mgr)
    l = load_logic_circuit("$tmp/temp2.sdd")
    @test num_edges(l) == 147
    @test num_variables(l) == 7
  end

  mktempdir() do tmp
    circuit = fully_factorized_circuit(LogicCircuit,10)
    @test_nowarn save_as_dot("$tmp/temp.dot", circuit)
  end

  mktempdir() do tmp
    circuit = little_3var_constants()
    @test_nowarn save_as_dot("$tmp/temp.dot", circuit)
  end

  mktempdir() do tmp
    a = compile(PlainLogicCircuit, Lit(1))
    t = compile(PlainLogicCircuit, true)
    na = compile(PlainLogicCircuit, Lit(-1))
    nb = compile(PlainLogicCircuit, Lit(-2))
    a′ = compile(PlainLogicCircuit, Lit(1))
    f = compile(PlainLogicCircuit, false)
    na′ = compile(PlainLogicCircuit, Lit(-1))
    nb′ = compile(PlainLogicCircuit, Lit(-2))
    and1 = a & nb
    and2 = na & t
    and3 = a′ & f
    and4 = na′ & nb′
    or = disjoin(and1, and2, and3, and4)
    @test_nowarn save_as_tex("$tmp/temp.tex", or)
    @test_nowarn save_as_dot2tex("$tmp/temp.tex", or)
  end

  # TODO add a test to load and save and load an .sdd file
  # currently we have no sdd+vtree in the model zoo to do this

end


text:
printstyled("  safe_isfile\n", color=:light_green)
@test safe_isfile("runtests.jl") == true
@test safe_isfile("foo.jl") == false

printstyled("  safe_isdir\n", color=:light_green)
@test safe_isdir("SampleFiles") == true
@test safe_isdir("Roms") == false

printstyled("  ls\n", color=:light_green)

@test any([occursin("test", i) for i in ls()])
coreutils_flist = [ "test_calculus.jl",
                    "test_get_svn.jl",
                    "test_ls.jl",
                    "test_FastIO.jl",
                    "test_poly.jl",
                    "test_read_utils.jl",
                    "test_time.jl",
                    "test_ts_timekeeping.jl",
                    "test_typ2code.jl" ]

S = [
      "CoreUtils/",
      "CoreUtils/test_ls.jl",
      "CoreUtils/poo",
      "SampleFiles/UW/*W",
      "SampleFiles/UW/02*o",
      "CoreUtils/test_*"
    ]
S_expect =  [
              coreutils_flist,
              ["test_ls.jl"],
              String[],
              ["00012502123W", "99011116541W"],
              ["02062915175o", "02062915205o"],
              coreutils_flist,
            ]

# Test that ls returns the same files as `ls -1`
for (n,v) in enumerate(S)
  files = String[splitdir(i)[2] for i in ls(v)]
  # if Sys.iswindows() == false
    expected = S_expect[n]
    @test sort(files) == sort(expected)
  # end
  [@test isfile(f) for f in ls(v)]
end
println("arch = ", Sys.ARCH)
arch_allowed = Sys.ARCH == :x86_64
# Arm64 = :aarch64

# Test that ls invokes find_regex under the right circumstances
if arch_allowed
  @test change_sep(ls(S[5])) == change_sep(regex_find("SampleFiles/", r"02.*o$"))
end

if has_restricted && arch_allowed
  println("extended ls tests")
  T = path .* [
                "/SampleFiles/*",
                "/SampleFiles/Restricted/*.cnt",
                "/SampleFiles/Restricted/2014092709*cnt"
              ]
  T_expect =  [540, 63, 60]

  # Test that ls finds the same number of files as bash `ls -1`
  for (n,v) in enumerate(T)
    files = ls(v)
    @test (isempty(files) == false)
    if n == 1
      @test length(files) ≥ T_expect[n]
    else
      @test length(files) == T_expect[n]
    end
    [@test isfile(f) for f in files]
  end

  # Test that ls invokes find_regex under the right circumstances
  @test change_sep(ls(T[1])) == change_sep(regex_find("SampleFiles", r".*$"))
  @test change_sep(ls(T[3])) == change_sep(regex_find("SampleFiles", r"Restricted/2014092709.*cnt$"))
else
  printstyled("  extended ls tests skipped. (files not found; is this Appveyor?)\n", color=:green)
end

if Sys.iswindows()
  @test safe_isfile("http://google.com") == false
  @test safe_isdir("http://google.com") == false
end


text:
import Dolo: MSM

"""
Computes the residuals of the arbitrage equations. The general form of the arbitrage equation is

    `0 = E_t [f(m, s, x, M, S, X; p)]`

where `m` are current exogenous variables, `s` are current states,
`x` are current controls, `M` are next period's exogenous variables, `S` are next period's states, `X` are next period's controls, and `p` are the model parameters. This function evaluates the right hand side of the arbitrage equation for the given inputs.

If the list of current controls `x` is provided as a two-dimensional array (`ListOfPoints`), it is transformed to a one-dimensional array (`ListOfListOfPoints`).

# Arguments
* `model::NumericModel`: Model object that describes the current model environment.
* `dprocess`: Discretized exogenous process.
* `s::ListOfPoints`: List of state variable values.
* `x::ListOfListOfPoints`: List of control variable values associated with each exogenous shock.
* `p::Vector{Float64}`: Model parameters.
* `dr`: Current guess for the decision rule.
# Returns
* `res`: Residuals of the arbitrage equation associated with each exogenous shock.
"""
# function euler_residuals_ti!(res::MSM{Point{n_x}}, model, dprocess::AbstractDiscretizedProcess,s::ListOfPoints{d}, x::MSM{Point{n_x}}, p::SVector, dr) where n_x where d

#     N = length(s)

#     for i in 1:length(res)
#         m = node(Point, dprocess, i)
#         for (w, M, j) in get_integration_nodes(Point, dprocess,i)
#             # Update the states
#             # TODO: replace views here
#             S = Dolo.transition(model, m, s, x.views[i], M, p)
#             X = dr(i, j, S)
#             res.views[i][:] += w*Dolo.arbitrage(model, m, s, x.views[i], M, S, X, p)
#         end
#     end
#     return res
# end

# function euler_residuals_ti(model, dprocess::AbstractDiscretizedProcess,s::ListOfPoints{d}, x::Vector{Vector{Point{n_x}}}, p::SVector, dr) where n_x where d
#     # res = deepcopy(x)
#     # for i_m=1:length(res)
#     #     res[i_m][:] *= 0.0
#     # end
#     xx = MSM(x)
#     res = copy(xx)
#     reset!(res)
#     euler_residuals_ti!(res, model, dprocess, s, xx, p, dr)
#     return vecvec(res)
# end

struct Euler{ID, n_m, n_s, n_x, Gx, Ge}

    model::AbstractModel
    p:: SVector  ## SVector
    grid::ProductGrid{Gx, Ge}
    grid_endo # obsolete
    grid_exo # obsolete
    dprocess
    s0::ListOfPoints{n_s}
    x0::MSM{SVector{n_x, Float64}}
    bounds::Union{Nothing,Tuple{MSM{SVector{n_x, Float64}},MSM{SVector{n_x, Float64}}}}
    dr::CachedDecisionRule
    cdr::CachedDecisionRule

    function Euler(model;
        discretization=Dict(),
        interpolation=:cubic,
        dr0=Dolo.ConstantDecisionRule(model),
        ignore_constraints=false
    )

        p = SVector(model.calibration[:parameters]...)

        grid, dprocess = discretize(model; discretization...)
        grid_endo = grid.endo
        grid_exo = grid.exo

        x = SVector(model.calibration[:controls]...)

        N_m = max(Dolo.n_nodes(grid_exo),1)
        N_s = Dolo.n_nodes(grid_endo)

        n_m = length(model.symbols[:exogenous])
        n_s = length(model.symbols[:states])
        n_x = length(model.symbols[:controls])

        xxx = [ [dr0(j, Dolo.node(grid_endo, i)) for i=1:N_s] for j=1:N_m]

        x0 = MSM(xxx)

        s0 = Dolo.nodes(grid_endo)

        dr = Dolo.CachedDecisionRule(dprocess, grid_endo, xxx)
        cdr = Dolo.CachedDecisionRule(dprocess, grid_endo, xxx)

        m0 = SVector(model.calibration[:exogenous]...)

        grid = ProductGrid(grid_exo, grid_endo)
        Gx = typeof(grid.exo)
        Ge = typeof(grid.endo)

        if (:controls_lb in keys(model.factories)) & !ignore_constraints
            if Dolo.n_nodes(grid_exo)>0
                lb = [Dolo.controls_lb(model, m, s0, p) for m in Dolo.nodes(grid_exo)]
                ub = [Dolo.controls_ub(model, m, s0, p) for m in Dolo.nodes(grid_exo)]
            else
                lb = [Dolo.controls_lb(model, m0, s0, p) ]
                ub = [Dolo.controls_ub(model, m0, s0, p) ]
            end
            bounds = (MSM(lb), MSM(ub))
        else
            bounds = nothing
        end

        ID = id(model)

        new{ID, n_m, n_s, n_x, Gx, Ge}(model, p, grid, grid_endo, grid_exo, dprocess, s0, x0, bounds, dr, cdr)

    end

end

repsvec(v1::Point{d1}, v2::Point{d}) where d1 where d= SVector(v1..., v2[d1+1:end]...)

function euler_residuals(F::Euler{id, n_m, n_s, n_x}, s::Vector{Point{n_s}}, x::MSM{Point{n_x}}, dr,
    dprocess, parms::SVector; exo=nothing, keep_J_S=false, out=nothing) where id where n_s where n_x where n_m #, jres=nothing, S_ij=nothing)
    #
    # if set_dr ==true
    #   set_values!(dr,x)
    # end

    model = F.model

    if out === nothing
        res = zeros_like(x)::MSM{Point{n_x}}
    else
        res = out
    end

    res.data .*= 0.0 # just to be on the safe side...

    N_s = length(s) # Number of gris points for endo_var

    n_ms = length(x.views)  # number of exo states today
    n_mst = n_inodes(dprocess,1)  # number of exo states tomorrow
    d = length(s[1])

    # TODO: allocate properly...

    if keep_J_S
        jres = zeros((n_ms,n_mst,N_s,n_x,n_x))
        fut_S = zeros((n_ms,n_mst,N_s,n_s))
        J_ij = Vector{SMatrix{n_x,n_x,Float64,n_x*n_x}}[to_LOJ(jres[i,j,:,:,:]) for i=1:size(jres,1), j=1:size(jres,2)]
        S_ij = Vector{Point{d}}[to_LOP(fut_S[i,j,:,:]) for i=1:size(fut_S,1), j=1:size(fut_S,2)]
    end

    for i_ms in 1:n_ms
        m = node(Point{n_m}, dprocess,i_ms)
        if !(exo === nothing)
            m = repsvec(exo[1], m)   # z0
        end
        xx = x.views[i_ms]
        for I_ms in 1:n_mst
            M = inode(Point{n_m}, dprocess, i_ms, I_ms)
            if !(exo === nothing)
                M = repsvec(exo[2], M)   # z1
            end
            w = iweight(dprocess, i_ms, I_ms)
            S = transition(model, m, s, xx, M, parms)
            X = dr(i_ms, I_ms, S)
            if keep_J_S==true
                rr, rr_XM = arbitrage(model,Val{(0,6)},m,s,xx,M,S,X,parms)
                J_ij[i_ms,I_ms][:] = w*rr_XM
                S_ij[i_ms,I_ms][:] = S
                res.views[i_ms][:] += w*rr
            else
                rr = arbitrage(model, m, s, xx, M, S, X, parms)
                res.views[i_ms][:] .+= w*rr
            end
        end
    end

    # TODO: make this function type-stable
    if keep_J_S
        return res,J_ij,S_ij
    else
        return res
    end
end

function (F::Euler)(x0::MSM, x1::MSM, z0::Point{n_z}, z1::Point{n_z}; exo=nothing, kwargs...) where n_z
    F(x0, x1; exo=(z0,z1), kwargs...)
end

function (F::Euler)(x0::MSM, x1::MSM; exo=nothing, set_future=true, ignore_constraints=false, out=nothing, p=F.p)

    if out === nothing
        rr = copy(x0)
    else
        rr = out
    end

    if set_future
        set_values!(F.dr, x1)
    end

    euler_residuals(F,F.s0,x0,F.dr,F.dprocess,p;exo=exo,keep_J_S=false,out=rr)

    if (F.bounds!==nothing) & !ignore_constraints
        lb, ub = F.bounds
        for n=1:length(rr.data)
            rr.data[n] = PhiPhi0(rr.data[n], x0.data[n], lb.data[n], ub.data[n])
        end
    end

    return rr

end

function df_A(F, z0, z1; set_future=false, inplace=false, ret_res=false, exo=nothing)

    if set_future
        set_values!(F.dr, z1)
    end

    if !inplace
        fun  = z->F(z, z1; set_future=false, ignore_constraints=true, exo=exo)
        f0,J = Dolo.DiffFun(fun, z0, 1e-8)

    else
        n_x = length(z0.data[1])

        f0 = z0*0
        fi = z0*0
        xi = z0*0

        N = length(xi.data)
        cata = zeros(SMatrix{n_x, n_x, Float64, n_x*n_x}, N)
        J = MSM(cata, f0.sizes)
        out = f0, fi, xi, J

        fun!  = (r,z)->F(z, z1; set_future=false, ignore_constraints=true, out=r)

        f00, J0 = Dolo.DiffFun!(fun!, z0, 1e-8, out)

        @assert f00 === f0
        @assert J === J0

    end
    if (F.bounds!==nothing)
        lb, ub = F.bounds
        for n=1:length(f0.data)
            z,z_f,z_x  = PhiPhi(f0.data[n], z0.data[n], lb.data[n], ub.data[n])
            f0.data[n] = z
            J.data[n]= z_f*J.data[n] + z_x
        end
    end

    if ret_res
        return (f0,J)
    else
        return J
    end
end

function df_B(F, x0, x1; set_future=false, exo=nothing)

    ddr_filt = F.cdr

    # res = deepcopy(z0)

    if set_future
        set_values!(F.dr, x1)
    end

    rr,J_ij,S_ij  = euler_residuals(F, F.s0, x0 , F.dr, F.dprocess, F.p; keep_J_S=true, exo=exo)

    if (F.bounds!==nothing)
        lb, ub = F.bounds
        PhiPhi!(rr.views, x1.views, lb.views, ub.views, J_ij)
    end

    L = LinearThing(J_ij, S_ij, ddr_filt)

    return L

end

function DiffFun2(fun, e::SVector{n_p}, ε=1e-6) where n_p

    # TODO: this could be done without allocation
    r0 = fun(e)
    n_x = length(r0.data[1])

    diffs = []
    for i=1:n_p
        ei = SVector([(j==i) ? e[i]+ε : e[j] for j=1:n_p ]...)
        d_i = (fun(ei)-r0)/ε
        push!(diffs, d_i)
    end

    res = [
        SMatrix{n_x, n_p, Float64, n_x*n_p}(cat([diffs[i].data[n] for i=1:n_p]...;dims=2))
        for n=1:length(r0.data)
    ]

    return MSM(res, r0.sizes)

end

function df_e(F, x0, x1, e1, e2; set_future=false)

    ε = 1e-8

    # diff w.r.t. e1
    # diff w.r.t. e2

    tt = F(x0,x1; exo=(e1,e2), set_future=false)

    # Awful hack
    fun1 = e -> F(x0,x1; exo=(e,e2))

    J1 = Dolo.DiffFun2(fun1, e1, 1e-8)

    fun2 = e -> F(x0,x1; exo=(e1,e))
    J2 = Dolo.DiffFun2(fun2, e2, 1e-8)

    return J1, J2

end

function prediv!(L::LinearThing,x::MSM)
    for i=1:size(L.M_ij,1)
        N = length(x.views[i])
        for j=1:size(L.M_ij,2)
            for n=1:N
                L.M_ij[i,j][n] = x.views[i][n] \ L.M_ij[i,j][n]
            end
        end

    end
    return L
end

function premult!(L::LinearThing,x::MSM)
    for i=1:size(L.M_ij,1)
        N = length(x.views[i])
        for j=1:size(L.M_ij,2)
            for n=1:N
                L.M_ij[i,j][n] = x.views[i][n] * L.M_ij[i,j][n]
            end
        end

    end
    return L
end

function mult!(L::LinearThing,x::T) where T<:Number
    for i=1:size(L.M_ij,1)
        for j=1:size(L.M_ij,2)
            L.M_ij[i,j][:] *= x
        end

    end
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

@recipe function f(problem::LearningProblem)
  sdata   = sourcedata(problem)
  tdata   = targetdata(problem)
  svars   = join(sort(collect(name.(variables(sdata)))), ", ")
  tvars   = join(sort(collect(name.(variables(tdata)))), ", ")

  title --> "Learning Problem"
  legend --> true

  @series begin
    seriescolor --> :green
    label --> "source data ($svars)"
    marker --> :xcross
    domain(sdata)
  end
  @series begin
    seriescolor --> :blue
    label --> "target data ($tvars)"
    marker --> :xcross
    domain(tdata)
  end
end


text:
# ******************************************************************************************
# Notices:
#
# Copyright © 2022 United States Government as represented by the Administrator of the
# National Aeronautics and Space Administration.  All Rights Reserved.
#
# Disclaimers
#
# No Warranty: THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY KIND,
# EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY
# THAT THE SUBJECT SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY
# WARRANTY THAT THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT DOCUMENTATION,
# IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE. THIS AGREEMENT DOES NOT, IN ANY MANNER,
# CONSTITUTE AN ENDORSEMENT BY GOVERNMENT AGENCY OR ANY PRIOR RECIPIENT OF ANY RESULTS,
# RESULTING DESIGNS, HARDWARE, SOFTWARE PRODUCTS OR ANY OTHER APPLICATIONS RESULTING FROM
# USE OF THE SUBJECT SOFTWARE.  FURTHER, GOVERNMENT AGENCY DISCLAIMS ALL WARRANTIES AND
# LIABILITIES REGARDING THIRD-PARTY SOFTWARE, IF PRESENT IN THE ORIGINAL SOFTWARE, AND
# DISTRIBUTES IT "AS IS."
#
# Waiver and Indemnity:  RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED
# STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT.
# IF RECIPIENT'S USE OF THE SUBJECT SOFTWARE RESULTS IN ANY LIABILITIES, DEMANDS, DAMAGES,
# EXPENSES OR LOSSES ARISING FROM SUCH USE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON,
# OR RESULTING FROM, RECIPIENT'S USE OF THE SUBJECT SOFTWARE, RECIPIENT SHALL INDEMNIFY AND
# HOLD HARMLESS THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL
# AS ANY PRIOR RECIPIENT, TO THE EXTENT PERMITTED BY LAW.  RECIPIENT'S SOLE REMEDY FOR ANY
# SUCH MATTER SHALL BE THE IMMEDIATE, UNILATERAL TERMINATION OF THIS AGREEMENT.
# ******************************************************************************************

const TODO_COUNTER = Threads.Atomic{Int64}(0)
const DONE_COUNTER = Threads.Atomic{Int64}(0)

add_todo(n::Int64=1) = Threads.atomic_add!(TODO_COUNTER, n)

function add_done(n::Int64=1)
    Threads.atomic_add!(DONE_COUNTER, n)
    Threads.atomic_sub!(TODO_COUNTER, n)
end

"""
Cell information pertaining to proof and processing status.
"""
Base.@kwdef mutable struct CellStatus
    depth::Int64 = 0
    detached::Bool = false
    hash::UInt64 = hash(0)
    member::Bool = false
    proven::Bool = false
    pid::Int64 = 0
end

"""
    BinaryRefinery <: AbstractRefinery

Binary refinement process.
"""
Base.@kwdef struct BinaryRefinery <: AbstractRefinery
    network::Network
    solver = ReluVal(max_iter=100)
    val::Float64 = 0.0                # boundary value
    output = HalfSpace([1.0], val)    # a⋅x ≤ b (region of interest)
    coutput = HalfSpace([-1.0], -val) # a⋅x ≥ b (complement region)
    tol::Float64 = 0.1                # maximum unproven length scale
    ks::Vector{Int64} = []            # detachment depths
end

"""
    IntervalRefinery <: AbstractRefinery

Interval refinement process.
"""
Base.@kwdef struct IntervalRefinery <: AbstractRefinery
    network::Network
    solver = ReluVal(max_iter=100)
    interval::Vector{Float64} = [-1.0, 1.0]
    output = Hyperrectangle(low=interval[1:1], high=interval[2:2])
    tol::Float64 = 0.1
    ks::Vector{Int64} = []
end

"""
Determine whether condition can be proven true or false over entire cell
or if further refinement is needed.
"""
function needs_refinement(cell::Cell, r::AbstractRefinery)
    # Checks status.
    cell.data.proven && return false

    # Checks tolerance.
    w = cell.boundary.widths
    maximum(w) < r.tol && return false

    # Determines corresponding input set.
    l = cell.boundary.origin
    h = l + w
    input = Hyperrectangle(low=convert(Vector{Float64}, l), high=convert(Vector{Float64}, h))

    # Attempts verification of condition.
    problem = Problem(r.network, input, r.output)
    result = solve(r.solver, problem)
    if result.status == :holds
        cell.data.member = true
        cell.data.proven = true
        return false
    end

    if r isa BinaryRefinery
        # Attempts verification of complement.
        problem = Problem(r.network, input, r.coutput)
        result = solve(r.solver, problem)
        if result.status == :holds
            cell.data.member = false
            cell.data.proven = true
            return false
        end
    end

    return true
end

"""
Generate data for newly-created cell.
"""
function child_data(c::Cell, idx::Tuple)
    return CellStatus(depth = c.data.depth + 1, hash = hash((c.data.hash, idx)))
end

"""
    refine!(cell::Cell, r::AbstractRefinery)

Recursively build miniminal k-d tree that defines volumes proven to satisfy given condition
and its complement. This version is single-process.
"""
function refine!(cell::Cell, r::AbstractRefinery)
    add_todo()
    yield()
    if needs_refinement(cell, r)
        split!(cell, child_data)
        for c in children(cell)
           refine!(c, r)
        end
        merge!(cell)
    end
    add_done()
    yield()
end

"""
Merge children of cell if all are proven and belong to the same set, or if all are unproven
and attached (not awaiting computation by another process). If `recursive`, all descendents
are merged.
"""
function merge!(cell::Cell; recursive::Bool=false)
    if !isleaf(cell)
        cs = children(cell)
        provens = [recursive ? merge!(c; recursive=true) : c.data.proven for c in cs]
        members = (c -> c.data.member).(cs)
        detachs = (c -> c.data.detached).(cs)
        if all(provens) && reduce(==, members) || !any(provens) && !any(detachs) && all(isleaf.(cs))
            cell.data.member = members[1]
            cell.data.proven = provens[1]
            cell.children = nothing
        end
    end
    return cell.data.proven
end

"""
Verify that tree is well-formed after recombination; i.e., no cells are marked as detached
and all children recognize the correct parent.
"""
function verify(cell::Cell)
    for c in allcells(cell)
        if c.data.detached
            return false
        end
        if !isleaf(c)
            if !all((x->x.parent == c).(children(c)))
                return false
            end
        end
    end
    return true
end

"""
Generate root from region limits.
"""
function get_root(limits::NTuple{2, Vector{<:Real}})
    lows, highs = limits
    widths = highs - lows
    root = Cell(SVector(lows...), SVector(widths...), CellStatus())
    return root
end

"""
Perform PVV analysis, using multiple processes if available.
"""
function analyze(r::AbstractRefinery, limits::Tuple{Vector, Vector}; progress::Bool=true, multiproc::Bool=true)
    TODO_COUNTER[] = 1
    tree = get_root(limits)

    @sync begin
        # Refinement process
        @async begin
            try
                if nprocs() == 1 || !multiproc
                    refine!(tree, r)
                else
                    refine_multiprocess!(tree, r)
                end
            finally
                TODO_COUNTER[] = 0
            end
        end

        # Progress meter
        if progress
            @async begin
                p = ProgressUnknown("Cells awaiting processing:")
                while TODO_COUNTER[] > 0
                    ProgressMeter.update!(p, TODO_COUNTER[]; ignore_predictor=true)
                    sleep(0.1)
                end
                ProgressMeter.update!(p, 0)
                ProgressMeter.finish!(p)
            end
        end
    end

    try
        print_metrics(tree)
    catch
        println("Unable to calculate metrics.")
    end

    return tree
end


text:
using LazySets: CachedMinkowskiSumArray

"""
    check_property(S, property, options)

Interface to property checking algorithms for an LTI system.

### Input

- `S`        -- LTI system, discrete or continuous
- `property` -- property
- `options`  -- additional options

### Notes

A dictionary with available algorithms is available via
`available_algorithms_check`.
"""
function check_property(S::IVP{<:AbstractDiscreteSystem},
                        property::Predicate,
                        options::TwoLayerOptions
                       )::Int
    # optional matrix conversion
    S = matrix_conversion(S, options)

    # list containing the arguments passed to any reachability function
    args = []

    # coefficients matrix
    A = S.s.A
    push!(args, A)

    # determine analysis mode (sparse/dense) for lazy_expm mode
    if A isa SparseMatrixExp
        push!(args, Val(options[:assume_sparse]))
    end

    n = statedim(S)
    blocks = options[:blocks]
    partition = convert_partition(options[:partition])
    T = options[:T]
    N = (T == Inf) ? nothing : ceil(Int, T / options[:δ])

    # Cartesian decomposition of the initial set
    if length(partition) == 1 && length(partition[1]) == n &&
            options[:block_options_init] == LinearMap
        info("- Skipping decomposition of X0")
        Xhat0 = [S.x0]
    else
        info("- Decomposing X0")
        @timing begin
            Xhat0 = array(decompose(S.x0, options[:partition],
                                    options[:block_options_init]))
        end
    end

    # shortcut if only the initial set is required
    if N == 1
        if length(blocks) == 1
            Xhat0_mod = Xhat0[blocks[1]]
        else
            Xhat0_mod = CartesianProductArray(Xhat0[blocks])
        end
        return property(Xhat0_mod) ? 0 : 1
    end
    push!(args, Xhat0)

    # inputs
    if !options[:assume_homogeneous] && inputdim(S) > 0
        U = inputset(S)
    else
        U = nothing
    end
    push!(args, U)

    # overapproximation function for states
    block_options_iter = options[:block_options_iter]
    if block_options_iter isa AbstractVector ||
            block_options_iter isa Dict{Int, Any}
        # individual overapproximation options per block
        overapproximate_fun = (i, X) -> overapproximate(X, block_options_iter[i])
    else
        # uniform overapproximation options for each block
        overapproximate_fun = (i, X) -> overapproximate(X, block_options_iter)
    end

    # overapproximate function for inputs
    lazy_inputs_interval = options[:lazy_inputs_interval]
    if lazy_inputs_interval == lazy_inputs_interval_always
        overapproximate_inputs_fun = (k, i, x) -> overapproximate_fun(i, x)
    else
        # first set in a series
        function _f(k, i, x::LinearMap{NUM}) where {NUM}
            @assert k == 1 "a LinearMap is only expected in the first iteration"
            return CachedMinkowskiSumArray(LazySet{NUM}[x])
        end
        # further sets of the series
        function _f(k, i, x::MinkowskiSum{NUM, <:CachedMinkowskiSumArray}) where NUM
            if has_constant_directions(block_options_iter, i)
                # forget sets if we do not use epsilon-close approximation
                forget_sets!(x.X)
            end
            push!(array(x.X), x.Y)
            if lazy_inputs_interval(k)
                # overapproximate lazy set
                y = overapproximate_fun(i, x.X)
                return CachedMinkowskiSumArray(LazySet{NUM}[y])
            end
            return x.X
        end
        function _f(k, i, x)
            # other set types
            if lazy_inputs_interval(k)
                # overapproximate lazy set
                return overapproximate_fun(i, x.X)
            end
            return x
        end
        overapproximate_inputs_fun = _f
    end
    push!(args, overapproximate_inputs_fun)

    # ambient dimension
    push!(args, n)

    # number of computed sets
    push!(args, N)

    # add mode-specific block(s) argument
    algorithm = options[:algorithm]
    if algorithm == "explicit"
        push!(args, blocks)
        push!(args, partition)
        algorithm_backend = "explicit_blocks"
    else
        error("Unsupported algorithm: ", algorithm)
    end

    # add eager/lazy checking option
    push!(args, options[:eager_checking])

    # add property
    push!(args, property)

    info("- Computing successors")

    # progress meter
    progress_meter = (N != nothing) ?
        Progress(N, 1, "Computing successors ") :
        nothing
    push!(args, progress_meter)

    # call the adequate function with the given arguments list
    answer =
        @timing available_algorithms_check[algorithm_backend]["func"](args...)

    # return the result
    return answer
end

function check_property(problem::IVP{<:AbstractContinuousSystem},
                        property::Predicate,
                        options::TwoLayerOptions
                       )::Int
    info("Time discretization...")
    Δ = @timing discretize(problem, options[:δ],
        algorithm=options[:discretization], exp_method=options[:exp_method],
        sih_method=options[:sih_method])
    return check_property(Δ, property, options)
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/node_bb.jl
# Defines storage for a node in the B&B tree & utilities functions
#############################################################################

"""
$(TYPEDEF)

Stores information associated with each node in Branch & Bound tree.

$(TYPEDFIELDS)
"""
struct NodeBB
    "Lower bounds of variable box."
    lower_variable_bounds::Vector{Float64}
    "Upper bounds of variable box."
    upper_variable_bounds::Vector{Float64}
    "Lower bound of problem solution on nodeBB"
    lower_bound::Float64
    "Upper bound of problem solution on nodeBB"
    upper_bound::Float64
    "Depth of node in B&B tree."
    depth::Int64
    "Unique id for each node."
    id::Int64
end

# Constructors
NodeBB() = NodeBB(Float64[], Float64[], -Inf, Inf, 0, 1)
NodeBB(x::NodeBB) = NodeBB(x.lower_variable_bounds, x.upper_variable_bounds,
                           x.lower_bound, x.upper_bound, x.depth, x.id)

# Copy utilities
Base.copy(x::NodeBB) = NodeBB(copy(x.lower_variable_bounds),
                              copy(x.upper_variable_bounds),
                              x.lower_bound, x.upper_bound, x.depth, x.id)

# using alternative name as to not interfere with ordering...
function uninitialized(x::NodeBB)
    flag = isempty(x.lower_variable_bounds)
    flag &= isempty(x.upper_variable_bounds)
    flag &= x.lower_bound === -Inf
    flag &= x.upper_bound === Inf
    flag &= x.depth === 0
    flag &= x.id === 1
    return flag
end

# Access functions for broadcasting data easily
lower_variable_bounds(x::NodeBB) = x.lower_variable_bounds
upper_variable_bounds(x::NodeBB) = x.upper_variable_bounds
lower_variable_bounds(x::NodeBB, id::Int64, nid::Int64) = x.lower_variable_bounds[id:nid]
upper_variable_bounds(x::NodeBB, id::Int64, nid::Int64) = x.upper_variable_bounds[id:nid]
lower_bound(x::NodeBB) = x.lower_bound
upper_bound(x::NodeBB) = x.upper_bound
depth(x::NodeBB) = x.depth

# Iterations Functions
Base.isless(x::NodeBB, y::NodeBB) = x.lower_bound < y.lower_bound
Base.length(x::NodeBB) = length(x.lower_variable_bounds)
function Base.isempty(x::NodeBB)
    for i = 1:length(x)
        @inbounds lower = x.lower_variable_bounds[i]
        @inbounds upper = x.upper_variable_bounds[i]
        (lower > upper) && (return true)
    end
    return false
end

"""
$(FUNCTIONNAME)

Checks that node `x` and `y` have equal domains withing a tolerance of `atol`.
"""
function same_box(x::NodeBB, y::NodeBB, r::Float64)
    (isempty(x.lower_variable_bounds) ⊻ isempty(y.lower_variable_bounds)) && (return false)
    (isempty(x.upper_variable_bounds) ⊻ isempty(y.upper_variable_bounds)) && (return false)
    for i = 1:length(x)
        ~isapprox(x.lower_variable_bounds[i], y.lower_variable_bounds[i], atol=r) && (return false)
        ~isapprox(x.upper_variable_bounds[i], y.upper_variable_bounds[i], atol=r) && (return false)
    end
    return true
end

# Compute middle & diameter
diam(x::NodeBB) = x.upper_variable_bounds - x.lower_variable_bounds
mid(x::NodeBB) = 0.5*(x.upper_variable_bounds + x.lower_variable_bounds)


text:
using EarthMoversDistance: emd!, libemd_emd, emd, emd_flows

W1 = [0.4, 0.2, 0.2, 0.1, 0.1]
W2 = [0.6, 0.2, 0.2]
Cost1 = ones(Float64, 5, 3)

@testset "Calling the C function emd directly" begin
    flows = zeros(Float64, 5, 3)
    distance = libemd_emd(length(W1), W1, length(W2), W2, Cost1, flows)
    @test distance == 1.0
    @test isapprox(flows[1, 1], 0.4)
    @test isapprox(flows[2, 1], 0.2)
    @test isapprox(flows[3, 2], 0.2)
    @test isapprox(flows[4, 3], 0.1)
    @test isapprox(flows[5, 3], 0.1)
    @test Cost1 == ones(Float64, 5, 3)
end

@testset "emd" begin
    distance = emd(W1, W2, Cost1)
    @test distance == 1.0
    @test Cost1 == ones(Float64, 5, 3)
end

@testset "emd_flows" begin
    distance, flows = emd_flows(W1, W2, Cost1)
    @test distance == 1.0
    @test isapprox(flows[1, 1], 0.4)
    @test isapprox(flows[2, 1], 0.2)
    @test isapprox(flows[3, 2], 0.2)
    @test isapprox(flows[4, 3], 0.1)
    @test isapprox(flows[5, 3], 0.1)
    @test Cost1 == ones(Float64, 5, 3)
end

@testset "python-emd example1" begin
    # This corresponds to the example in the file: https://github.com/pdinges/python-emd/blob/master/example1.py

    # features1 = [Feature(100, 40, 22), Feature(211, 20, 2),
    #             Feature(32, 190, 150), Feature(2, 100, 100)]
    # weights1  = [0.4, 0.3, 0.2, 0.1]
    features1 = [[100, 40, 22], [211, 20, 2], [32, 190, 150], [2, 100, 100]]
    weights1  = [0.4, 0.3, 0.2, 0.1]

    #features2 = [Feature(0, 0, 0), Feature(50, 100, 80), Feature(255, 255, 255)]
    #weights2  = [0.5, 0.3, 0.2]
    features2 = [[0, 0, 0], [50, 100, 80], [255, 255, 255]]
    weights2  = [0.5, 0.3, 0.2]

    # The costs are the euclidean distances between all pairs
    cost = zeros(Float64, length(features1), length(features2))
    for i in 1:length(features1)
        for j in 1:length(features2)
            cost[i, j] = sqrt(sum((features1[i] .- features2[j]).^2))
        end
    end

    distance = emd(weights1, weights2, cost)
    @show distance
    @test isapprox(distance, 160.542; atol = 1e-3)
end


text:
# Quasistatic analysis for beam, including several timesteps

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Postprocess
using JuliaFEM.Abaqus: create_surface_elements

Logging.configure(level=Logging.DEBUG)

function create_body(mesh, name, E, nu, rho)
    body = Problem(mesh, Elasticity, name, 3)
    update!(body, "youngs modulus", E)
    update!(body, "poissons ratio", nu)
    update!(body, "density", rho)
    return body
end

""" Create boundary condition from surface set. """
function create_bc_from_surface_set(mesh, name, u1, u2, u3)
    bc = Problem(Dirichlet, string(name), 3, "displacement")
    bc.elements = create_surface_elements(mesh, name)
    update!(bc, "displacement 1", u1)
    update!(bc, "displacement 2", u2)
    update!(bc, "displacement 3", u3)
    return bc
end

""" Create boundary condition from node set. """
function create_bc_from_node_set(mesh, name, u1, u2, u3)
    bc = Problem(Dirichlet, string(name), 3, "displacement")
    bc.elements = [Element(Poi1, [nid]) for nid in mesh.node_sets[name]]
    update!(bc, "geometry", mesh.nodes)
    update!(bc, "displacement 1", u1)
    update!(bc, "displacement 2", u2)
    update!(bc, "displacement 3", u3)
    return bc
end

function create_bc(mesh, name, u1=0.0, u2=0.0, u3=0.0)
    name = Symbol(name)
    if haskey(mesh.surface_sets, name)
        return create_bc_from_surface_set(mesh, name, u1, u2, u3)
    elseif haskey(mesh.node_sets, name)
        return create_bc_from_node_set(mesh, name, u1, u2, u3)
    else
        error("Mesh does not contain node or surface set $name")
    end
end

function create_interface(mesh, slave_surface::String, master_surface::String)
    interface = Problem(Mortar, "interface between $slave_surface and $master_surface", 3, "displacement")
    interface.properties.dual_basis = true
    slave_elements = create_surface_elements(mesh, Symbol(slave_surface))
    master_elements = create_surface_elements(mesh, Symbol(master_surface))
    nslaves = length(slave_elements)
    nmasters = length(master_elements)
    info("$nslaves slaves, $nmasters masters")
    update!(slave_elements, "master elements", master_elements)
    interface.elements = [slave_elements; master_elements]
    return interface
end

function create_interface(mesh, slave::Problem, master::Problem)
    slave_surface = slave.name * "_TO_" * master.name
    master_surface = master.name * "_TO_" * slave.name
    return create_interface(mesh, slave_surface, master_surface)
end

""" Convert Mesh object from quadratic to linear. """
function to_linear!(mesh)
    mapping = Dict(:Tet10 => :Tet4, :Tri6 => :Tri3)
    nnodes = Dict(:Tet4 => 4, :Tri3 => 3)
    for elid in keys(mesh.elements)
        eltype = mesh.element_types[elid]
        if haskey(mapping, eltype)
            mesh.element_types[elid] = mapping[eltype]
            nnodes_new = nnodes[mapping[eltype]]
            mesh.elements[elid] = mesh.elements[elid][1:nnodes_new]
        end
    end
end

# start of simulation
#mesh = abaqus_read_mesh("beam_3d_1st_order_tetra_30mm.inp")
mesh = abaqus_read_mesh("beam_3d_2nd_order_tetra_30mm.inp")
#to_linear!(mesh)
info("element sets = ", collect(keys(mesh.element_sets)))
info("surface sets = ", collect(keys(mesh.surface_sets)))

# parts
beam = create_body(mesh, "beam", 210.0e3, 0.3, 7.85e-9)

# boundary conditions
bc1 = create_bc(mesh, "fixed")

# load
load = Problem(Elasticity, "pressure load", 3)
load.elements = create_surface_elements(mesh, :load)
update!(load, "surface pressure", 0.0 => 0.0)
update!(load, "surface pressure", 1.0 => 20.0)

# give some information about surface load
nload = length(load.elements)
info("$nload elements in load surface")
area = 0.0
time = 0.0
for element in load.elements
    for ip in get_integration_points(element)
        detJ = element(ip, time, Val{:detJ})
        area += ip.weight*detJ
    end
end
info("load surface area: $area")

# solution
isfile("results_quasistatic.h5") && rm("results_quasistatic.h5")
isfile("results_quasistatic.xmf") && rm("results_quasistatic.xmf")
solver = Solver(Linear, beam, bc1, load)
solver.xdmf = Xdmf("results")

for time in [0.0, 1.0]
    solver.time = time
    solver()
end

# calculate stresses in integration points and use least-squares fitting to
# extrapolate results to nodes

""" Return stress tensor. """
function get_stress_tensor(element, ip, time)
    haskey(element, "displacement") || return nothing
    gradu = element("displacement", ip, time, Val{:Grad})
    eps = 0.5*(gradu' + gradu)
    E = element("youngs modulus", ip, time)
    nu = element("poissons ratio", ip, time)
    mu = E/(2.0*(1.0+nu))
    la = E*nu/((1.0+nu)*(1.0-2.0*nu))
    S = la*trace(eps)*I + 2.0*mu*eps
    return S
end

""" Return stress vector in "ABAQUS" order 11, 22, 33, 12, 23, 13. """
function get_stress(element, ip, time)
    S = get_stress_tensor(element, ip, time)
    return [S[1,1], S[2,2], S[3,3], S[1,2], S[2,3], S[1,3]]
end

""" Return principal stresses. """
function get_stress_principal(element, ip, time)
    S = get_stress_tensor(element, ip, time)
    return sort(eigvals(S))
end

""" Make least squares fit for some field to nodes. """
function lsq_fit(elements, field, time)
    A = SparseMatrixCOO()
    b = SparseMatrixCOO()
    volume = 0.0
    for element in elements
        gdofs = get_connectivity(element)
        # increase integration order by 1 from default
        for ip in get_integration_points(element, 1)
	        detJ = element(ip, time, Val{:detJ})
	        w = ip.weight*detJ
	        N = element(ip, time)
	        f = field(element, ip, time)
	        add!(A, gdofs, gdofs, w*kron(N', N))
	        for i=1:length(f)
                add!(b, gdofs, w*f[i]*N, i)
            end
            volume += w
	    end
    end
    info("Mass matrix for least-squares fit is assembled. Total volume to fit: $volume")
    A = sparse(A)
    b = sparse(b)
    A = 1/2*(A + A')

    SparseArrays.droptol!(A, 1.0e-6)
    SparseArrays.dropzeros!(A)
    nz = get_nonzero_rows(A)
    F = ldltfact(A[nz,nz])

    x = zeros(size(b)...)
    x[nz, :] = F \ b[nz, :]

    nodal_values = Dict(i => vec(x[i,:]) for i in nz)
    return nodal_values
end

xdmf = get(solver.xdmf)

using JuliaFEM: new_dataitem, new_child, set_attribute, add_child, save!

""" Store values to xml/h5. """
function add_field_to_xdmf!(xdmf, field_name, xmlpath, xdmfpath, data::Dict; field_type="Vector")
    node_ids = sort(collect(keys(data)))
    info("Storing field $field_name to xdmfpath $xdmfpath, xmlpath $xmlpath")
    for (i, nid) in enumerate(node_ids)
        println("$nid -> $(data[nid])")
        if i > 2
            println("...")
            break
        end
    end
    datavec = hcat([data[nid] for nid in node_ids]...)
    dataitem = new_dataitem(xdmf, xdmfpath, datavec)
    frame = read(xdmf, xmlpath)
    attribute = new_child(frame, "Attribute")
    set_attribute(attribute, "Name", field_name)
    set_attribute(attribute, "Center", "Node")
    set_attribute(attribute, "AttributeType", field_type)
    add_child(attribute, dataitem)
end

# calculate and update principle stresses to xdmf file, also calculate minimum
# and maximum principle components over time

minp = Dict()
maxp = Dict()

for (i, time) in enumerate([0.0, 1.0])
    S = lsq_fit(beam.elements, get_stress, time)
    #Sp = lsq_fit(beam.elements, get_stress_principal)

    # Calculate principal stresses in nodes
    Sp = Dict()
    for (nid, s) in S
        # order is: 11, 22, 33, 12, 23, 13
        stress_tensor = [
            s[1] s[4] s[6]
            s[4] s[2] s[5]
            s[6] s[5] s[3]]
        principle = sort(eigvals(stress_tensor))
        Sp[nid] = principle

        # minimum and maximum of principle stress over time
        if haskey(minp, nid)
            minp[nid] = min(minp[nid], minimum(principle))
        else
            minp[nid] = minimum(principle)
        end
        if haskey(maxp, nid)
            maxp[nid] = max(maxp[nid], maximum(principle))
        else
            maxp[nid] = maximum(principle)
        end
    end

    xml_path = "/Domain/Grid/Grid[$i]"
    xdmf_path = "/Results/Time $time/Nodal Fields/Stress"
    add_field_to_xdmf!(xdmf, "Stress", xml_path, xdmf_path, S; field_type="Tensor6")
    xdmf_path = "/Results/Time $time/Nodal Fields/Principal Stress"
    add_field_to_xdmf!(xdmf, "Principal Stress", xml_path, xdmf_path, Sp; field_type="Vector")
end

save!(xdmf)
close(xdmf.hdf)

# calculate mean stress and stress amplitude

Sa = Dict()
Sm = Dict()
node_ids = keys(minp)
for nid in node_ids
    Sa[nid] = (maxp[nid] - minp[nid]) / 2
    Sm[nid] = (maxp[nid] + minp[nid]) / 2
end

# create a new xdmf result output and save mean stress and stress amplitude there

using JuliaFEM: new_element, get_all_elements, get_element_type,
                filter_by_element_type, get_element_id

isfile("results_cycle.h5") && rm("results_cycle.h5")
isfile("results_cycle.xmf") && rm("results_cycle.xmf")
xdmf = Xdmf("results_cycle")
domain = new_child(xdmf, "Domain")
frame = new_child(domain, "Grid")

# 1. save geometry
X_ = solver("geometry", solver.time)
node_ids = sort(collect(keys(X_)))
X = hcat([X_[nid] for nid in node_ids]...)
data_node_ids = new_dataitem(xdmf, "/Node IDs", node_ids)
geometry = new_element("Geometry", Dict("Type" => "XYZ"))
add_child(geometry, new_dataitem(xdmf, "/Geometry", X))
add_child(frame, geometry)

# 2. save topology
nid_mapping = Dict(j=>i for (i, j) in enumerate(node_ids))
all_elements = get_all_elements(solver)
nelements = length(all_elements)
debug("Saving topology: $nelements elements total.")
element_types = unique(map(get_element_type, all_elements))

xdmf_element_mapping = Dict(
    "Poi1" => "Polyvertex",
    "Seg2" => "Polyline",
    "Tri3" => "Triangle",
    "Quad4" => "Quadrilateral",
    "Tet4" => "Tetrahedron",
    "Pyramid5" => "Pyramid",
    "Wedge6" => "Wedge",
    "Hex8" => "Hexahedron",
    "Seg3" => "Edge_3",
    "Tri6" => "Tri_6",
    "Quad8" => "Quad_8",
    "Tet10" => "Tet_10",
    "Pyramid13" => "Pyramid_13",
    "Wedge15" => "Wedge_15",
    "Hex20" => "Hex_20")

for element_type in element_types
    elements = filter_by_element_type(element_type, all_elements)
    nelements = length(elements)
    info("Xdmf save: $nelements elements of type $element_type")
    sort!(elements, by=get_element_id)
    element_ids = map(get_element_id, elements)
    element_conn = map(element -> [nid_mapping[j]-1 for j in get_connectivity(element)], elements)
    element_conn = hcat(element_conn...)
    element_code = split(string(element_type), ".")[end]
    dataitem = new_dataitem(xdmf, "/Topology/$element_code/Element IDs", element_ids)
    dataitem = new_dataitem(xdmf, "/Topology/$element_code/Connectivity", element_conn)
    topology = new_element("Topology")
    set_attribute(topology, "TopologyType", xdmf_element_mapping[element_code])
    set_attribute(topology, "NumberOfElements", length(elements))
    add_child(topology, dataitem)
    add_child(frame, topology)
end

add_field_to_xdmf!(xdmf, "Stress (Amplitude)", "/Domain/Grid", "/Results/Stress (Amplitude)", Sa; field_type="Scalar")
add_field_to_xdmf!(xdmf, "Stress (Mean)", "/Domain/Grid", "/Results/Stress (Mean", Sm; field_type="Scalar")

save!(xdmf)
close(xdmf.hdf)