text:
export Thunk, delayed, delayedmap

let counter=0
    global next_id
    next_id() = counter+=1
end

global _thunk_dict = Dict{Int, Any}()

# A thing to run
mutable struct Thunk
    f::Function
    inputs::Tuple
    id::Int
    get_result::Bool # whether the worker should send the result or only the metadata
    meta::Bool
    persist::Bool # don't `free!` result after computing
    cache::Bool   # release the result giving the worker an opportunity to
                  # cache it
    cache_ref::Nullable
    affinity::Nullable{Vector{Pair{Processor, Int}}}
    function Thunk(f, xs...;
                   id::Int=next_id(),
                   get_result::Bool=false,
                   meta::Bool=false,
                   persist::Bool=false,
                   cache::Bool=false,
                   cache_ref::Nullable{Any}=Nullable{Any}(),
                   affinity=Nullable(),
                  )
        thunk = new(f,xs,id,get_result,meta,persist, cache, cache_ref, affinity)
        _thunk_dict[id] = thunk
        thunk
    end
end

function affinity(t::Thunk)
    if !isnull(t.affinity)
        @logmsg("$t has (cached) affinity: $(get(t.affinity))")
        return get(t.affinity)
    end

    if t.cache && !isnull(t.cache_ref)
        aff_vec = affinity(get(t.cache_ref))
    else
        aff = Dict{Processor,Int}()
        for inp in inputs(t)
           #if haskey(state.cache, inp)
           #    as = affinity(state.cache[inp])
           #    for a in as
           #        proc, sz = a
           #        aff[proc] = get(aff, proc, 0) + sz
           #    end
           #else
                if isa(inp, Union{Chunk, Thunk})
                    for a in affinity(inp)
                        proc, sz = a
                        aff[proc] = get(aff, proc, 0) + sz
                    end
                end
           #end
        end
        aff_vec = collect(aff)
    end
    @logmsg("$t has affinity: $aff_vec")
    aff_vec
   #if length(aff) > 1
   #    return sort!(aff_vec, by=last,rev=true)
   #else
   #    return aff_vec
   #end
end

function delayed(f; kwargs...)
    (args...) -> Thunk(f, args...; kwargs...)
end

delayedmap(f, xs...) = map(delayed(f), xs...)

persist!(t::Thunk) = (t.persist=true; t)
cache_result!(t::Thunk) = (t.cache=true; t)

# @generated function compose{N}(f, g, t::NTuple{N})
#     if N <= 4
#       ( :(()->f(g())),
#         :((a)->f(g(a))),
#         :((a,b)->f(g(a,b))),
#         :((a,b,c)->f(g(a,b,c))),
#         :((a,b,c,d)->f(g(a,b,c,d))), )[N+1]
#     else
#         :((xs...) -> f(g(xs...)))
#     end
# end

# function Thunk(f::Function, t::Tuple{Thunk})
#     g = compose(f, t[1].f, t[1].inputs)
#     @logmsg(string("FUSING ", f, "*", t[1].f))
#     Thunk(g, t[1].inputs)
# end

# this gives a ~30x speedup in hashing
Base.hash(x::Thunk, h::UInt) = hash(x.id, hash(h, 0x7ad3bac49089a05f))
Base.isequal(x::Thunk, y::Thunk) = x.id==y.id

function Base.show(io::IO, p::Thunk)
    write(io, "*$(p.id)*")
end

inputs(x::Thunk) = x.inputs
inputs(x) = ()

istask(x::Thunk) = true
istask(x) = false


text:
include("HelmholtzForwardTest.jl");
include("ShiftedLaplacianTest.jl");
include("AttenuationTest.jl");
println("==== fictitious source test (2D) ===")
include("testFictitiousSource2D.jl")

ext: jl


text:
module DecodeTest

using Test
using ErlPort.Exceptions
using ErlPort.ErlTerms.Decode

function run()::Nothing

    # XXX Do we need to test REFERENCE_EXT, PORT_EXT, PID_EXT, NEW_REFERENCE_EXT, FUN_EXT, NEW_FUN_EXT, EXPORT_EXT

    # data for use by tests
    badsizecompdata = vcat(b"\x83P\0\0\0\x16\x78\xda\xcb\x66\x10\x49\xc1\2\0\x5d",
                           b"\x60\x08\x50")
    compdata1 = b"\x83P\0\0\0\x17\x78\xda\xcb\x66\x10\x49\xc1\2\0\x5d\x60\x08\x50"
    compdata2 = vcat(b"\x83P\0\0\0\x17\x78\xda\xcb\x66\x10\x49\xc1\2\0\x5d\x60",
                     b"\x08\x50tail")

    tuplepart = b"d\0\1ad\0\1bd\0\1cd\0\1dd\0\1ed\0\1fd\0\1gd\0\1hd\0\1id\0\1jd\0\1kd\0\1ld\0\1md\0\1nd\0\1od\0\1p"
    largetuple1 = vcat(b"i\0\0\1\0",
                       tuplepart, tuplepart, tuplepart, tuplepart, tuplepart, tuplepart,
                       tuplepart, tuplepart, tuplepart, tuplepart, tuplepart, tuplepart,
                       tuplepart, tuplepart, tuplepart, tuplepart)

    smallbigintmax = vcat(b"n\xff\0", fill(UInt8(255), 255))
    largebigintmin = vcat(b"o\0\0\1\0\0", fill(UInt8(0), 255), b"\1")

    # tests for supporting functions
    @testset "tests for supporting functions" begin
        @test size1unpack(b"\x0a")             == 0x0a
        @test size1unpack(b"\xaa")             == 0xaa
        @test size2unpack(b"\x0a\x0b")         == 0x0a0b
        @test size2unpack(b"\xaa\xbb")         == 0xaabb
        @test size4unpack(b"\x0a\x0b\x0c\x0d") == 0x0a0b0c0d
        @test size4unpack(b"\xaa\xbb\xcc\xdd") == 0xaabbccdd

        expected = [107,0,20,100,100,100,100,100,100,100,100,100,100,
                    100,100,100,100,100,100,100,100,100,100]
        @test_throws NotImplemented decompressterm(compdata1)
    end

    # basic decode errors
    @testset "basic decode errors" begin
        @test_throws UnsupportedType decode("")
        @test_throws IncompleteData decode(b"")
        @test_throws UnknownProtocolVersion decode(b"\0")
        @test_throws IncompleteData decode(b"\x83")
    end

    # decode empty list (that is in fact a nil) (NIL_EXT)
    @testset "decode empty list" begin
        @test_throws IncompleteData decodenil(b"")
        #@test_throws IncompleteData decodenil(b"\0") # XXX decode functions don't check term tag
        @test decode_with_tail(b"\x83j") == ([], b"")
        @test decode(b"\x83j") == []
        @test decodenil(b"j") == ([], b"")
        @test decodenil(b"jtail") == ([], b"tail")
    end

    # decode atoms (ATOM_EXT)
    @testset "decode atoms" begin
        @test_throws IncompleteData decode(b"\x83d")
        @test_throws IncompleteData decode(b"\x83d\0")
        @test_throws IncompleteData decode(b"\x83d\0\1")
        @test decode(b"\x83d\0\0") == Symbol("")
        @test_throws IncompleteData decode(b"\x83d\0\0tail")
        @test decode(b"\x83d\0\4test") == :test
        @test_throws IncompleteData decode(b"\x83d\0\4testtail")
        @test decode(b"\x83d\0\4true") == true
        @test decode(b"\x83d\0\5false") == false
        @test decode(b"\x83d\0\x09undefined") == nothing

        @test decode_with_tail(b"\x83d\0\0") == (Symbol(""), [])
        @test decode_with_tail(b"\x83d\0\0tail") == (Symbol(""), b"tail")
        @test decode_with_tail(b"\x83d\0\4test") == (:test, b"")
        @test decode_with_tail(b"\x83d\0\4testtail") == (:test, b"tail")
        @test decode_with_tail(b"\x83d\0\4true") == (true, b"")
        @test decode_with_tail(b"\x83d\0\5false") == (false, b"")
        @test decode_with_tail(b"\x83d\0\x09undefined") == (nothing, b"")

        @test decodeatom(b"d\0\0") == (Symbol(""), [])
        @test decodeatom(b"d\0\0tail") == (Symbol(""), b"tail")
        @test decodeatom(b"d\0\4test") == (:test, b"")
        @test decodeatom(b"d\0\4testtail") == (:test, b"tail")
        @test decodeatom(b"d\0\4true") == (true, b"")
        @test decodeatom(b"d\0\5false") == (false, b"")
        @test decodeatom(b"d\0\x09undefined") == (nothing, b"")
    end

    # decode predefined atoms (ATOM_CACHE_REF)
    @testset "decode predefined atoms" begin
        #atomcachereftag = UInt8(82)
    end

    # decode small atoms (SMALL_ATOM_EXT)
    @testset "decode small atoms" begin
        #smallatomtag = UInt8(115)
    end

    # decode UTF-8 atoms (ATOM_UTF8_EXT)
    @testset "decode UTF-8 atoms" begin
        #atomutf8tag = UInt8(118)
    end

    # decode small, UTF-8 atoms (SMALL_ATOM_UTF8_EXT)
    @testset "decode small, UTF-8 atoms" begin
        #smallatomutf8tag = UInt8(119)
    end

    # decode string list (STRING_EXT)
    @testset "decode string list" begin
        @test_throws IncompleteData decodestring(b"")
        @test_throws IncompleteData decodestring(b"\0")
        @test_throws IncompleteData decodestring(b"\0\0")

        @test decode(b"\x83k\0\1\1") == b"\1"
        @test decode(b"\x83k\0\6string") == b"string"

        @test decode_with_tail(b"\x83k\0\1\1") == (b"\1", b"")
        @test decode_with_tail(b"\x83k\0\6string") == (b"string", b"")

        @test decodestring(b"k\0\0") == (b"", b"")
        @test decodestring(b"k\0\0tail") == (b"", b"tail")
        @test decodestring(b"k\0\6string") == (b"string", b"")
        @test decodestring(b"k\0\x1f\xc3\xa1rv\xc3\xadzt\xc5\xb1r\xc5\x91\x20t\xc3\xbck\xc3\xb6rf\xc3\xbar\xc3\xb3g\xc3\xa9p") == (b"árvíztűrő tükörfúrógép", b"")
    end

    # decode list (LIST_EXT)
        @testset "decode list" begin
        @test decode(b"\x83j") == []
        @test decode(b"\x83l\0\0\0\4a\1d\0\1aa\3d\0\x09undefinedj") == [1,:a,3,nothing]
        @test decode(b"\x83\x6c\0\0\0\0\x6a") == []

        @test decode_with_tail(b"\x83j") == ([], b"")
        @test decode_with_tail(b"\x83l\0\0\0\4a\1d\0\1aa\3d\0\x09undefinedj") == ([1,:a,3,nothing], b"")

        @test decodelist(b"l\0\0\0\4a\1d\0\1aa\3d\0\x09undefinedj") == ([1,:a,3,nothing], b"")
        @test decodelist(b"l\0\0\0\4a\1d\0\1aa\3d\0\x09undefinedjtail") == ([1,:a,3,nothing], b"tail")
    end

    # decode improper list (LIST_EXT continued)
    # XXX how should we handle improper lists in Julia?

    # decode small tuple (SMALL_TUPLE_EXT)
    @testset "decode small tuple" begin
        @test_throws IncompleteData decodesmalltuple(b"")
        @test_throws IncompleteData decodesmalltuple(b"\0")
        @test_throws IncompleteData decodesmalltuple(b"h\1")
        @test decode(b"\x83h\2a\1a\2") == (1,2)
        @test decode(b"\x83h\3a\x0aa\x14a\x1e") == (10,20,30)
        @test decode(b"\x83h\4d\0\1ad\0\1bd\0\1cd\0\1d") == (:a,:b,:c,:d)
        @test_throws IncompleteData decode(b"\x83h\2a\1a\2tail")
        @test decode_with_tail(b"\x83h\2a\1a\2") == ((1,2), b"")
        @test decode_with_tail(b"\x83h\3a\x0aa\x14a\x1e") == ((10,20,30), b"")
        @test decode_with_tail(b"\x83h\4d\0\1ad\0\1bd\0\1cd\0\1d") == ((:a,:b,:c,:d), b"")
        @test decode_with_tail(b"\x83h\2a\1a\2tail") == ((1,2), b"tail")
        @test decodesmalltuple(b"h\0") == ((), b"")
        @test decodesmalltuple(b"h\0tail") == ((), b"tail")
        # it's a real tuple that has only 1 element
        @test decodesmalltuple(b"h\1a\1") == ((1,), b"")
        @test decodesmalltuple(b"h\2a\1a\2tail") == ((1,2), b"tail")
    end

    # decode large tuple (LARGE_TUPLE_EXT)
    @testset "decode large tuple" begin
        (lt1, tail1) = decode_with_tail(vcat(b"\x83", largetuple1, b"tail"))
        @test tail1 == b"tail"

        (lt2, tail2) = decodelargetuple(largetuple1)
        @test (length(lt2), tail2) == (256, b"")
        @test lt2[1]   == :a
        @test lt2[256] == :p
    end

    # decode opaque object

    # decode small integer (SMALL_INTEGER_EXT)
    @testset "decode small integer" begin
        @test_throws IncompleteData decode(b"\x83a")
        @test decode(b"\x83a\0") == 0
        @test_throws IncompleteData decode(b"\x83a\0tail")
        @test decode(b"\x83a\xff") == 255
        @test decode_with_tail(b"\x83a\0") == (0, b"")
        @test decode_with_tail(b"\x83a\0tail") == (0, b"tail")
        @test decode_with_tail(b"\x83a\xff") == (255, b"")
        @test decode_with_tail(b"\x83a\xfftail") == (255, b"tail")
        @test_throws IncompleteData decodesmallint(b"a")
        @test decodesmallint(b"a\0") == (0, b"")
        @test decodesmallint(b"a\0tail") == (0, b"tail")
        @test decodesmallint(b"a\xff") == (255, b"")
        @test decodesmallint(b"a\xfftail") == (255, b"tail")
    end

    # decode integer (INTEGER_EXT)
    @testset "decode integer" begin
        @test_throws IncompleteData decode(b"\x83b")
        @test_throws IncompleteData decode(b"\x83b\0")
        @test_throws IncompleteData decode(b"\x83b\0\0")
        @test_throws IncompleteData decode(b"\x83b\0\0\0")
        @test decode(b"\x83b\0\0\0\0") == 0
        @test_throws IncompleteData decode(b"\x83b\0\0\0\0tail")
        @test decode(b"\x83b\x7f\xff\xff\xff") == 2147483647
        @test_throws IncompleteData decode(b"\x83b\x7f\xff\xff\xfftail")
        @test decode(b"\x83b\xff\xff\xff\xff") == -1
        @test_throws IncompleteData decode(b"\x83b\xff\xff\xff\xfftail")

        @test decode_with_tail(b"\x83b\0\0\0\0") == (0, [])
        @test decode_with_tail(b"\x83b\0\0\0\0tail") == (0, b"tail")
        @test decode_with_tail(b"\x83b\x7f\xff\xff\xff") == (2147483647, [])
        @test decode_with_tail(b"\x83b\x7f\xff\xff\xfftail") == (2147483647, b"tail")
        @test decode_with_tail(b"\x83b\xff\xff\xff\xff") == (-1, [])
        @test decode_with_tail(b"\x83b\xff\xff\xff\xfftail") == (-1, b"tail")

        @test_throws IncompleteData decodeint(b"b")
        @test_throws IncompleteData decodeint(b"b\0")
        @test_throws IncompleteData decodeint(b"b\0\0")
        @test_throws IncompleteData decodeint(b"b\0\0\0")
        @test decodeint(b"b\0\0\0\0") == (0, [])
        @test decodeint(b"b\0\0\0\0tail") == (0, b"tail")
        @test decodeint(b"b\x7f\xff\xff\xff") == (2147483647, [])
        @test decodeint(b"b\x7f\xff\xff\xfftail") == (2147483647, b"tail")
        @test decodeint(b"b\xff\xff\xff\xff") == (-1, [])
        @test decodeint(b"b\xff\xff\xff\xfftail") == (-1, b"tail")
    end

    # decode binary (BINARY_EXT)
    @testset "decode binary" begin
        @test_throws IncompleteData decode(b"\x83m")
        @test_throws IncompleteData decode(b"\x83m\0")
        @test_throws IncompleteData decode(b"\x83m\0\0")
        @test_throws IncompleteData decode(b"\x83m\0\0\0")
        @test_throws IncompleteData decode(b"\x83m\0\0\0\1")
        @test decode(b"\x83m\0\0\0\0") == b""
        @test_throws IncompleteData decode(b"\x83m\0\0\0\0tail")
        @test decode(b"\x83m\0\0\0\4data") == b"data"
        @test_throws IncompleteData decode(b"\x83m\0\0\0\4datatail")

        @test decode_with_tail(b"\x83m\0\0\0\0") == (b"", b"")
        @test decode_with_tail(b"\x83m\0\0\0\0tail") == (b"", b"tail")
        @test decode_with_tail(b"\x83m\0\0\0\4data") == (b"data", b"")
        @test decode_with_tail(b"\x83m\0\0\0\4datatail") == (b"data", b"tail")

        @test decodebin(b"m\0\0\0\0") == (b"", b"")
        @test decodebin(b"m\0\0\0\0tail") == (b"", b"tail")
        @test decodebin(b"m\0\0\0\4data") == (b"data", b"")
        @test decodebin(b"m\0\0\0\4datatail") == (b"data", b"tail")
        @test_throws IncompleteData decodebin(b"m")
        @test_throws IncompleteData decodebin(b"m\0")
        @test_throws IncompleteData decodebin(b"m\0\0")
        @test_throws IncompleteData decodebin(b"m\0\0\0")
        @test_throws IncompleteData decodebin(b"m\0\0\0\1")
        @test decodebin(b"m\0\0\0\x1f\xc3\xa1rv\xc3\xadzt\xc5\xb1r\xc5\x91\x20t\xc3\xbck\xc3\xb6rf\xc3\xbar\xc3\xb3g\xc3\xa9p") == (b"árvíztűrő tükörfúrógép", b"")
    end

    # decode bitstring (BIT_BINARY_EXT)
    @testset "decode bitstring" begin
        # XXX follow-up
        #@test_throws IncompleteData decode(b"\x83M")
        #@test_throws IncompleteData decode(b"\x83M\0")
        #@test_throws IncompleteData decode(b"\x83M\0\0")
        #@test_throws IncompleteData decode(b"\x83M\0\0\0")
        #@test_throws IncompleteData decode(b"\x83M\0\0\0\0")
        #@test_throws IncompleteData decode(b"\x83M\0\0\0\1\0")
        #@test_throws IncompleteData decode(b"\x83M\0\0\0\1\0")
        #@test decode(b"\x83M\0\0\0\0\0") == (b"", b"")
        #@test decode(b"\x83M\0\0\0\0\3") == (b"", b"")
        #@test decode(b"\x83M\0\0\0\0\5tail") == (b"", b"tail")
        #@test decode(b"\x83M\0\0\0\1\1\x80") == (b"\1", b"")
        #@test decode(b"\x83M\0\0\0\1\1\x80tail") == (b"\1", b"tail")
    end

    # decode new float (NEW_FLOAT_EXT)
    @testset "decode float" begin
        @test_throws IncompleteData decode(b"\x83\x63")
        @test_throws IncompleteData decode(b"\x83\x63\0")
        @test_throws IncompleteData decode(b"\x83\x63\0\0")
        @test_throws IncompleteData decode(b"\x83\x63\0\0\0")
        @test_throws IncompleteData decode(b"\x83\x63\0\0\0\0")
        @test_throws IncompleteData decode(b"\x83\x63\0\0\0\0\0")
        @test_throws IncompleteData decode(b"\x83\x63\0\0\0\0\0\0")
        @test_throws IncompleteData decode(b"\x83\x63\0\0\0\0\0\0\0")
        @test decode(b"\x83\x63\0\0\0\0\0\0\0\0") == 0.0
        @test_throws IncompleteData decode(b"\x83\x63\0\0\0\0\0\0\0\0tail")
        @test decode(b"\x83\x63\x3f\xf8\0\0\0\0\0\0") == 1.5
        @test_throws IncompleteData decode(b"\x83\x63\x3f\xf8\0\0\0\0\0\0tail")

        @test decode_with_tail(b"\x83\x63\0\0\0\0\0\0\0\0") == (0.0, b"")
        @test decode_with_tail(b"\x83\x63\0\0\0\0\0\0\0\0tail") == (0.0, b"tail")
        @test decode_with_tail(b"\x83\x63\x3f\xf8\0\0\0\0\0\0") == (1.5, b"")
        @test decode_with_tail(b"\x83\x63\x3f\xf8\0\0\0\0\0\0tail") == (1.5, b"tail")

        @test decodefloat(b"\x63\0\0\0\0\0\0\0\0") == (0.0, b"")
        @test decodefloat(b"\x63\0\0\0\0\0\0\0\0tail") == (0.0, b"tail")
        @test decodefloat(b"\x63\x3f\xf8\0\0\0\0\0\0") == (1.5, b"")
        @test decodefloat(b"\x63\x3f\xf8\0\0\0\0\0\0tail") == (1.5, b"tail")
        @test_throws IncompleteData decodefloat(b"\x63")
        @test_throws IncompleteData decodefloat(b"\x63\0")
        @test_throws IncompleteData decodefloat(b"\x63\0\0")
        @test_throws IncompleteData decodefloat(b"\x63\0\0\0")
        @test_throws IncompleteData decodefloat(b"\x63\0\0\0\0")
        @test_throws IncompleteData decodefloat(b"\x63\0\0\0\0\0")
        @test_throws IncompleteData decodefloat(b"\x63\0\0\0\0\0\0")
        @test_throws IncompleteData decodefloat(b"\x63\0\0\0\0\0\0\0")
    end

    # decode new float (NEW_FLOAT_EXT)
    @testset "decode new float" begin
        @test_throws IncompleteData decode(b"\x83F")
        @test_throws IncompleteData decode(b"\x83F\0")
        @test_throws IncompleteData decode(b"\x83F\0\0")
        @test_throws IncompleteData decode(b"\x83F\0\0\0")
        @test_throws IncompleteData decode(b"\x83F\0\0\0\0")
        @test_throws IncompleteData decode(b"\x83F\0\0\0\0\0")
        @test_throws IncompleteData decode(b"\x83F\0\0\0\0\0\0")
        @test_throws IncompleteData decode(b"\x83F\0\0\0\0\0\0\0")
        @test decode(b"\x83F\0\0\0\0\0\0\0\0") == 0.0
        @test_throws IncompleteData decode(b"\x83F\0\0\0\0\0\0\0\0tail")
        @test decode(b"\x83F?\xf8\0\0\0\0\0\0") == 1.5
        @test_throws IncompleteData decode(b"\x83F?\xf8\0\0\0\0\0\0tail")
        @test isnan(decode(b"\x83d\0\3nan"))

        @test decode_with_tail(b"\x83F\0\0\0\0\0\0\0\0") == (0.0, b"")
        @test decode_with_tail(b"\x83F\0\0\0\0\0\0\0\0tail") == (0.0, b"tail")
        @test decode_with_tail(b"\x83F?\xf8\0\0\0\0\0\0") == (1.5, b"")
        @test decode_with_tail(b"\x83F?\xf8\0\0\0\0\0\0tail") == (1.5, b"tail")

        @test decodenewfloat(b"F\0\0\0\0\0\0\0\0") == (0.0, b"")
        @test decodenewfloat(b"F\0\0\0\0\0\0\0\0tail") == (0.0, b"tail")
        @test decodenewfloat(b"F?\xf8\0\0\0\0\0\0") == (1.5, b"")
        @test decodenewfloat(b"F?\xf8\0\0\0\0\0\0tail") == (1.5, b"tail")
        @test_throws IncompleteData decodenewfloat(b"F")
        @test_throws IncompleteData decodenewfloat(b"F\0")
        @test_throws IncompleteData decodenewfloat(b"F\0\0")
        @test_throws IncompleteData decodenewfloat(b"F\0\0\0")
        @test_throws IncompleteData decodenewfloat(b"F\0\0\0\0")
        @test_throws IncompleteData decodenewfloat(b"F\0\0\0\0\0")
        @test_throws IncompleteData decodenewfloat(b"F\0\0\0\0\0\0")
        @test_throws IncompleteData decodenewfloat(b"F\0\0\0\0\0\0\0")
    end

    # decode small big integer (SMALL_BIG_EXT)
    @testset "decode small big integer" begin
        @test_throws IncompleteData decode(b"\x83n")
        @test_throws IncompleteData decode(b"\x83n\0")
        @test_throws IncompleteData decode(b"\x83n\1\0")
        @test_throws IncompleteData decode(b"\x83n\2\0\0")

        @test decode(b"\x83n\1\0\0") == 0
        @test_throws IncompleteData decode(b"\x83n\1\0\0tail")

        @test decode_with_tail(b"\x83n\1\0\0") == (0, b"")
        @test decode_with_tail(b"\x83n\1\0\0tail") == (0, b"tail")

        @test decodesmallbigint(b"n\1\0\1") == (1, b"")
        @test decodesmallbigint(b"n\1\1\1") == (-1, b"")
        @test decodesmallbigint(b"n\2\0\1\2") == (513, b"")
        # TODO ??? @test decodesmallbigint(smallbigintmax) == (256^255-1, b"")
    end

    # decode large big integer (LARGE_BIG_EXT)
    @testset "decode large big integer" begin
        @test_throws IncompleteData decode(b"\x83o")
        @test_throws IncompleteData decode(b"\x83o\0")
        @test_throws IncompleteData decode(b"\x83o\0\0")
        @test_throws IncompleteData decode(b"\x83o\0\0\0")
        @test_throws IncompleteData decode(b"\x83o\0\0\0\0")
        @test_throws IncompleteData decode(b"\x83o\0\0\0\1\0")
        @test_throws IncompleteData decode(b"\x83o\0\0\0\2\0\0")
        @test decode(b"\x83o\0\0\0\1\0\0") == 0
        @test_throws IncompleteData decode(b"\x83o\0\0\0\1\0\0tail")

        @test decode_with_tail(b"\x83o\0\0\0\1\0\0") == (0, b"")
        @test decode_with_tail(b"\x83o\0\0\0\1\0\0tail") == (0, b"tail")

        @test decodelargebigint(b"o\0\0\0\1\0\1") == (1, b"")
        @test decodelargebigint(b"o\0\0\0\1\1\1") == (-1, b"")
        @test decodelargebigint(b"o\0\0\0\2\1\1\2") == (-513, b"")
        # TODO ??? @test decodelargebigint(largebigintmin) == (256^255, b"")
    end

    # decode map (MAP_EXT)
    @testset "decode map" begin
        @test_throws IncompleteData decode(b"\x83t")
        @test_throws IncompleteData decode(b"\x83t\0")

        @test decode(b"\x83t\0\0\0\0") == Dict()
        @test decode(b"\x83t\0\0\0\1d\0\1aa\2") == Dict(:a => 2)
        @test decode(b"\x83t\0\0\0\1d\0\1at\0\0\0\1d\0\1aa\2") == Dict(:a => Dict(:a => 2))
        @test decode(b"\x83t\0\0\0\2d\0\1aa\2k\0\1bj") == Dict(:a => 2, b"b" => [])

        @test decode_with_tail(b"\x83t\0\0\0\0") == (Dict(), b"")
        @test decode_with_tail(b"\x83t\0\0\0\1d\0\1aa\2") == (Dict(:a => 2), b"")
        @test decode_with_tail(b"\x83t\0\0\0\1d\0\1at\0\0\0\1d\0\1aa\2") == (Dict(:a => Dict(:a => 2)), b"")
        @test decode_with_tail(b"\x83t\0\0\0\2d\0\1aa\2k\0\1bj") == (Dict(:a => 2, b"b" => []), b"")

        @test decodemap(b"t\0\0\0\0tail") == (Dict(), b"tail")
        @test decodemap(b"t\0\0\0\1d\0\1aa\2tail") == (Dict(:a => 2), b"tail")
        @test decodemap(b"t\0\0\0\2d\0\1aa\2k\0\1bjtail") == (Dict(:a => 2, b"b" => []), b"tail")
    end

    # decode compressed term
    @testset "decode compressed term" begin
        #@test_throws IncompleteData decode(b"\x83P")
        #@test_throws IncompleteData decode(b"\x83P\0")
        #@test_throws IncompleteData decode(b"\x83P\0\0")
        #@test_throws IncompleteData decode(b"\x83P\0\0\0")
        #@test_throws IncompleteData decode(b"\x83P\0\0\0\0")
        #@test_throws InvalidCompressedTag decode(badsizecompdata)
        #@test decode(compdata1) == ([100,100,100,100,100,100,100,100,100,100,
        #                             100,100,100,100, 100,100,100,100,100,100],
        #                            UInt8[])

        @test_throws NotImplemented decode(b"\x83P")
        @test_throws NotImplemented decode(b"\x83P\0")
        @test_throws NotImplemented decode(b"\x83P\0\0")
        @test_throws NotImplemented decode(b"\x83P\0\0\0")
        @test_throws NotImplemented decode(b"\x83P\0\0\0\0")
        @test_throws NotImplemented decode(badsizecompdata)
        @test_throws NotImplemented decode(compdata1)
        @test_throws NotImplemented decode(compdata2)
    end

    return nothing
end # run()

end # module


text:
"""
    Path{D,T}
Alias of `MutableLinkedList{NTuple{D,T}}`. Used to represent the path to be
taken by an agent in a `D` dimensional space.
"""
const Path{D,T} = MutableLinkedList{NTuple{D,T}}

struct AStar{D,P,M,T,C<:CostMetric{D}} <: GridPathfinder{D,P,M}
    agent_paths::Dict{Int,Path{D,T}}
    dims::NTuple{D,T}
    neighborhood::Vector{CartesianIndex{D}}
    admissibility::Float64
    walkmap::BitArray{D}
    cost_metric::C

    function AStar{D,P,M,T,C}(
        agent_paths::Dict,
        dims::NTuple{D,T},
        neighborhood::Vector{CartesianIndex{D}},
        admissibility::Float64,
        walkmap::BitArray{D},
        cost_metric::C,
    ) where {D,P,M,C,T}
        @assert all(dims .> 0) "Invalid pathfinder dimensions: $(dims)"
        T <: Integer && @assert size(walkmap) == dims "Walkmap must be same dimensions as grid"
        @assert admissibility >= 0 "Invalid value for admissibility: $admissibility ≱ 0"
        if cost_metric isa PenaltyMap{D}
            @assert size(cost_metric.pmap) == size(walkmap) "Penaltymap dimensions must be same as walkable map"
        elseif cost_metric isa DirectDistance{D}
            if M
                @assert length(cost_metric.direction_costs) >= D "DirectDistance direction_costs must have as many values as dimensions"
            else
                @assert length(cost_metric.direction_costs) >= 1 "DirectDistance direction_costs must have non-zero length"
            end
        end
        new(agent_paths, dims, neighborhood, admissibility, walkmap, cost_metric)
    end
end

"""
    Pathfinding.AStar(space; kwargs...)
Enables pathfinding for agents in the provided `space` (which can be a [`GridSpace`](@ref) or
[`ContinuousSpace`](@ref)) using the A* algorithm. This struct must be passed into any
pathfinding functions.

For [`ContinuousSpace`](@ref), a walkmap or instance of [`PenaltyMap`](@ref) must be provided
to specify the level of discretisation of the space.

## Keywords
- `diagonal_movement = true` specifies if movement can be to diagonal neighbors of a
  tile, or only orthogonal neighbors. Only available for [`GridSpace`](@ref)
- `admissibility = 0.0` allows the algorithm to aprroximate paths to speed up pathfinding.
  A value of `admissibility` allows paths with at most `(1+admissibility)` times the optimal
  length.
- `walkmap = trues(size(space))` specifies the (un)walkable positions of the
  space. If specified, it should be a `BitArray` of the same size as the corresponding
  `GridSpace`. By default, agents can walk anywhere in the space. An example usage can
  be found in [Maze Solver](@ref)
- `cost_metric = DirectDistance{D}()` is an instance of a cost metric and specifies the
  metric used to approximate the distance between any two points. See [`CostMetric`](@ref).
  An example usage can be found in [Mountain Runners](@ref).
"""
function AStar(
    dims::NTuple{D,T};
    periodic::Bool = false,
    diagonal_movement::Bool = true,
    admissibility::Float64 = 0.0,
    walkmap::BitArray{D} = trues(dims),
    cost_metric::CostMetric{D} = DirectDistance{D}(),
) where {D,T}
    neighborhood = diagonal_movement ? moore_neighborhood(D) : vonneumann_neighborhood(D)
    return AStar{D,periodic,diagonal_movement,T,typeof(cost_metric)}(
        Dict{Int,Path{D,T}}(),
        dims,
        neighborhood,
        admissibility,
        walkmap,
        cost_metric,
    )
end

AStar(
    space::GridSpace{D,periodic};
    diagonal_movement::Bool = true,
    admissibility::Float64 = 0.0,
    walkmap::BitArray{D} = trues(size(space.s)),
    cost_metric::CostMetric{D} = DirectDistance{D}(),
) where {D,periodic} =
    AStar(size(space); periodic, diagonal_movement, admissibility, walkmap, cost_metric)

function AStar(
    space::ContinuousSpace{D,periodic};
    walkmap::Union{BitArray{D},Nothing} = nothing,
    admissibility::Float64 = 0.0,
    cost_metric::CostMetric{D} = DirectDistance{D}(),
) where {D,periodic}
    @assert walkmap isa BitArray{D} || cost_metric isa PenaltyMap "Pathfinding in ContinuousSpace requires either walkmap to be specified or cost_metric to be a PenaltyMap"
    isnothing(walkmap) && (walkmap = BitArray(trues(size(cost_metric.pmap))))
    AStar(size(space); periodic, diagonal_movement = true, admissibility, walkmap, cost_metric)
end

moore_neighborhood(D) = [
    CartesianIndex(a)
    for a in Iterators.product([-1:1 for φ in 1:D]...) if a != Tuple(zeros(Int, D))
]

function vonneumann_neighborhood(D)
    hypercube = CartesianIndices((repeat([-1:1], D)...,))
    [β for β ∈ hypercube if LinearAlgebra.norm(β.I) == 1]
end

function Base.show(io::IO, pathfinder::AStar{D,P,M}) where {D,P,M}
    periodic = P ? "periodic, " : ""
    moore = M ? "diagonal, " : "orthogonal, "
    s =
        "A* in $(D) dimensions, $(periodic)$(moore)ϵ=$(pathfinder.admissibility), " *
        "metric=$(pathfinder.cost_metric)"
    print(io, s)
end

struct GridCell
    f::Int
    g::Int
    h::Int
end

GridCell(g::Int, h::Int, admissibility::Float64) =
    GridCell(round(Int, g + (1 + admissibility) * h), g, h)

GridCell() = GridCell(typemax(Int), typemax(Int), typemax(Int))

ordering(cell) = cell.f

"""
    find_path(pathfinder::AStar{D}, from::NTuple{D,Int}, to::NTuple{D,Int})
Calculate the shortest path from `from` to `to` using the A* algorithm.
If a path does not exist between the given positions, an empty linked list is returned.
"""
function find_path(pathfinder::AStar{D}, from::Dims{D}, to::Dims{D}) where {D}
    if !all(1 .<= from .<= size(pathfinder.walkmap)) ||
        !all(1 .<= to .<= size(pathfinder.walkmap)) ||
        !pathfinder.walkmap[from...] ||
        !pathfinder.walkmap[to...]
        return # nothing
    end
    parent = Dict{Dims{D},Dims{D}}()

    open_list = PriorityQueue{Dims{D},GridCell}(Base.By(ordering))
    closed_list = Set{Dims{D}}()

    enqueue!(
        open_list,
        from,
        GridCell(0, delta_cost(pathfinder, from, to), pathfinder.admissibility)
    )

    while !isempty(open_list)
        cur, cell = dequeue_pair!(open_list)
        cur == to && break
        push!(closed_list, cur)

        nbors = get_neighbors(cur, pathfinder)
        for nbor in Iterators.filter(n -> inbounds(n, pathfinder, closed_list), nbors)
            nbor_cell = haskey(open_list, nbor) ? open_list[nbor] : GridCell()
            new_g_cost = cell.g + delta_cost(pathfinder, cur, nbor)

            if new_g_cost < nbor_cell.g
                parent[nbor] = cur
                open_list[nbor] = GridCell(
                    new_g_cost,
                    delta_cost(pathfinder, nbor, to),
                    pathfinder.admissibility,
                )
            end
        end
    end

    agent_path = Path{D,Int64}()
    cur = to
    while true
        haskey(parent, cur) || break
        pushfirst!(agent_path, cur)
        cur = parent[cur]
    end
    cur == from || return # nothing
    return agent_path
end

@inline get_neighbors(cur, pathfinder::AStar{D,true}) where {D} =
    (mod1.(cur .+ β.I, size(pathfinder.walkmap)) for β in pathfinder.neighborhood)
@inline get_neighbors(cur, pathfinder::AStar{D,false}) where {D} =
    (cur .+ β.I for β in pathfinder.neighborhood)
@inline inbounds(n, pathfinder, closed) =
    all(1 .<= n .<= size(pathfinder.walkmap)) && pathfinder.walkmap[n...] && n ∉ closed

Base.isempty(id::Int, pathfinder::AStar) =
    !haskey(pathfinder.agent_paths, id) || isempty(pathfinder.agent_paths[id])

Agents.is_stationary(
    agent::A,
    pathfinder::AStar,
) where {A<:AbstractAgent} = isempty(agent.id, pathfinder)

"""
    Pathfinding.penaltymap(pathfinder)
Return the penalty map of a [`Pathfinding.AStar`](@ref) if the
[`Pathfinding.PenaltyMap`](@ref) metric is in use, `nothing` otherwise.

It is possible to mutate the map directly, for example
`Pathfinding.penaltymap(pathfinder)[15, 40] = 115`
or `Pathfinding.penaltymap(pathfinder) .= rand(50, 50)`. If this is mutated,
a new path needs to be planned using [`plan_route!`](@ref).
"""
function penaltymap(pathfinder::AStar)
    if pathfinder.cost_metric isa PenaltyMap
        return pathfinder.cost_metric.pmap
    else
        return nothing
    end
end

"""
    Pathfinding.kill_agent!(agent, model, pathfinder)
The same as `kill_agent!(agent, model)`, but also removes the agent's path data
from `pathfinder`.
"""
function Agents.kill_agent!(
    agent::A,
    model::ABM{S,A},
    pathfinder::AStar,
) where {S,A<:AbstractAgent}
    delete!(pathfinder.agent_paths, agent.id)
    delete!(model.agents, agent.id)
    Agents.remove_agent_from_space!(agent, model)
end


text:
const QNBlock = Pair{QN,Int64}

const QNBlocks = Vector{QNBlock}

qn(qnblock::QNBlock) = qnblock.first

# Get the dimension of the specified block
blockdim(qnblock::QNBlock) = qnblock.second

# Get the dimension of the specified block
blockdim(qnblocks::QNBlocks, b::Integer) = blockdim(qnblocks[b])
blockdim(qnblocks::QNBlocks, b::Block{1}) = blockdim(qnblocks[only(b)])

# Get the QN of the specified block
qn(qnblocks::QNBlocks, b::Integer) = qn(qnblocks[b])
qn(qnblocks::QNBlocks, b::Block{1}) = qn(qnblocks[only(b)])

nblocks(qnblocks::QNBlocks) = length(qnblocks)

function dim(qnblocks::QNBlocks)
  dimtot = 0
  for (_, blockdim) in qnblocks
    dimtot += blockdim
  end
  return dimtot
end

function -(qnb::QNBlock)
  return QNBlock(-qn(qnb), blockdim(qnb))
end

function (qn1::QNBlock + qn2::QNBlock)
  qn(qn1) != qn(qn2) && error("Cannot add qn blocks with different qns")
  return QNBlock(qn(qn1), blockdim(qn1) + blockdim(qn2))
end

function -(qns::QNBlocks)
  qns_new = copy(qns)
  for i in 1:length(qns_new)
    qns_new[i] = -qns_new[i]
  end
  return qns_new
end

"""
A QN Index is an Index with QN block storage instead of
just an integer dimension. The QN block storage is a
vector of pairs of QNs and block dimensions.
The total dimension of a QN Index is the sum of the
dimensions of the blocks of the Index.
"""
const QNIndex = Index{QNBlocks}

# Trait for the symmetry type (QN or not QN)
struct HasQNs <: SymmetryStyle end

symmetrystyle(::QNIndex) = HasQNs()
symmetrystyle(::HasQNs, ::HasQNs) = HasQNs()
symmetrystyle(::NonQN, ::NonQN) = NonQN()
symmetrystyle(::HasQNs, ::NonQN) = HasQNs()
symmetrystyle(::NonQN, ::HasQNs) = HasQNs()

hasqns(::QNIndex) = true

function have_same_qns(qnblocks::QNBlocks)
  qn1 = qn(qnblocks, 1)
  for n in 2:nblocks(qnblocks)
    !have_same_qns(qn1, qn(qnblocks, n)) && return false
  end
  return true
end

function have_same_mods(qnblocks::QNBlocks)
  qn1 = qn(qnblocks, 1)
  for n in 2:nblocks(qnblocks)
    !have_same_mods(qn1, qn(qnblocks, n)) && return false
  end
  return true
end

"""
    Index(qnblocks::Vector{Pair{QN, Int64}}; dir::Arrow = Out,
                                             tags = "",
                                             plev::Integer = 0)

Construct a QN Index from a Vector of pairs of QN and block
dimensions.

Note: in the future, this may enforce that all blocks have the
same QNs (which would allow for some optimizations, for example
when constructing random QN ITensors).

# Example
```
Index([QN("Sz", -1) => 1, QN("Sz", 1) => 1]; tags = "i")
```
"""
function Index(qnblocks::QNBlocks; dir::Arrow=Out, tags="", plev=0)
  # TODO: make this a debug check?
  #have_same_qns(qnblocks) || error("When creating a QN Index, the QN blocks must have the same QNs")
  #have_same_mods(qnblocks) || error("When creating a QN Index, the QN blocks must have the same mods")
  return Index(rand(index_id_rng(), IDType), qnblocks, dir, tags, plev)
end

"""
    Index(qnblocks::Vector{Pair{QN, Int64}}, tags; dir::Arrow = Out,
                                                   plev::Integer = 0)

Construct a QN Index from a Vector of pairs of QN and block
dimensions.

# Example
```
Index([QN("Sz", -1) => 1, QN("Sz", 1) => 1], "i"; dir = In)
```
"""
function Index(qnblocks::QNBlocks, tags; dir::Arrow=Out, plev::Integer=0)
  return Index(qnblocks; dir=dir, tags=tags, plev=plev)
end

"""
    Index(qnblocks::Pair{QN, Int64}...; dir::Arrow = Out,
                                        tags = "",
                                        plev::Integer = 0)

Construct a QN Index from a list of pairs of QN and block
dimensions.

# Example
```
Index(QN("Sz", -1) => 1, QN("Sz", 1) => 1; tags = "i")
```
"""
function Index(qnblocks::QNBlock...; dir::Arrow=Out, tags="", plev=0)
  return Index([qnblocks...]; dir=dir, tags=tags, plev=plev)
end

dim(i::QNIndex) = dim(space(i))

nblocks(i::QNIndex) = nblocks(space(i))
# Define to be 1 for non-QN Index
nblocks(i::Index) = 1

# Get the Block that the index value falls in
# For example:
# qns = [QN(0,2) => 2, QN(0,2) => 2]
# block(qns, 1) == Block(1)
# block(qns, 2) == Block(1)
# block(qns, 3) == Block(2)
# block(qns, 4) == Block(2)
function block(qns::QNBlocks, n::Int)
  tdim = 0
  for b in 1:nblocks(qns)
    tdim += blockdim(qns, Block(b))
    (n <= tdim) && return Block(b)
  end
  error("qn: QN Index value out of range")
  return Block(0)
end

function block(iv::Pair{<:Index})
  i = ind(iv)
  v = val(iv)
  return block(space(i), v)
end

# Get the QN of the block
qn(i::QNIndex, b::Block{1}) = qn(space(i), b)

qn(ib::Pair{<:Index,Block{1}}) = qn(first(ib), last(ib))

# XXX: deprecate the Integer version
qn(i::QNIndex, b::Integer) = qn(i, Block(b))

# Get the QN of the block the IndexVal lies in
qn(iv::Pair{<:Index}) = qn(ind(iv), block(iv))

flux(i::QNIndex, b::Block{1}) = dir(i) * qn(i, b)

flux(ib::Pair{<:Index,Block{1}}) = flux(first(ib), last(ib))

flux(iv::Pair{<:Index}) = flux(ind(iv), block(iv))

qnblocks(i::QNIndex) = space(i)

# XXX: deprecate the Integer version
blockdim(i::QNIndex, b::Block) = blockdim(space(i), b)
blockdim(i::QNIndex, b::Integer) = blockdim(i, Block(b))
function blockdim(i::Index, b::Union{Block,Integer})
  return error(
    "`blockdim(i::Index, b)` not currently defined for non-QN Index $i of type `$(typeof(i))`. In the future this may be defined for `b == Block(1)` or `b == 1` as `dim(i)` and error otherwise.",
  )
end

dim(i::QNIndex, b::Block) = blockdim(space(i), b)

NDTensors.eachblock(i::Index) = (Block(n) for n in 1:nblocks(i))

# Return the first block of the QNIndex with the flux q
function block(::typeof(first), ind::QNIndex, q::QN)
  for b in eachblock(ind)
    if flux(ind => b) == q
      return b
    end
  end
  error("No block found with QN equal to $q")
  return Block(0)
end

# Find the first block that matches the pattern f,
# for example `f(blockind) = qn(blockind) == target_qn`.
# `f` accepts a pair of `i => Block(n)` where `n`
# runs over `nblocks(i)`.
function findfirstblock(f, i::QNIndex)
  for b in ITensors.eachblock(i)
    if f(i => b)
      return b
    end
  end
  error("No block of Index $i matching the specified pattern.")
  return Block(0)
end

# XXX: call this simply `block` and return a Block{1}
# Deprecate this
"""
    qnblocknum(ind::QNIndex, q::QN)

Given a QNIndex `ind` and QN `q`, return the
number of the block (from 1,...,nblocks(ind))
of the QNIndex having QN equal to `q`. Assumes
all blocks of `ind` have a unique QN.
"""
function qnblocknum(ind::QNIndex, q::QN)
  for b in 1:nblocks(ind)
    if flux(ind => Block(b)) == q
      return b
    end
  end
  error("No block found with QN equal to $q")
  return 0
end

blockdim(ind::QNIndex, q::QN) = blockdim(ind, block(first, ind, q))

# XXX: deprecate in favor of blockdim
"""
    qnblockdim(ind::QNIndex, q::QN)

Given a QNIndex `ind` and QN `q`, return the
dimension of the block of the QNIndex having
QN equal to `q`. Assumes all blocks of `ind`
have a unique QN.
"""
qnblockdim(ind::QNIndex, q::QN) = blockdim(ind, qnblocknum(ind, q))

(dir::Arrow * qnb::QNBlock) = QNBlock(dir * qn(qnb), blockdim(qnb))

function (dir::Arrow * qn::QNBlocks)
  # XXX use:
  # dir .* qn
  qnR = copy(qn)
  for i in 1:nblocks(qnR)
    qnR[i] = dir * qnR[i]
  end
  return qnR
end

(qn1::QNBlock * qn2::QNBlock) = QNBlock(qn(qn1) + qn(qn2), blockdim(qn1) * blockdim(qn2))

# TODO: rename tensorproduct with ⊗ alias
function outer(qn1::QNBlocks, qn2::QNBlocks)
  qnR = ITensors.QNBlocks(undef, nblocks(qn1) * nblocks(qn2))
  for (i, t) in enumerate(Iterators.product(qn1, qn2))
    qnR[i] = prod(t)
  end
  return qnR
end

# TODO: rename tensorproduct with ⊗ alias
function outer(i1::QNIndex, i2::QNIndex; dir=nothing, tags="", plev::Integer=0)
  if isnothing(dir)
    if ITensors.dir(i1) == ITensors.dir(i2)
      dir = ITensors.dir(i1)
    else
      dir = Out
    end
  end
  newspace = dir * ((ITensors.dir(i1) * space(i1)) ⊗ (ITensors.dir(i2) * space(i2)))
  return Index(newspace; dir=dir, tags=tags, plev=plev)
end

# TODO: rename tensorproduct with ⊗ alias
function outer(i::QNIndex; dir=nothing, tags="", plev::Integer=0)
  if isnothing(dir)
    dir = ITensors.dir(i)
  end
  newspace = dir * (ITensors.dir(i) * space(i))
  return Index(newspace; dir=dir, tags=tags, plev=plev)
end

# TODO: add ⊕ alias
function directsum(
  i::Index{Vector{Pair{QN,Int}}}, j::Index{Vector{Pair{QN,Int}}}; tags="sum"
)
  dir(i) ≠ dir(j) && error(
    "To direct sum two indices, they must have the same direction. Trying to direct sum indices $i and $j.",
  )
  return Index(vcat(space(i), space(j)); dir=dir(i), tags=tags)
end

isless(qnb1::QNBlock, qnb2::QNBlock) = isless(qn(qnb1), qn(qnb2))

function permuteblocks(i::QNIndex, perm)
  qnblocks_perm = space(i)[perm]
  return replaceqns(i, qnblocks_perm)
end

function combineblocks(qns::QNBlocks)
  perm = sortperm(qns)
  qnsP = qns[perm]
  qnsC = [qnsP[1]]
  comb = Vector{Int}(undef, nblocks(qns))

  # Which block this is, after combining
  block_count = 1
  comb[1] = block_count
  for i in 2:nblocks(qnsP)
    if qn(qnsP[i]) == qn(qnsP[i - 1])
      qnsC[block_count] += qnsP[i]
    else
      push!(qnsC, qnsP[i])
      block_count += 1
    end
    comb[i] = block_count
  end
  return qnsC, perm, comb
end

function splitblocks(qns::QNBlocks)
  idim = dim(qns)
  split_qns = similar(qns, idim)
  for n in 1:idim
    b = block(qns, n)
    split_qns[n] = qn(qns, b) => 1
  end
  return split_qns
end

# Make a new Index with the specified qn blocks
replaceqns(i::QNIndex, qns::QNBlocks) = setspace(i, qns)

function setblockdim!(i::QNIndex, newdim::Integer, n::Integer)
  qns = space(i)
  qns[n] = qn(qns[n]) => newdim
  return i
end

function setblockqn!(i::QNIndex, newqn::QN, n::Integer)
  qns = space(i)
  qns[n] = newqn => blockdim(qns[n])
  return i
end

function deleteat!(i::QNIndex, pos)
  deleteat!(space(i), pos)
  return i
end

function resize!(i::QNIndex, n::Integer)
  resize!(space(i), n)
  return i
end

function combineblocks(i::QNIndex)
  qnsR, perm, comb = combineblocks(space(i))
  iR = replaceqns(i, qnsR)
  return iR, perm, comb
end

removeqns(i::QNIndex) = setdir(setspace(i, dim(i)), Neither)

function addqns(i::Index, qns::QNBlocks; dir::Arrow=Out)
  @assert dim(i) == dim(qns)
  return setdir(setspace(i, qns), dir)
end

function addqns(i::QNIndex, qns::QNBlocks)
  @assert dim(i) == dim(qns)
  @assert nblocks(qns) == nblocks(i)
  iqns = space(i)
  j = copy(i)
  jqn = space(j)
  for n in 1:nblocks(i)
    @assert blockdim(iqns, n) == blockdim(qns, n)
    iqn_n = qn(iqns, n)
    qn_n = qn(qns, n)
    newqn = iqn_n
    for nqv in 1:nactive(qn_n)
      qv = qn_n[nqv]
      newqn = addqnval(newqn, qv)
    end
    jqn[n] = newqn => blockdim(iqns, n)
  end
  return j
end

# Check that the QNs are all the same
function hassameflux(i1::QNIndex, i2::QNIndex)
  dim_i1 = dim(i1)
  dim_i1 ≠ dim(i2) && return false
  for n in 1:dim_i1
    flux(i1 => n) ≠ flux(i2 => n) && return false
  end
  return true
end

hassameflux(::QNIndex, ::Index) = false
hassameflux(::Index, ::QNIndex) = false

# Split the blocks into blocks of size 1 with the same QNs
splitblocks(i::Index) = setspace(i, splitblocks(space(i)))

function mutable_storage(::Type{Order{N}}, ::Type{IndexT}) where {N,IndexT<:QNIndex}
  return SizedVector{N,IndexT}(undef)
end

isfermionic(i::Index) = false

isfermionic(i::QNIndex) = any(q -> isfermionic(qn(q)), space(i))

function show(io::IO, i::QNIndex)
  idstr = "$(id(i) % 1000)"
  if length(tags(i)) > 0
    print(
      io, "(dim=$(dim(i))|id=$(idstr)|\"$(tagstring(tags(i)))\")$(primestring(plev(i)))"
    )
  else
    print(io, "(dim=$(dim(i))|id=$(idstr))$(primestring(plev(i)))")
  end
  println(io, " <$(dir(i))>")
  for (n, qnblock) in enumerate(space(i))
    print(io, " $n: $qnblock")
    n < length(space(i)) && println(io)
  end
end

function HDF5.write(parent::Union{HDF5.File,HDF5.Group}, name::AbstractString, B::QNBlocks)
  g = create_group(parent, name)
  attributes(g)["type"] = "QNBlocks"
  attributes(g)["version"] = 1
  write(g, "length", length(B))
  dims = [block[2] for block in B]
  write(g, "dims", dims)
  for n in 1:length(B)
    write(g, "QN[$n]", B[n][1])
  end
end

function HDF5.read(
  parent::Union{HDF5.File,HDF5.Group}, name::AbstractString, ::Type{QNBlocks}
)
  g = open_group(parent, name)
  if read(attributes(g)["type"]) != "QNBlocks"
    error("HDF5 group or file does not contain QNBlocks data")
  end
  N = read(g, "length")
  dims = read(g, "dims")
  B = QNBlocks(undef, N)
  for n in 1:length(B)
    B[n] = QNBlock(read(g, "QN[$n]", QN), dims[n])
  end
  return B
end


text:
using SparseTransforms
using Test

function test_codes()
    println("Testing LDPC code generation correctness")
    @test do_test(20, 25, 100)
end


text:
# Generate initial monolothic solution
function uhqF₀(u₀,h₀,q₀,F₀,X,Y,dΩ; mass_matrix_solver=Gridap.Algebra.BackslashSolver())
  a((u,p,r,u2),(v,q,s,v2))=∫(v⋅u+q*p+s*r+v2⋅u2)dΩ
  b((v,q,s,v2))=∫( v⋅u₀+ q*h₀ + s*q₀ + v2⋅F₀ )dΩ
  fes=FESolver(mass_matrix_solver)
  solve(fes,AffineFEOperator(a,b,X,Y))
end

"""
  Solves the nonlinear rotating shallow water equations
  T : [0,T] simulation interval
  N : number of time subintervals
  θ : Theta-method parameter [0,1)
  τ : APVM method stabilization parameter (dt/2 is typically a reasonable value)
"""
function shallow_water_theta_method_full_newton_time_stepper(
      nls::Gridap.Algebra.NonlinearSolver,
      model, order, degree, h₀, u₀, f₀, topography, g, θ, T, N, τ;
      mass_matrix_solver::Gridap.Algebra.LinearSolver=Gridap.Algebra.BackslashSolver(),
      matrix_type::Type{<:AbstractSparseMatrix}=SparseMatrixCSC{Float64,Int},
      vector_type::Type{<:AbstractVector}=Vector{Float64},
      write_diagnostics=true,
      write_diagnostics_freq=1,
      dump_diagnostics_on_screen=true,
      am_i_root=true,
      write_solution=false,
      write_solution_freq=N/10,
      output_dir="nswe_ncells_$(num_cells(model))_order_$(order)_theta_method_full_newton")

  Ω  = Triangulation(model)
  n  = get_normal_vector(Ω)
  dΩ = Measure(Ω,degree)
  dω = Measure(Ω,degree,ReferenceDomain())

  # Setup the trial and test spaces
  R, S, U, V, P, Q = setup_mixed_spaces(model, order)

  if (write_diagnostics)
    # assemble the mass matrices
    H1MM, _, L2MM, H1MMchol, RTMMchol, L2MMchol =
      setup_and_factorize_mass_matrices(dΩ, R, S, U, V, P, Q;
                                        mass_matrix_solver=mass_matrix_solver)
  end

  Y = MultiFieldFESpace([V,Q,S,V])
  X = MultiFieldFESpace([U,P,R,U])

  fes=FESolver(mass_matrix_solver)

  a1(u,v)=∫(v⋅u)dΩ
  l1(v)=∫(v⋅u₀)dΩ
  assem=SparseMatrixAssembler(matrix_type,vector_type,U,V)
  un=solve(fes,AffineFEOperator(a1,l1,U,V,assem))

  a2(u,v)=∫(v*u)dΩ
  l2(v)=∫(v*h₀)dΩ
  assem=SparseMatrixAssembler(matrix_type,vector_type,P,Q)
  hn=solve(fes,AffineFEOperator(a2,l2,P,Q,assem))

  a3(u,v)=∫(v*u)dΩ
  l3(v)=∫(v*f₀)dΩ
  assem=SparseMatrixAssembler(matrix_type,vector_type,R,S)
  fn=solve(fes,AffineFEOperator(a3,l3,R,S,assem))

  unv,hnv,fnv=get_free_dof_values(un,hn,fn)

  b = interpolate(topography,P)

  # Compute:
  #     - Initial potential vorticity (q₀)
  #     - Initial volume flux (F₀)
  #     - Initial full solution
  q₀=clone_fe_function(R,fn)
  compute_potential_vorticity!(q₀,H1MM,H1MMchol,dΩ,R,S,hn,un,fn,n,assem)
  F₀=clone_fe_function(V,un)
  compute_mass_flux!(F₀,dΩ,V,RTMMchol,un*hn)
  ΔuΔhqF=uhqF₀(un,hn,q₀,F₀,X,Y,dΩ; mass_matrix_solver=mass_matrix_solver)
  Δu,Δh,q,F = ΔuΔhqF

  h_tmp = copy(hnv)
  w_tmp = copy(fnv)

  function run_simulation(pvd=nothing)
    dt  = T/N
    τ   = dt/2 # APVM stabilization parameter
    hc  = CellField(h₀,Ω)
    uc  = CellField(u₀,Ω)

    if (write_diagnostics)
      ϕ  = clone_fe_function(Q,hn)
      wn = clone_fe_function(S,fn)
      diagnostics_file = joinpath(output_dir,"nswe_diagnostics_theta_method_full_newton.csv")
      if (am_i_root)
        initialize_csv(diagnostics_file,
                       "time", "mass", "vorticity", "kinetic", "potential", "power")
      end
    end

    for step=1:N
      function residual((Δu,Δh,qvort,F),(v,q,s,v2))
         one_m_θ = (1-θ)
         uiΔu    = un     + one_m_θ*Δu
         hiΔh    = hn     + one_m_θ*Δh
         hbiΔh   = hn + b + one_m_θ*Δh
         ∫((1.0/dt)*v⋅(Δu)-(∇⋅(v))*(g*hbiΔh + 0.5*uiΔu⋅uiΔu)+
             (qvort-τ*(uiΔu⋅∇(qvort)))*(v⋅⟂(F,n)) +   # eq1
           (1.0/dt)*q*(Δh))dΩ + ∫(q*(DIV(F)))dω +  # eq2
         ∫(s*qvort*hiΔh + ⟂(∇(s),n)⋅uiΔu - s*fn +   # eq3
             v2⋅(F-hiΔh*uiΔu))dΩ                      # eq4
       end

       function jacobian((Δu,Δh,qvort,F),(du,dh,dq,dF),(v,q,s,v2))
        one_m_θ = (1-θ)
        uiΔu  = un + one_m_θ*Δu
        hiΔh  = hn + one_m_θ*Δh
        uidu  = one_m_θ*du
        hidh  = one_m_θ*dh
        ∫((1.0/dt)*v⋅du +  (dq    - τ*(uiΔu⋅∇(dq)+uidu⋅∇(qvort)))*(v⋅⟂(F ,n))
                        +  (qvort - τ*(           uiΔu⋅∇(qvort)))*(v⋅⟂(dF,n))
                        -  (∇⋅(v))*(g*hidh +uiΔu⋅uidu)   +  # eq1
          (1.0/dt)*q*dh)dΩ + ∫(q*(DIV(dF)))dω             +  # eq2
          ∫(s*(qvort*hidh+dq*hiΔh) + ⟂(∇(s),n)⋅uidu       +  # eq3
            v2⋅(dF-hiΔh*uidu-hidh*uiΔu))dΩ                   # eq4
       end

       # Solve fully-coupled monolithic nonlinear problem
       # Use previous time-step solution, ΔuΔhqF, as initial guess
       # Overwrite solution into ΔuΔhqF

       assem = SparseMatrixAssembler(matrix_type,Vector{Float64},X,Y)
       op=FEOperator(residual,jacobian,X,Y,assem)
       solver=FESolver(nls)

       solve!(ΔuΔhqF,solver,op)

       # Update current solution
       unv .= unv .+ get_free_dof_values(Δu)
       hnv .= hnv .+ get_free_dof_values(Δh)

       if (write_diagnostics && write_diagnostics_freq>0 && mod(step, write_diagnostics_freq) == 0)
        compute_diagnostic_vorticity!(wn, dΩ, S, H1MMchol, un, n)
        compute_bernoulli_potential!(ϕ,dΩ,Q,L2MMchol,un⋅un,hn,g)
        dump_diagnostics_shallow_water!(h_tmp, w_tmp,
                                        model, dΩ, dω, S, L2MM, H1MM,
                                        hn, un, wn, ϕ, F, g, step, dt,
                                        diagnostics_file,
                                        dump_diagnostics_on_screen,
                                        am_i_root)
      end
      if (write_solution && write_solution_freq>0 && mod(step, write_solution_freq) == 0)
        if (!write_diagnostics || write_diagnostics_freq != write_solution_freq)
          compute_diagnostic_vorticity!(wn, dΩ, S, H1MMchol, un, n)
        end
        pvd[dt*Float64(step)] = new_vtk_step(Ω,joinpath(output_dir,"n=$(step)"),hn,un,wn)
      end
    end
    hn, un, num_free_dofs(Y)
  end
  if (am_i_root && (write_diagnostics || write_solution))
    rm(output_dir,force=true,recursive=true)
    mkdir(output_dir)
  end
  if (write_solution)
     pvdfile=joinpath(output_dir,
        "nswe_ncells_$(num_cells(model))_order_$(order)_theta_method_full_newton")
     paraview_collection(run_simulation,pvdfile)
   else
    run_simulation()
  end
end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


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

using Markdown
using InteractiveUtils

# ╔═╡ ca3d7c7e-e2ba-11ea-10a7-035583e2d4ad

# ╔═╡ 920f088e-e2ba-11ea-3e4b-a1ba9219af80
begin
	n = 10;
	m = 3
	A = rand(1:5,n,m)
	b = rand(100:20:200,n)
end

# ╔═╡ cda7f448-e2ba-11ea-0ad8-99a2231fe6cc
begin
	doc"""

		betahat(X, y)

	Estimates coefficients of a linear regression model using the `X\y`.

	# Arguments
	- `X`: `n` by `k` numeric array.
	- `y`: `n` by 1 numeric array.

	# Examples
	```julia-repl
	julia> n = 10
	julia> m = 3
	julia> A = rand(1:5,n,m)
	julia> b = rand(100:20:200,n)
	julia> betahat(X = A, y = b)
	```

	"""

	function betahat(X, y)
		return X\y
	end
end

# ╔═╡ 18ef23ae-e2bb-11ea-3302-79415e4dfc19
betahat(A, b)

# ╔═╡ 7541a2da-e2bd-11ea-2740-65d1f4f062f0
b

# ╔═╡ Cell order:
# ╠═ca3d7c7e-e2ba-11ea-10a7-035583e2d4ad
# ╠═920f088e-e2ba-11ea-3e4b-a1ba9219af80
# ╠═cda7f448-e2ba-11ea-0ad8-99a2231fe6cc
# ╠═18ef23ae-e2bb-11ea-3302-79415e4dfc19
# ╠═7541a2da-e2bd-11ea-2740-65d1f4f062f0

max_stars_repo_name: jeksterslabds/jeksterslablinreg

max_issues_repo_name: jeksterslabds/jeksterslablinreg

max_forks_repo_name: jeksterslabds/jeksterslablinreg


text:
@testset "$V{$T}" for V in [RefStdVector, GCStdVector, STL.SharedStdVector], T in sort!(collect(STL.StdVector_types); by=string)
    # These segfault; converting STL.SharedStdString to Ptr{StdString} doesn't work
    T ≡ RefStdVector{STL.SharedStdString} && continue
    T ≡ STL.SharedStdVector{STL.SharedStdString} && continue
    # @show V T

    # TODO: Need to free the generated RefStdString objects
    function mkT(i)
        return T ≡ Bool ? i % Bool :
               T <: Ptr ? T(i) :
               T ≡ RefStdString ? RefStdString(string(i)) :
               T ≡ STL.SharedStdString ? STL.SharedStdString(string(i)) :
               T ≡ RefStdVector{RefStdString} ? (x = T(1); x[0] = RefStdString(string(i)); x) :
               T ≡ RefStdVector{STL.SharedStdString} ? (x = T(1); x[0] = STL.SharedStdString(string(i)); x) :
               T ≡ STL.SharedStdVector{RefStdString} ? (x = T(1); x[0] = RefStdString(string(i)); x) :
               T ≡ STL.SharedStdVector{STL.SharedStdString} ? (x = T(1); x[0] = STL.SharedStdString(string(i)); x) : i
    end

    vec = V{T}()
    @test eltype(vec) ≡ T

    @test length(vec) == 0
    @test isempty(vec)

    if vec isa RefStdVector
        free(vec)
    end

    vec = V{T}(10)
    @test length(vec) == 10

    for i in 0:9
        vec[i] = mkT(i)
    end
    for i in 0:9
        @test vec[i] == mkT(i)
    end

    elts = convert(Vector{T}, vec)
    @test elts == [mkT(i) for i in 0:9]

    elts = [elt for elt in vec]
    @test typeof(elts) == Vector{T}
    @test elts == [mkT(i) for i in 0:9]

    resize!(vec, 20)
    @test length(vec) == 20

    for i in 0:9
        @test vec[i] == mkT(i)
    end

    empty!(vec)
    @test isempty(vec)

    if vec isa RefStdVector
        free(vec)
    end
    GC.gc(true)
end


text:
function solvejobs!(sp, jobsd, to::TimerOutput, stopatinf::Bool)
    infeasibility_detected = false
    for (node, jobs) in jobsd
        for job in jobs
            if !stopatinf || SOI.allfeasible(job.parent.pool)
                solvejob!(sp, job, node, to)
                if !SOI.allfeasible(job.parent.pool)
                    infeasibility_detected = true
                end
            end
        end
    end
    infeasibility_detected
end

function gencuts(pathsd, sp, to::TimerOutput, ztol)
    for (parent, paths) in pathsd
        for path in paths
            SOI.addcut!(sp, parent, path.pool, to, ztol)
        end
    end
end

function applycuts(pathsd, sp)
    for (node, _) in pathsd
        SOI.applycuts!(sp, node)
    end
end

"""
    struct Algorithm{ST<:AbstractPathSampler}
        K::Int
        verbose::Int
        pathsampler::ST
        ztol::Float64
        stopatinf::Bool
        mergepaths::Bool
        forwardcuts::Bool
        backwardcuts::Bool
    end

SDDP algorithm exploring `K` paths per iteration stages.
The paths will be selected according to `pathsampler` and equivalent paths might be merged if their difference is smaller than `ztol` and `mergepaths` is true.
The parameter `ztol` is also used to check whether a new cut is useful.
When a scenario is infeasible and `stopatinf` is true then no other scenario with the same ancestor is run. Note that since the order in which the different scenarios is run is not deterministic, this might introduce nondeterminism even if the sampling is deterministic.
By default, the cuts are added backward. However, if `forwardcuts` is set to `true` and `backwardcuts` is set to `false` the cuts are added forward.
"""
struct Algorithm{ST<:AbstractPathSampler} <: SOI.AbstractAlgorithm
    K::Int
    pathsampler::ST
    ztol::Float64
    stopatinf::Bool
    mergepaths::Bool
    forwardcuts::Bool
    backwardcuts::Bool
end
function Algorithm(; K::Int=25, pathsampler::AbstractPathSampler=ProbaPathSampler(true), ztol=1e-6, stopatinf=false, mergepaths=true, forwardcuts=false, backwardcuts=true)
    Algorithm(K, pathsampler, ztol, stopatinf, mergepaths, forwardcuts, backwardcuts)
end

struct Paths{NodeT, TT, SolT} <: SOI.AbstractPaths
    paths::Vector{Vector{Tuple{NodeT, Vector{SDDPPath{TT, SolT}}}}}
end
function SOI.npaths(paths::Paths)
    # At the master node, all paths are still concentrated in the same one
    node_paths = paths.paths[1]
    @assert length(node_paths) == 1
    paths_vec = node_paths[1][2]
    @assert length(paths_vec) == 1
    Ks = paths_vec[1].K
    @assert length(Ks) == 1
    return Ks[1]
end

function SOI.forward_pass!(sp::SOI.AbstractStochasticProgram, algo::Algorithm, to::TimerOutput, result::SOI.Result, verbose)
    master = SOI.get(sp, SOI.MasterNode())
    NodeT = typeof(master)
    @timeit to "solve" mastersol = SOI.get(sp, SOI.Solution(), master)
    infeasibility_detected = SOI.getstatus(mastersol) == :Infeasible

    num_stages = SOI.get(sp, SOI.NumberOfStages())
    PathT = SDDPPath{SOI.get(sp, SOI.TransitionType()), typeof(mastersol)}
    TT = Tuple{NodeT, Vector{PathT}}
    pathsd = Vector{Vector{TT}}(undef, num_stages)
    if infeasibility_detected
        pathsd[1] = TT[]
    else
        pathsd[1] = TT[(master, [SDDPPath{SOI.get(sp, SOI.TransitionType())}(mastersol, [SOI.getnodeobjectivevalue(mastersol)], [1.], [algo.K])])]
    end
    endedpaths = PathT[]

    for t in 2:num_stages
        verbose >= 3 && println("Stage $t/$num_stages")

        # Merge paths
        if algo.mergepaths
            pathsd[t-1] = mergesamepaths(pathsd[t-1], to, algo.ztol)
        end

        # Make jobs
        jobsd = childjobs(sp, pathsd[t-1], algo.pathsampler, t, num_stages, endedpaths)

        # Solve Jobs (parallelism possible here)
        infeasibility_detected |= solvejobs!(sp, jobsd, to, algo.stopatinf)

        if algo.forwardcuts
            gencuts(pathsd[t-1], sp, to, algo.ztol)
            # The cut is added after so that they are added by group and duplicate detection in CutPruners works better
            applycuts(pathsd[t-1], sp)
        end

        # Jobs -> Paths
        pathsd[t] = jobstopaths(jobsd, sp)
    end

    if infeasibility_detected
        z_UB = Inf # FIXME assumes minimization
        σ = 0
    else
        for (_, paths) in pathsd[end]
            append!(endedpaths, paths)
        end
        z_UB, σ = meanstdpaths(endedpaths, algo.K)
    end

    # update stats
    result.upperbound = z_UB
    result.σ_UB = σ
    result.paths = Paths(pathsd)
    result.status = SOI.getstatus(mastersol)             # TODO status in backward pass
    result.lowerbound = SOI.getobjectivevalue(mastersol) # TODO lowerbound in backward pass
end

function SOI.backward_pass!(sp::SOI.AbstractStochasticProgram, algo::Algorithm, to::TimerOutput, result::SOI.Result, verbose)
    if algo.backwardcuts
        # We could ask a node if it has new fcuts/ocuts before resolving, that way the last stage won't be a particular case anymore
        num_stages = SOI.get(sp, SOI.NumberOfStages())
        for t in num_stages:-1:2
            # The last stages do not need to be resolved since it does not have any new cut
            if t != num_stages
                # Make jobs
                endedpaths = SDDPPath{SOI.get(sp, SOI.TransitionType())}[]
                jobsd = childjobs(sp, result.paths.paths[t-1], algo.pathsampler, t, num_stages, endedpaths) # FIXME shouldn't need pathsampler here
                # Solve Jobs (parallelism possible here)
                solvejobs!(sp, jobsd, to, algo.stopatinf)
            end

            gencuts(result.paths.paths[t-1], sp, to, algo.ztol)
            # The cut is added after so that they are added by group and duplicate detection in CutPruners works better
            applycuts(result.paths.paths[t-1], sp)
        end
    end
end


text:
"""
only works when t's all the same for outputs
"""
function plotgp(t::Array{Array{<:AbstractFloat,1},1} , gp::GaussianProcess;
                samps=false, test=nothing, jitter=1e-6, save=nothing)

    μ, K = posterior(t, gp)

    dist = MvNormal(μ, K + jitter * I)

    μ_arr = split_outputs(μ, t)

    K_arr = split_outputs(sqrt.(diag(K + jitter * I)), t)

    p = plot(t, μ_arr, legend=true, layout=gp.D)
    plot!(t, μ_arr + 2 * K_arr, legend=false, lc="red", ls=:dot)
    plot!(t, μ_arr - 2 * K_arr, legend=false, lc="red", ls=:dot)
    scatter!(gp.data.X, hcat(gp.data.Y...),
        markershape=:o,
        markercolor="blue",
        markersize=1,
        markerstrokewidth=0)

    if test !== nothing
        scatter!(test.X, hcat(test.Y...),
        markershape=:o,
        markercolor="green",
        markersize=1,
        markerstrokewidth=0)
    end

    if samps
        N = 20
        for i = 1:N
            s = rand(dist)
            s = split_outputs(s, t)
            plot!(t, s, legend=false, lc="gray", lw=0.2)
        end

    end

    if save !== nothing
        savefig("/Users/magnus/Documents/phd/code/repos/VolterraGP/bin/plots/$save.svg")
    end
    display(p)
end


text:
using Distributions
using ForwardDiff
using Optim
using Pathfinder
using Random
using Test

include("test_utils.jl")

@testset "MvNormal functions" begin
    @testset "fit_mvnormals" begin
        n = 10
        logp(x) = logp_banana(x)
        ∇logp(x) = ForwardDiff.gradient(logp, x)
        θ₀ = 10 * randn(n)
        optimizer = Optim.LBFGS()
        θs, logpθs, ∇logpθs = Pathfinder.maximize_with_trace(logp, ∇logp, θ₀, optimizer)
        Σs = Pathfinder.lbfgs_inverse_hessians(θs, ∇logpθs; history_length=optimizer.m)
        dists = @inferred Pathfinder.fit_mvnormals(θs, ∇logpθs; history_length=optimizer.m)
        @test dists isa Vector{<:MvNormal{Float64,<:Pathfinder.WoodburyPDMat}}
        @test Σs ≈ getproperty.(dists, :Σ)
        @test θs .+ Σs .* ∇logpθs ≈ getproperty.(dists, :μ)
    end

    @testset "rand_and_logpdf" begin
        ndraws = 20
        @testset "MvNormal" begin
            n = 10
            Σ = rand_pd_mat(Float64, 10)
            μ = randn(n)
            dist = MvNormal(μ, Σ)

            rng = MersenneTwister(42)
            x, logpx = @inferred Pathfinder.rand_and_logpdf(rng, dist, ndraws)
            rng = MersenneTwister(42)
            x2 = rand(rng, dist, ndraws)
            logpx2 = logpdf(dist, x2)
            @test x ≈ x2
            @test logpx ≈ logpx2
        end

        @testset "MvNormal{T,Pathfinder.WoodburyPDMat{T}}" begin
            n = 10
            ndraws = 20
            nhist = 4
            A = rand_pd_diag_mat(Float64, 10)
            D = rand_pd_mat(Float64, 2nhist)
            B = randn(n, 2nhist)
            Σ = Pathfinder.WoodburyPDMat(A, B, D)
            μ = randn(n)
            dist = MvNormal(μ, Σ)

            rng = MersenneTwister(42)
            x, logpx = @inferred Pathfinder.rand_and_logpdf(rng, dist, ndraws)
            rng = MersenneTwister(42)
            x2 = rand(rng, dist, ndraws)
            logpx2 = logpdf(dist, x2)
            @test x ≈ x2
            @test logpx ≈ logpx2
        end

        @testset "Normal" begin
            σ = rand() * 10
            μ = randn()
            dist = Normal(μ, σ)

            rng = MersenneTwister(42)
            x, logpx = @inferred Pathfinder.rand_and_logpdf(rng, dist, ndraws)
            rng = MersenneTwister(42)
            x2 = rand(rng, dist, ndraws)
            logpx2 = logpdf.(dist, x2)
            @test x ≈ x2'
            @test logpx ≈ logpx2
        end
    end
end


text:
module FlowChart

end # module


text:
module System
	using CBinding

	export @sys
	export @c_str, isqualifiedwith, unqualifiedtype, bitstype

	macro sys(exprs...) return sys(__module__, __source__, exprs...) end

	function sys(mod::Module, loc::LineNumberNode, expr::Expr)
		Base.is_expr(expr, :using, 1) || Base.is_expr(expr, :import, 1) || error("Expected a `using ...` or `import ...` statement in @sys usage")

		path = expr.args[1]
		if Base.is_expr(path, :(:)) && length(path.args) >= 1
			path = path.args[1]
		end
		name = path.args[1]

		sym = gensym(name)
		pushfirst!(path.args, :Base, :Main, sym)

		return quote
			@eval Base.Main module $(sym)
				import $(name)
			end
			$(expr)
		end
	end

	function __init__()
		push!(Base.LOAD_PATH, joinpath(dirname(@__DIR__), "pkgs"))
	end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 5955f292-0ee4-4574-a184-0a1e12b68101
begin
	using SymPy,LinearAlgebra
end

# ╔═╡ b78cb61c-c5bc-11ec-26a4-1bd5b2742132
md"# SymPy"

# ╔═╡ 79b66a3a-1dad-4308-a5cb-ff0fc82d40dd
begin
	n=2
	m=3
	A=Array{Sym,2}(undef,m,n)
	for i in 1:m
		for j in 1:n
		A[i,j]=Sym("a$i$j")
		end
	end
	A
end

# ╔═╡ faf1d18c-900c-4fa6-a4bd-5982c8486975
begin
	B=[Sym("a$i$j") for i=1:3 for j=1:2].reshape(3,2)
	B
end

# ╔═╡ 60c1690d-2ab7-4090-9175-0ac1a5198702
begin
	B.T
end

# ╔═╡ 044ae63d-dc68-4690-b3ae-b9c712bd71f6
begin
	Sym[1 -1; 3 4; 0 2]
end

# ╔═╡ b5462599-e657-4e39-9119-dca930f50b97
Sym[1, 2, 3]

# ╔═╡ 91529d3f-4e45-4c05-9ac3-09947549f936
begin
	M = Sym[1 2 3; 3 2 1]
	N = Sym[0, 1, 1]
    M*N

end

# ╔═╡ 140fa211-495f-4150-886d-ab693b206cfc
M.row(0)

# ╔═╡ b98ef5d4-5d47-4ee5-9e22-b235d41bb89d
M.col(0)

# ╔═╡ 14aee5b2-8496-4adf-b450-f1e69d3386de
sympy.eye(4)

# ╔═╡ 241347f0-d3d1-41d9-b09c-84ea4ca4fa03
sympy.zeros(2, 3)

# ╔═╡ db255a23-7f51-4b07-9e9b-288f968c716e
sympy.diag(1, 2, 3)

# ╔═╡ bab2ccc7-82b6-42fe-bc81-b31fb9fb98a5
let
	M = Sym[1  0  1 3; 2 3 4 7; -1 -3 -3 -4]
	M.rref()
end

# ╔═╡ 8e7b1930-cbd3-402e-ad42-78f0a126bacf
let
	M = Sym[1 1 2; 2 1 3; 3 1 4]
	M.columnspace()
end

# ╔═╡ f135e32a-320a-4e64-b02b-243910562607
let
	M = Sym[3 -2  4 -2; 5  3 -3 -2; 5 -2  2 -2; 5 -2 -3  3]
	#M.eigenvals()
	M.eigenvects()
end

# ╔═╡ 64c8015a-aab2-4edb-ad45-82bca24c4c2f
function makeSymArray(m::Number, n::Number, elname::String="a")
	"""
	 生成m行, n列的矩阵
	"""
    A = [Sym("$elname$i$j") for i=1:m for j=1:n].reshape(m,n)
    return A
end

# ╔═╡ c4adfad5-3558-4260-a1c2-261cd93d6662
begin
	C=makeSymArray(3,2,"c")
	C

end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"

[compat]
SymPy = "~1.1.4"
"""

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

julia_version = "1.7.1"
manifest_format = "2.0"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

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

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.14.0"

[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

[[deps.CommonEq]]
git-tree-sha1 = "d1beba82ceee6dc0fce8cb6b80bf600bbde66381"
uuid = "3709ef60-1bee-4518-9f2f-acd86f176c50"
version = "0.2.0"

[[deps.CommonSolve]]
git-tree-sha1 = "68a0743f578349ada8bc911a5cbd5a2ef6ed6d1f"
uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
version = "0.2.0"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "96b0bc6c52df76506efc8a441c6cf1adcb1babc4"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.42.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.Conda]]
deps = ["Downloads", "JSON", "VersionParsing"]
git-tree-sha1 = "6e47d11ea2776bc5627421d59cdcc1296c058071"
uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d"
version = "1.7.0"

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

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

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

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "91b5dcf362c5add98049e6c29ee756910b03051d"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.3"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[deps.LaTeXStrings]]
git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.3.0"

[[deps.Latexify]]
deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
git-tree-sha1 = "4f00cc36fede3c04b8acf9b2e2763decfdcecfa6"
uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
version = "0.15.13"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "58f25e56b706f95125dcb796f39e1fb01d913a71"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.10"

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

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

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

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

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

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.3"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "d3538e7f8a790dc8903519090857ef8e1283eecd"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.5"

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

[[deps.PyCall]]
deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"]
git-tree-sha1 = "1fc929f47d7c151c839c5fc1375929766fb8edcc"
uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
version = "1.93.1"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

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

[[deps.RecipesBase]]
git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.2.1"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

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

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

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.1.4"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.SymPy]]
deps = ["CommonEq", "CommonSolve", "Latexify", "LinearAlgebra", "Markdown", "PyCall", "RecipesBase", "SpecialFunctions"]
git-tree-sha1 = "1763d267a68a4e58330925b7ce8b9ea2ec06c882"
uuid = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
version = "1.1.4"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

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

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

[[deps.VersionParsing]]
git-tree-sha1 = "58d6e80b4ee071f5efd07fda82cb9fbe17200868"
uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
version = "1.3.0"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═b78cb61c-c5bc-11ec-26a4-1bd5b2742132
# ╠═5955f292-0ee4-4574-a184-0a1e12b68101
# ╠═79b66a3a-1dad-4308-a5cb-ff0fc82d40dd
# ╠═faf1d18c-900c-4fa6-a4bd-5982c8486975
# ╠═60c1690d-2ab7-4090-9175-0ac1a5198702
# ╠═c4adfad5-3558-4260-a1c2-261cd93d6662
# ╠═044ae63d-dc68-4690-b3ae-b9c712bd71f6
# ╠═b5462599-e657-4e39-9119-dca930f50b97
# ╠═91529d3f-4e45-4c05-9ac3-09947549f936
# ╠═140fa211-495f-4150-886d-ab693b206cfc
# ╠═b98ef5d4-5d47-4ee5-9e22-b235d41bb89d
# ╠═14aee5b2-8496-4adf-b450-f1e69d3386de
# ╠═241347f0-d3d1-41d9-b09c-84ea4ca4fa03
# ╠═db255a23-7f51-4b07-9e9b-288f968c716e
# ╠═bab2ccc7-82b6-42fe-bc81-b31fb9fb98a5
# ╠═8e7b1930-cbd3-402e-ad42-78f0a126bacf
# ╠═f135e32a-320a-4e64-b02b-243910562607
# ╠═64c8015a-aab2-4edb-ad45-82bca24c4c2f
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: phpsmarter/plutonote

max_issues_repo_name: phpsmarter/plutonote

max_forks_repo_name: phpsmarter/plutonote


text:
## A class to hold them all!
# struct rod_problem
#     ns::Vector{Int}
#     bc::Vector{rod_bc}
# end

"""
    discretize_span(sspan, ns)

Return the vertex and edge locations (`sv` and `se`) for `sspan` divided into `ns` segments.

# Examples
```julia-repl
julia> sv, se = discretize_span([0,1], 10)
```
"""
function discretize_span(sspan, ns)
    ds = (sspan[2] - sspan[1])/ns
    sv = [sspan[1]; LinRange(sspan[1]+ds/2, sspan[2]-ds/2, ns); sspan[2]]
    se = LinRange(sspan[1], sspan[2], ns+1)
    sv, se
end

"""
    straight_rod(p1, p2, d, N)

Returns a rod with endpoints `p1`, `p2` and director `d` with `N` segments.

If `d` has a component parallel to `p2 - p1`, it is projected out.

"""
function straight_rod(p1::Array{T,1}, p2::Array{T,1}, d::Array{T,1}, N::Int) where T

    # X coordinates
    L = norm(p2-p1)
    t = (p2-p1)/L
    sv, ~ = discretize_span([0,L], N)
    xmat = p1 .+ t*sv'

    # Directors
    d_projected = d - dot(d, t)*t
    d_projected /= norm(d_projected)
    dmat = repeat(d_projected, 1, N+1)

    basic_rod(xmat, dmat)
end

# Helper for continuous_rod
using OrdinaryDiffEq
function bishop_system(vec, p, s)
    # κ1_fun, κ2_fun, τ_fun = p
    #κ1, κ2, τ = κ1_fun(s), κ2_fun(s), τ_fun(s)
    κ1, κ2, τ = p[1](s), p[2](s), p[3](s)

    x = @view vec[1:3]
    t = @view vec[4:6]
    d1 = @view vec[7:9]
    d2 = cross(t, d1)

    [t; κ1*d1 + κ2*d2; -κ1*t + τ*d2]
end

"""
    continuous_rod(x0, t0, d0, κ1, κ2, τ, s)

Returns two functions x and d defining the centerline and directors of a rod.

These are determined by solving the ODE which defines a Bishop frame with
curvatures `κ1(s)`, `κ2(s)`, and `τ(s)`. The integration starts from initial conditions
`x0`, `t0`, and `d0`, and covers length `s = [s1, s2]`

# Examples
```julia-repl
julia> xfun, dfun = continuous_rod([0,0,0], [1,0,0], [0,0,1],
                                      s -> 2, s -> 1, s -> s*(1-s), [0,1])
```
"""
function continuous_rod(x0, t0, d0, κ1::Function, κ2::Function, τ::Function, s)
    p = (κ1, κ2, τ)
    prob = ODEProblem(bishop_system, [x0;t0;d0], s, p)
    sol = solve(prob, Tsit5(), abs_tol = 1e-8, rel_tol = 1e-8)

    xfun(s) = sol(s)[1:3]
    dfun(s) = sol(s)[7:9]
    xfun, dfun
end

"""
    discrete_rod(xfun, dfun, s, N)

Return a rod with centerline and directors approximating xfun(s) and dfun(s)

Directors and centerline are sampled from the continous curve, and then
any tangential component of the directors is removed by projection.

# Examples
```julia-repl
julia> r = discrete_rod([0,0,0], [1,0,0], [0,0,1], s -> 2, s -> 1, s -> s*(1-s), [0,1])
```
"""
function discrete_rod(xfun::Function, dfun::Function, s, N::Int)

    # Sample Directors
    sv, se = discretize_span(s, N)
    x = Matrix{Float64}(undef, 3, N+2)
    d = Matrix{Float64}(undef, 3, N+1)
    for i in 1:N+2; x[:,i] = xfun(sv[i]); end
    for i in 1:N+1; d[:,i] = dfun(se[i]); end

    # Project and normalize directors
    t = tangents(x)
    d = d .- dot3(d,t).*t
    d = d./norm3(d)

    basic_rod(x,d)
end

"""
    elliptical_stiffness(E, ν, a, b)

Return the elastic properties of a rod with elliptical section.

# Arguments
- `E`: Elastic Modulus
- `ν`: Poisson Ratio
- `a`, `b`: Radii of the cross section
"""
function elliptical_stiffness(E, ν, a, b)
     G = E/(2*(1 + ν))
     A = π*a*b
     k = E*A
     B = E*A/4*[a^2 b^2]
     β = G*A*(a^2 + b^2)/4
     scalar_rod_properties(k, B, β)
end

"""
    random_rod(ns)

Return an admissible rod with ns segments. Centerline and directors are random.

"""
function random_rod(ns::Int)

    x = randn(ns+2, 3)
    d = randn(ns+1, 3)

    t = tangents(x)
    d = d .- dot3(d,t).*t
    norm3!(d)

    basic_rod(x,d)
end

"""
    random_delta(ns)

Return an admissible rod_delta for a rod with ns segments.

"""
function random_delta(ns::Int)

  x = randn(ns+2, 3)
  θ = randn(ns+1, 1)
  rod_delta(x,θ)

end

"""
   arc_rod(p1, p2, d, r, ns)

Return a rod spanning an arc from p1 to p2 with radius r and directors d.

Discretization is done with ns segments. The director should be orthogonal
to p2 - p1, and ther resulting arc will lie in the plane orthogonal to d.

"""

function arc_rod(p1, p2, d, r, ns)
    l = norm(p2 - p1)
    e1 = (p2 - p1)/l
    e2 = sign(r)*cross(e1,d)
    R = abs(r)

    θ = asin(l/(2*R))
    C = (p1 + p2)/2 - R*cos(θ)*e2

    sv, se = discretize_span([-θ, θ], ns)
    vertices = C' .+ R*sin.(sv).*e1' + R*cos.(sv).*e2'
    directors = repeat(d', ns+1)
    basic_rod(vertices, directors)
end

# ------------------------------------------------------------------------------
# Specification Layer Stuff
# ------------------------------------------------------------------------------

using LinearAlgebra
import Base.length

abstract type rod_bc end

struct free_bc <: rod_bc
    rod_num::Int
    rod_end::Int
end

struct static_clamp_bc <: rod_bc
    rod_num::Int
    rod_end::Int
end

num_free_parameters(bc::free_bc) = 7
num_free_parameters(bc::static_clamp_bc) = 0

# How to apply the BCs
function apply_bc(T::Type, bc::free_bc, params)
    return hcat(params[1:3], params[5:7]), params[4]
end

function apply_bc(T::Type, bc::static_clamp_bc, params)
    return zeros(T,3,2), zero(T)
end

# Compile unknowns and specification layer
function get_specification_layer(rods::Vector{Int}, bcs::Vector{rod_bc})

    # Define size of input vector
    nq = sum((4 .*rods.-7)) + sum(num_free_parameters, bcs)

    # Define a function mapping q to rod_deltas
    function param_2_rod(q::Vector{T}) where T

        # Distribute core degrees of freedom
        Δx_cores = Vector{Matrix{T}}()
        Δθ_cores = Vector{Matrix{T}}()
        index = 0
        for ns in rods
            ncore = 4*ns-7
            s_ind = index
            e_ind = index + ncore
            push!(Δx_cores, vcat(reshape(q[s_ind+2:4:e_ind], 1, ns-2),
                                 reshape(q[s_ind+3:4:e_ind], 1, ns-2),
                                 reshape(q[s_ind+4:4:e_ind], 1, ns-2)))
            push!(Δθ_cores, reshape(q[s_ind+1:4:e_ind], 1, ns-1))
            index += ncore
        end

        # Work on boundary degrees of freedom
        Δx_start = Vector{Matrix{T}}(undef, length(rods))
        Δθ_start = Vector{Vector{T}}(undef, length(rods))
        Δx_end = Vector{Matrix{T}}(undef, length(rods))
        Δθ_end = Vector{Vector{T}}(undef, length(rods))

        for bc in bcs
            nparams = num_free_parameters(bc)
            s_ind = index + 1
            e_ind = index + nparams
            Δx , Δθ = apply_bc(T, bc, q[s_ind:e_ind])
            if (bc.rod_end == 0)
                Δx_start[bc.rod_num] = Δx
                Δθ_start[bc.rod_num] = [Δθ]
            else
                Δx_end[bc.rod_num] = Δx
                Δθ_end[bc.rod_num] = [Δθ]
            end
            index += nparams
        end

        # Concatenate everything
        #deltas = rod_delta(hcat(Δx_start[1], Δx_cores[1], Δx_end[1]), hcat(Δθ_start[1], Δθ_cores[1], Δθ_end[1]))
        deltas = rod_delta(zeros(T,3,4), zeros(T,3,4))
                  #rod_delta(hcat(Δx_start[2], Δx_cores[2], Δx_end[2]), hcat(Δθ_start[2], Δθ_cores[2], Δθ_end[2]))]
        deltas

        #[rod_delta(hcat(Δx_start[i], Δx_cores[i], Δx_end[i]), hcat(Δθ_start[i], Δθ_cores[i], Δθ_end[i])) for i in 1:2]
    end

    return nq, param_2_rod
end

using Zygote:@adjoint
@adjoint basic_rod(x, d) = basic_rod(x, d), dr -> (dr.x, dr.d)
@adjoint rod_delta(Δx, Δθ) = rod_delta(Δx, Δθ), dr -> (dr.Δx, dr.Δθ)
@adjoint rod_strains(l, κ, τ) = rod_strains(l, κ, τ), ds -> (ds.l, ds.κ, ds.τ)

# ------------------------------------------------------------------------------
# Rotations via MRPs
# ------------------------------------------------------------------------------
using Rotations
using StaticArrays
using LinearAlgebra
function spin_angle(R::Rotation{3,T}, t) where T
    Rt = R*t
    gP = RodriguesParam(SVector{3}(cross(t, R*t)/(1 + dot(t,R*t))))
    gR = RodriguesParam(R)
    2*atan(norm(gR*gP))
end


text:
# NOTE: This is one-time defines to avoid name name conflicts with the user defined type names.

# `trigger` field  type symbol. Common to all AbstractComponent
const TRIGGER_TYPE_SYMBOL = gensym()

# `handshake` field type symbol. Commona to all AbstractComponent.
const HANDSHAKE_TYPE_SYMBOL = gensym()

# `callbacks` field type symbol. Commone to all AbstractComponent
const CALLBACKS_TYPE_SYMBOL = gensym()

# `id` field type symbol. Common to all AbstractComponent.
const ID_TYPE_SYMBOL = gensym()

# `modelargs` field type symbol. Common to all AbstractDynamicalSystem
const MODEL_ARGS_TYPE_SYMBOL = gensym()

# `modelkwargs` field type symbol. Common to all AbstractDynamicalSystem
const MODEL_KWARGS_TYPE_SYMBOL = gensym()

# `solverargs` field type symbol. Common to all AbstractDynamicalSystem
const SOLVER_ARGS_TYPE_SYMBOL = gensym()

# `solverkwargs` field type symbol. Common to all AbstractDynamicalSystem
const SOLVER_KWARGS_TYPE_SYMBOL = gensym()

# `alg` field type symbol. Common to all AbstractDynamicalSystem
const ALG_TYPE_SYMBOL = gensym()

# `integrator` field type symbol. Common to all AbstractDynamicalSystem
const INTEGRATOR_TYPE_SYMBOL = gensym()

# `input` field type symbol. Common to all AbstractSink
const INPUT_TYPE_SYMBOL = gensym()

# `plugin` field type symbol. Common to all AbstractSink
const PLUGIN_TYPE_SYMBOL = gensym()

# `timebuf` field type symbol. Common to all AbstractSink
const TIMEBUF_TYPE_SYMBOL = gensym()

# `databuf` field type symbol. Common to all AbstractSink
const DATABUF_TYPE_SYMBOL = gensym()

# `sinkcallback` field type symbol. Common to all AbstractSink
const SINK_CALLBACK_TYPE_SYMBOL = gensym()


text:
# Profiler interface

export
    @cuprofile

@enum(CUoutput_mode, CSV            = 0x00,
                     KEY_VALUE_PAIR = 0x01)

start_profiler() = @apicall(:cuProfilerStart, (Ptr{Void},), C_NULL)
stop_profiler() = @apicall(:cuProfilerStop, (Ptr{Void},), C_NULL)

macro cuprofile(ex)
    quote
        start_profiler()
        try
            $(esc(ex))
        finally
            stop_profiler()
        end
    end
end


text:
module OpenTelemetryProto

include("generated/OpentelemetryClients.jl")

end


text:
using LinearAlgebraicRepresentation
L = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL
# using Revise

# square=([[0; 0] [0; 1] [1; 0] [1; 1]], [[1, 2, 3, 4]],
# [[1,2], [1,3], [2,4], [3,4]])
# V,FV,EV = square
# model = V,([[1],[2],[3],[4]],EV,FV)
# table=L.apply(L.t(-0.5,-0.5), square)
# chair=L.Struct([L.t(0.75,0),L.s(0.35,0.35),table])
# structo=L.Struct([L.t(2,1),table,repeat([L.r(pi/2),chair],outer=4)...])
# structo1=L.Struct(repeat([structo,L.t(0,2.5)],outer=10));
# structo2=L.Struct(repeat([structo1,L.t(3,0)],outer=10));
#
# V,FV,EV = L.struct2lar(structo2);
# GL.VIEW([ GL.GLHulls(V,FV, GL.COLORS[1]) ]);
# #GL.VIEW([ GL.GLHulls(V,EV, GL.COLORS[1]) ]);  # BUG with colors?
#
# assembly = (L.struct2lar ∘ L.Struct)([L.sphere()(), L.t(3,0,-1), L.cylinder()()])
# GL.VIEW([ GL.GLHulls(L.cylinder()()...,GL.Point4d(1,1,1,0.2)) ]);
# GL.VIEW([ GL.GLHulls(assembly...,GL.Point4d(1,1,1,0.2)) ]);

cube = L.apply( L.t(-.5,-.5,0), L.cuboid([1,1,1]))
tableTop = L.Struct([ L.t(0,0,.85), L.s(1,1,.05), cube ])
tableLeg = L.Struct([ L.t(-.475,-.475,0), L.s(.1,.1,.89), cube ])
tablelegs = L.Struct( repeat([ tableLeg, L.r(0,0,pi/2) ],outer=4) )
table = L.Struct([ tableTop, tablelegs ])
table = L.struct2lar(table)
GL.VIEW([ GL.GLPol(table..., GL.Point4d(0,1,0,0.1)) ]);

cylndr = L.rod(.06, .5, 2*pi)([8,1])
chairTop = L.Struct([ L.t(0,0,0.5), L.s(0.5,0.5,0.04), cube ])
chairLeg = L.Struct([ L.t(-.22,-.22,0), L.s(.5,.5,1), L.r(0,0,pi/8), cylndr ])
chairlegs = L.Struct( repeat([ chairLeg, L.r(0,0,pi/2) ],outer=4) );
chair = L.Struct([ chairTop, chairlegs ]);
chair = L.struct2lar(chair)
GL.VIEW([ GL.GLPol(chair..., GL.Point4d(1,1,1,0.2)) ]);

theChair = L.Struct([ L.t(-.8,0,0), chair ])
fourChairs = L.Struct( repeat([L.r(0,0,pi/2), theChair],outer=4) );
fourSit = L.Struct([fourChairs,table]);
fourSit = L.struct2lar(fourSit)
GL.VIEW([ GL.GLPol(fourSit..., GL.Point4d(1,1,1,0.2)) ]);

singleRow=L.Struct(repeat([fourSit,L.t(0,2.5,0)],outer=10));
singleRow = L.struct2lar(singleRow)
GL.VIEW([ GL.GLPol(singleRow..., GL.Point4d(1,1,1,0.2)) ]);

refectory=L.Struct(repeat([singleRow,L.t(3,0,0)],outer=10));
refectory = L.struct2lar(refectory)
GL.VIEW([ GL.GLPol(refectory..., GL.Point4d(1,1,1,0.2)) ]);


text:
struct TCPHeader
    src_port::UInt16
    dst_port::UInt16
    seq_num::UInt32
    ack_num::UInt32
    data_offset::UInt8
    flags::UInt8
    window_size::UInt16
    checksum::UInt16
    urgent_pointer::UInt16
end

const FIN_MASK = 0x01
const SYN_MASK = 0x01 << 1
const RST_MASK = 0x01 << 2
const PSH_MASK = 0x01 << 3
const ACK_MASK = 0x01 << 4
const URG_MASK = 0x01 << 5

function Base.getproperty(h::TCPHeader, p::Symbol)
    p == :FIN && return (h.flags & FIN_MASK) != 0
    p == :SYN && return (h.flags & SYN_MASK) != 0
    p == :RST && return (h.flags & RST_MASK) != 0
    p == :PSH && return (h.flags & PSH_MASK) != 0
    p == :ACK && return (h.flags & ACK_MASK) != 0
    p == :URG && return (h.flags & URG_MASK) != 0
    getfield(h, p)
end

struct TCPPacket
    header::TCPHeader
    payload::UnsafeArray{UInt8, 1}
end

function decode_tcp(data::DenseVector{UInt8})
    p = Base.unsafe_convert(Ptr{UInt8}, data)
    GC.@preserve data begin
        h = unsafe_load(convert(Ptr{TCPHeader}, p))
    end
    h = TCPHeader(
        ntoh(h.src_port),
        ntoh(h.dst_port),
        ntoh(h.seq_num),
        ntoh(h.ack_num),
        h.data_offset >> 4,
        h.flags,
        ntoh(h.window_size),
        ntoh(h.checksum),
        ntoh(h.urgent_pointer),
    )
    TCPPacket(h, UnsafeArray{UInt8, 1}(p + h.data_offset *4, (Int(length(data) - h.data_offset * 4),)))
end

max_stars_repo_name: ancapdev/NetworkProtocols

max_issues_repo_name: ancapdev/NetworkProtocols

max_forks_repo_name: ancapdev/NetworkProtocols


text:
using Wkhtmltox
using Base.Test

#### conversion to pdf ####

exfile = joinpath(dirname(@__FILE__),"../examples/index.html")
# exfile = joinpath(dirname(@__FILE__),"../examples/example.html")
outfile = tempname() * ".pdf"

pdf_init(0)

ps = PdfSettings("out" => outfile,
                 "orientation" => "Landscape")

@test ps["out"] == outfile
@test ps["orientation"] == "Landscape"

ps["size.pageSize"] = "A4"
@test ps["size.pageSize"] == "A4"

os = PdfObject("page" => exfile,
               "web.javascript" => "true",
               "load.proxy" => "http://proxy-mkt.int.world.socgen:8080")
@test os["page"] == exfile

os["web.defaultEncoding"]

os2 = PdfObject()
@test os2["page"] == ""
os2["page"] = exfile
@test os2["page"] == exfile

conv = PdfConverter(ps)
push!(conv, os)

run(conv) == 1

isfile(outfile)

conv = nothing
pdf_deinit()

# wkhtmltopdf_set_progress_changed_callback(c, progress_changed);
# wkhtmltopdf_set_phase_changed_callback(c, phase_changed);
# wkhtmltopdf_set_error_callback(c, error);
# wkhtmltopdf_set_warning_callback(c, warning);

reload("Wkhtmltox")
reload("VegaLite")

module Try
end

module Try

using VegaLite
using Wkhtmltox

ts = sort(rand(10))
ys = Float64[ rand()*0.1 + cos(x) for x in ts]

v = data_values(time=ts, res=ys) +    # add the data vectors & assign to symbols 'time' and 'res'
      mark_line() +                   # mark type = line
      encoding_x_quant(:time) +       # bind x dimension to :time, quantitative scale
      encoding_y_quant(:res)

tmppath = "/tmp/vegalite.html"
open(io -> VegaLite.writehtml(io, v), tmppath, "w")

png_fn = "/tmp/vegalite.png"

img_init(0)
is = ImgSettings()
is = ImgSettings("in" => tmppath,
                 "out" => png_fn,
                 "fmt" => "png")  # png format output

is["in"]
is["in"] = "http://example.com"
is["out"]

is["fmt"]
is["fmt"]
is["screenHeight"] = "200"
is["screenHeight"]
is["smartWidth"] = "true"
is["smartWidth"]

pars = ["crop.left", "crop.top", "crop.width", "crop.height", "load.cookieJar",
        "transparent", "in", "out", "fmt", "screenWidth", "smartWidth",
        "quality", "web.background", "web.loadImages", "web.enableJavascript",
        "web.enableIntelligentShrinking", "web.minimumFontSize",
        "web.printMediaType", "web.defaultEncoding", "web.userStyleSheet",
        "web.enablePlugins", "load.username", "load.password", "load.jsdelay",
        "load.zoomFactor", "load.customHeaders", "load.repertCustomHeaders",
        "load.cookies", "load.post", "load.blockLocalFileAccess",
        "load.stopSlowScript", "load.debugJavascript", "load.loadErrorHandling",
        "load.proxy", "load.runScript", "header.fontSize", "header.fontName",
        "header.left", "header.center", "header.right", "header.line",
        "header.spacing", "header.htmlUrl"]

is["useGraphics"]

is["left"]

for p in pars
  println("$p => '$(is[p])'")
end

is["load.jsdelay"] = "1000"
is["load.jsdelay"]

is["web.defaultEncoding"]
is["web.defaultEncoding"] = "utf-8"

conv = ImgConverter(is)
run(conv)

out = Array(UInt8, 200_000)

function img_get_output(conv::Ptr{Wkhtmltox.Converter}, output::Vector{UInt8})
  ccall((:wkhtmltoimage_get_output, Wkhtmltox.libwkhtml),
        Int,
        (Ptr{Wkhtmltox.Converter}, Cstring),
        conv, convert(Cstring, pointer(output)))
end

img_get_output(conv.ptr, out)

tmppath = "/tmp/vegalite2.png"
open(io -> write(io, out[1:151930]), tmppath, "w")

out[1:10]
out[0]

[ convert(Int, x) for x in out[1:10]]

len = wkhtmltoimage_get_output(c, &data);
	printf("%ld len\n", len);

conv = nothing
img_deinit()

####################################

ccall((:wkhtmltoimage_extended_qt, Wkhtmltox.libwkhtml),
      Int, ())

img_init(1)
is = ImgSettings("in" => tmppath,
                 "out" => png_fn,
                 "fmt" => "png")  # png format output

conv = ImgConverter(is)
run(conv)
conv = nothing
img_deinit()

#########

reload("Wkhtmltox")
reload("VegaLite")

module Try
end

module Try

using VegaLite
using Wkhtmltox

pdf_init(0)
ps = PdfSettings("out" => "/tmp/example.pdf")   # orientation

ps["size.pageSize"]  # settings can be read like this

os = PdfObject()
# os["page"] = "https://github.com/fredo-dedup/Wkhtmltox.jl"
os["page"] = "/tmp/vegalite.html"
# os["page"] = "/tmp/example.html"
# os["page"] = "http://example.com"
# os["page"] = joinpath(dirname(@__FILE__),"../examples/example.html")

conv = PdfConverter(ps, os) # pass the pdf settings and source settings
run(conv)

conv = nothing
pdf_deinit()

end
#### conversion to image ####