text:
"""
    Web3

A module for Ethereum connectivity.

Parse ABI files, encode/decode ABI data, and make JSON-RPC calls.
"""
module Web3

#=
Parse ABI files and encode / decode ABI data

TUPLE AND FIXED LENGTH ARRAY: [fix1 ... fixN][var1 ... varN]
for static types, fixN is the encoded data and varN is empty
for dynamic types, fixN is the length as a uint256 and varN is the encoded data

DYNAMIC ARRAY: [LEN][ARRAY ENCODING...]

UINT<M> and INT<M>: encoded as 32-byte, big-endian values, i.e. as uint256 and int256 numbers

BYTES: [len as uint256][data]

STRING: a bytes representation of the UTF-8 data

ADDRESS: a uint160, encoded as a 32-byte, big-endian value

BOOL: a uint8, encoded as a 32-byte, big-endian value

BYTES<M>: the sequence of bytes, padded with trailing zeroes to a length of 32-bytes

FUNCTION CALL: [4-byte hash][args tuple]
hash is the first 4 high-order bytes of the hash of the signature name(type,type...)

RETURN VALUE: [return tuple]

EVENT: [contract address][topics[0]: signature][topics[n]: indexed args[n - 1]][data -- unindexed args]
topics are 32-bytes; dynamic values are represented as a hash and lose information
=#

using HTTP, JSON, Core

export readABI, Web3Connection, Contract, contracts
export FunctionCall, encodefunctioncall, decodefunctioncall
export FunctionResult, encodefunctionresult, decodefunctionresult
export Event, encodeevent, decodeevent
export clientversion, eth, utils, net
export computetypes, ContractContext, contract, connection, functions
export sha3_224, sha3_256, sha3_384, sha3_512
export setverbose
#export keccac224, keccac256, keccac384, keccac512

include("keccak.jl")

####################
# Web3
####################

"""
    Web3Connection

A JSON-RPC connection to an Ethereum node

```
    web3 = Web3Connection("http://localhost:8545")
    clientversion(web3) # -- return the client version
    jsonget("http://localhost:8545", :web3_clientVersion, []) # -- equivalent to clientversion(web3)
    rawjsonget("http://localhost:8545", :web3_clientVersion, []) # -- return full JSON object
```
"""
struct Web3Connection
    url
end

"""
    jsonget(url, methodname, params...)

Call a JSON-RPC method and return the result property of the JSON result
"""
#jsonget(url, method, params...) = rawjsonget(url, method, params...)["result"]
function jsonget(url, method, params...)
    json = rawjsonget(url, method, params...)
    if verbose println("JSON: ", repr(json)) end
    json["result"]
end

"""
    rawjsonget(url, method, params...)

Call a JSON-RPC method and return the full JSON result
"""
function rawjsonget(url, method, params...)
    req = JSON.json(Dict([
        :jsonrpc => "2.0"
        :method => method
        :params => params
    ]))
    if verbose println("\nREQUEST: $(req)\n") end
    resp = HTTP.request("POST", url, [], req)
    if resp.status == 200
        resultstr = String(resp.body)
        result = JSON.parse(resultstr)
        if verbose
            println("RESULT JSON: $resultstr")
            println("RESULT: $(repr(result))")
        end
        result
    else
        throw(ErrorException("Error, status = $(resp.status)"))
    end
end

"""
    apifunc

Create a call to an API func, given a function that converts inputs to JSON-ready inputs
"""
apifunc(apimethod, func) = (con::Web3Connection, args...; raw=false)-> (raw ? rawjsonget : jsonget)(con.url, apimethod, func(args...)...)

function hash end

const clientversion = apifunc(:web3_clientVersion, ()-> ())

const net = (
    version = apifunc(:net_version, ()-> ()),
    peercount = apifunc(:net_peerCount, ()-> ()),
    listening = apifunc(:net_listening, ()-> ()),
)

const eth = (
    protocolversion = apifunc(:eth_protocolVersion, ()-> ()),
    syncing = apifunc(:eth_syncing, ()-> ()),
    mining = apifunc(:eth_mining, ()-> ()),
    coinbase = apifunc(:eth_coinbase, ()-> ()),
    hashrate = apifunc(:eth_hashrate, ()-> ()),
    accounts = apifunc(:eth_accounts, ()-> ()),
    gasprice = apifunc(:eth_gasPrice, ()-> ()),
    blocknumber = apifunc(:eth_blockNumber, ()-> ()),
    pendingtransactions = apifunc(:eth_pendingTransactions, ()-> ()),
    newblockfilter = apifunc(:eth_newBlockFilter, ()-> ()),
    newpendingtransactionfilter = apifunc(:eth_newPendingTransactionFilter, ()-> ()),
    getwork = apifunc(:eth_getWork, ()-> ()),
    getbalance = apifunc(:eth_getBalance, (addr, ctx)-> (addr, ctx)),
    getstorageat = apifunc(:eth_getStorageAt, (addr, pos, ctx)-> (addr, pos, ctx)),
    getcode = apifunc(:eth_getCode, (addr, ctx)-> (addr, ctx)),
    getblockbyhash = apifunc(:eth_getBlockByHash, (hash, ctx)-> (hash, ctx)),
    getblockbynumber = apifunc(:eth_getBlockByNumber, (tag, ctx)-> (tag, ctx)),
    getblocktransactioncountbyhash = apifunc(:eth_getBlockTransactionCountByHash, (hash)-> (hash,)),
    getblocktransactioncountbynumber = apifunc(:eth_getBlockTransactionCountByNumber, (tag)-> (tag,)),
    gettransactioncount = apifunc(:eth_getTransactionCount, (addr, ctx)-> (addr, ctx)),
    gettransactionbyhash = apifunc(:eth_getTransactionByHash, (hash)-> (hash,)),
    gettransactionreceipt = apifunc(:eth_getTransactionReceipt, (hash)-> (hash,)),
    estimategas = apifunc(:eth_estimateGas, (dict)-> dict),
    call = apifunc(:eth_call, (dict)-> dict),
    sendtransaction = apifunc(:eth_sendTransaction, (from, to, gas, gasprice, value, data, nonce)->
                              Dict([:from => from
                                    :to => to
                                    :gas => gas
                                    :gasprice => gasprice
                                    :value => value
                                    :data => data
                                    :nonce => nonce])),
    getlogs = apifunc(:eth_getLogs, (dict)-> dict),
)

const db = (
    putstring = apifunc(:db_putString, (db, key, value)-> (db, key, value,)),
    getstring = apifunc(:db_getString, (db, key)-> (value,)),
    puthex = apifunc(:db_putHex, (db, key, value)-> (db, key, value,)),
    gethex = apifunc(:db_getHex, (db, key)-> (value,)),
)

const shh = (
    version = apifunc(:shh_version, ()-> ()),
    newidentity = apifunc(:shh_newIdentity, ()-> ()),
)

const utils = (
    keccak = hash,
)

# This is a separate function so that test code can override it
function hash(con::Web3Connection, str::String)
    #hex2bytes(jsonget(con.url, :web3_sha3, ("0x" * bytes2hex(Vector{UInt8}(str))))[3:end])
    keccak256(collect(UInt8, str))
end

"""
    resultbytes

Convert a JSON-RPC API call result to bytes
"""
resultbytes(func) = (args...)-> hex2bytes(func(args...)[3:end])

####################
# ABI
####################

## Int256 and UInt256 should be changed to primitive types
## See https://github.com/rfourquet/BitIntegers2.jl/blob/master/src/BitIntegers.jl

struct Int256
    big::Int128
    little::UInt128
    Int256(i::Integer) = new(i < 0 ? -1 : 0, i)
    Int256(big::Integer, little::Unsigned) = new(big, little)
end

struct UInt256
    big::UInt128
    little::UInt128
    UInt256(i::Unsigned) = new(0, i)
    UInt256(big::Unsigned, little::Unsigned) = new(big, little)
end

struct ABIType{T}
    ABIType(arg) = new{arg}()
end

struct Decl{ENCMODE, X, Y, Z}
    name::String
    typename::String
    components::Array
    indexed::Bool
end

struct ABIFunction
    constant
    hash # first 4 bytes of keccak hash
    inputs::Array{Decl}
    name
    outputs
    payable
    signature
    statemutability
    argtypes
end

struct ABIEvent
    name
    hash
    signature
    anonymous
    inputs::Array{Decl}
end

struct FunctionCall
    decl
    inputs::Array
end

struct FunctionResult
    decl
    result
end

struct Event
    decl
    parameters::Array
end

struct Contract{Name}
    id::String
    functions
    events
    function Contract(id::String)
        id = cleanaddress(id)
        new{Symbol(id)}(id, Dict{Union{String, Vector{UInt8}}, ABIFunction}(), Dict())
    end
end

function cleanaddress(str::String)
    if match(r"^0[xX]", str) != nothing
        str = str[3:end]
    end
    @assert match(r"[0-9a-fA-F]{20}", str) != nothing
    str
end

const NumDecl = Union{Decl{T, :int} where T, Decl{T, :uint} where T}
const FunctionABI = Union{ABIType{:function}, ABIType{:constructor}, ABIType{:fallback}};

"A dictionary of contract-address => Contract structures"
const contracts = Dict()

const fixedarraypattern = r".*\[([^\]]+)\]"
const bitspattern = r"^([^[0-9]+)([0-9]*)"

##############
# ENCODING
##############

"""
    encodefunctioncall(f::ABIFunction, inputs::Array) -> data
    encodefunctioncall(io::IO, f::ABIFunction, inputs::Array)
    encodefunctioncall(io::IOBuffer, f::ABIFunction, inputs::Array) -> data

Encode a call to a function
"""
encodefunctioncall(f::ABIFunction, inputs::Array) = encodefunctioncall(IOBuffer(), f, inputs)
encodefunctioncall(io::IO, f::ABIFunction, inputs::Array) = basicencodefunctioncall(io, f, inputs)
function encodefunctioncall(io::IOBuffer, f::ABIFunction, inputs::Array)
    basicencodefunctioncall(io, f, inputs)
    take!(io)
end

function basicencodefunctioncall(io::IO, f::ABIFunction, inputs::Array)
    if length(f.inputs) != length(inputs)
        throw("Wrong number of inputs to $(f.name), expecting $(length(f.inputs)) but got $(length(inputs))")
    end
    write(io, f.hash)
    encode(io, f.inputs, inputs)
end

"""
    encodefunctionresult(io::IO, f::ABIFunction, outputs::Array)

Encode the results of a function
"""
function encodefunctionresult(io::IO, f::ABIFunction, outputs::Array)
    basicencodefunctionresult(io, f, outputs)
end

function encodefunctionresult(io::IOBuffer, f::ABIFunction, outputs::Array)
    basicencodefunctionresult(io, f, outputs)
    io.data
end

function basicencodefunctionresult(io::IO, f::ABIFunction, outputs::Array)
    write(io, f.hash)
    encode(io, f.outputs, outputs)
end

"""
    encodeevent(io::IO, e::ABIEvent, inputs::Array)

Encode an event
"""
function encodeevent(io::IO, e::ABIEvent, inputs::Array)
    basicencodeevent(io, e, inputs)
end

function encodeevent(io::IOBuffer, e::ABIEvent, inputs::Array)
    basicencodeevent(io, e, inputs)
    io.data
end

function basicencodeevent(io::IO, e::ABIEvent, inputs::Array)
    write(io, e.hash)
    encode(io, e.inputs, inputs)
end

function encode(io::IO, decl::Union{Decl, Array}, value)
    encodehead(io, decl, value, length(value) * 32)
    encodetail(io, decl, value)
end

encint(io, i::Signed) = encints(Int128(i < 0 ? -1 : 0), Int128(i))
encint(io, i::Unsigned) = encints(io, UInt128(0), UInt128(i))
encint(io, i::Union{Int256, UInt256}) = encints(i.big, i.little)
encints(io, ints...) = vcat(hton.(ints)...)

# utilities
writeint(io, i::Signed) = writeints(io, Int128(i < 0 ? -1 : 0), Int128(i))
writeint(io, i::Unsigned) = writeints(io, UInt128(0), UInt128(i))
writeint(io, i::Union{Int256, UInt256}) = writeints(io, i.big, i.little)
writeints(io, ints...) = write(io, hton.(ints)...)
encodescalar(io::IO, decl::Decl{:scalar, :bool}, value) = writeint(io, value ? 1 : 0)
encodescalar(io::IO, decl::NumDecl, value) = writeint(io, value)

# scalar types
function encodehead(io::IO, decl::Decl{:scalar}, v, offset::Int)
    encodescalar(io, decl, v)
    offset
end
encodetail(io, ::Decl{:scalar}, v) = nothing

# dynamic types
function encodehead(io::IO, ::Union{Decl{:string}, Decl{:bytes}}, v, offset::Int)::Int
    writeint(io, offset)
    offset + ceil(length(IOBuffer(v).data) / 32) * 32 + 32
end
function encodetail(io, ::Union{Decl{:string}, Decl{:bytes}}, v)
    buf = IOBuffer(v)
    len = bytesavailable(buf)
    writeint(io, len)
    write(io, v)
    pad = 32 - len % 32
    if pad != 32
        for i = 1 : pad
            write(io, Int8(0))
        end
    end
end
function encodehead(io::IO, ::Decl{:dynamic}, v, offset::Int)
    writeint(io, offset)
    offset + length(v) * 64
end
function encodetail(io::IO, decl::Decl{:dynamic}, values)
    for v in values
        encodescalar(decl, v)
    end
end

# fixed-size array
function encodehead(io::IO, decl::Decl{:array, BASE, BITS, LENGTH}, values, offset::Int) where {BASE, BITS, LENGTH}
    writeint(io, offset)
    offset + LENGTH * 64
end
function encodetail(io::IO, decl::Decl{:array, BASE, BITS, LENGTH}, values) where {BASE, BITS, LENGTH}
    t = arraycomptype(decl)
    writeint(length)
    offset = 0
    for i in 1:LENGTH
        offset = encodehead(io, t, values[i], offset)
    end
end
function arraycomptype(decl::Decl{:array, BASE, BITS, LENGTH}) where {BASE, BITS, LENGTH}
    if decl.components != nothing
        Decl{:tuple, length(decl.components), :n, :n}(decl.name, "tuple", decl.components, false)
    else
        Decl{:scalar, BASE, BITS, :none}(decl.name, decl.typename, :nothing, false)
    end
end

# tuple
encodehead(io::IO, decl::Decl{:tuple}, values, offset) = encodehead(io, decl.components, values, offset)
encodetail(io::IO, decl::Decl{:tuple}, values) = encodetail(io, decl.components, values)
function encodehead(io::IO, decls::Array{Decl}, values, offset::Int)
    for i in 1:length(decls)
        offset = encodehead(io, decls[i], values[i], offset)
    end
    offset
end
function encodetail(io::IO, decls::Array{Decl}, values)
    for i in 1:length(decls)
        encodetail(io, decls[i], values[i])
    end
end

############
# DECODING
############

const signedTypes = Dict([t.size => t for t in (Int8, Int16, Int32, Int64, Int128)])
const unsignedTypes = Dict([t.size => t for t in (UInt8, UInt16, UInt32, UInt64, UInt128)])

function readint(io::IO)
    big = read(io, Int128)
    Int256(ntoh(big), ntoh(read(io, UInt128)))
end

function readuint(io::IO)
    big = read(io, UInt128)
    Int256(ntoh(big), ntoh(read(io, UInt128)))
end

readlength(io::IO) = (read(io, UInt128);ntoh(read(io, UInt128)))

"""
    decodefunctioncall(io::IO, con::Contract)

Decode a function call
"""
function decodefunctioncall(io::IO, con::Contract)
    decl = con.functions[read(io, 4)]
    FunctionCall(decl, decode(io, decl.inputs))
end

"""
    decodefunctionresult(io::IO, con::Contract)

Decode a function call result
"""
function decodefunctionresult(io::IO, con::Contract)
    hash = read(io, 4)
    f = con.functions[hash]
    FunctionResult(f, decode(io, f.outputs))
end

"""
    decodeevent(io::IO, con::Contract)

Decode an event in a transaction log
"""
function decodeevent(io::IO, con::Contract)
    decl = con.events[read(io, 4)]
    Event(decl, decode(io, decl.inputs))
end

# general
decode(io::IO, decl) = decodetail(io, decl, decodehead(io, decl))

# scalar types
decodehead(io::IO, decl::Decl{:scalar}) = decodescalar(io, decl)
decodetail(io::IO, ::Decl{:scalar}, head) = head
function decodescalar(io::IO, ::Decl{T, :int, SIZE}) where {T, SIZE}
    big = ntoh(read(io, Int128))
    little = ntoh(read(io, UInt128))
    SIZE <= 128 ? smallint(SIZE, big < 0 ? -little : little) : Int256(big, little)
end
function decodescalar(io::IO, ::Decl{T, :uint, SIZE}) where {T, SIZE}
    big = ntoh(read(io, UInt128))
    little = ntoh(read(io, UInt128))
    SIZE <= 128 ? smalluint(SIZE, little) : UInt256(big, little)
end

smallint(size, value) = signedTypes[Int(2^ceil(log2(floor((33 + 7) / 8))))](value)
smalulint(size, value) = unsignedTypes[Int(2^ceil(log2(floor((33 + 7) / 8))))](value)

function decode(io::IO, ::Decl{T, :bool}) where T
    big = read(io, UInt128)
    little = read(io, UInt128)
    litle == 0 ? false : true
end

# dynamic types
decodehead(io::IO, decl::Decl{:dynamic, :bytes}) = readlength(io)
decodetail(io::IO, decl::Decl{:dynamic, :bytes}, head) = readbytes(io, head)
readbytes(io::IO, len) = read(io, len)

function decode(io::IO, decl::Decl{:string})
    if verbose println("Reading string at $(position(io))") end
    offset = readlength(io)
    pos = position(io)
    if verbose println("Offset: $(offset)") end
    seek(io, offset)
    count = readlength(io)
    if verbose println("Length: $(count)") end
    result = String(read(io, count))
    seek(io, pos)
    result
end

# array types
function decode(io::IO, decl::Decl{:array, BASE, BITS, LENGTH}) where {BASE, BITS, LENGTH}
    [decodescalar(io, decl) for i in 1:LENGTH]
end
function decode(io::IO, ::Decl{:array, :bytes, BITS, LENGTH}) where {BITS, LENGTH}
    lens = [readlength(io) for i in 1:LENGTH]
    [readbytes(io, len) for len in lens]
end
function decode(io::IO, ::Decl{:array, :string, BITS, LENGTH}) where {BITS, LENGTH}
    lens = [readlength(io) for i in 1:LENGTH]
    [String(readbytes(io, len)) for len in lens]
end
function decode(io::IO, decl::Decl{:array, :tuple, BITS, LENGTH}) where {BITS, LENGTH}
    decodetail(io::IO, decl.components, decodehead(io::IO, decl.components))
end

# tuple
decode(io::IO, decl::Decl{:tuple}) = decodetail(io, decl.components)
#decode(io::IO, decls::Array) = decodetail(io, decls, decodehead(io, decls))
#decodehead(io::IO, decls::Array) = [decodehead(io, head) for head in decls]
#decodetail(io::IO, decls::Array, heads) = [decodetail(io, decls[i], heads[i]) for i in 1:length(decls)]
function decode(io::IO, decls::Array)
    result = []
    for decl in decls
        push!(result, decode(io, decl))
    end
    result
end

####################
# DECL PARSING
####################

rows(array) = [array[row, :] for row in 1:size(a)[1]]

# Make conversion mapping given specs:
#   (Solidity-prefix, Julia type, byte length)
conversions(rows) = vcat([["$stype$(bytes * 8)" => jtype for bytes in rng] for (stype, jtype, rng) in rows]...)

const soliditytojulia = Dict(conversions([
    # Solidity-prefix, Julia type, bytes
    ("int", Int8, [8])
    ("uint", UInt8, [8])
    ("int", Int16, [16])
    ("uint", UInt16, [16])
    ("int", Int32, 3:4)
    ("uint", UInt32, 3:4)
    ("int", Int64, 5:8)
    ("uint", UInt64, 5:8)
    ("int", Int128, 9:16)
    ("uint", UInt128, 9:16)
    ("int", BigInt, 17:32)
    ("uint", BigInt, 17:32)
]))

parseABI(connection::Web3Connection, json) = parseABI(connection, ABIType(Symbol(get(json, "type", "function"))), json)
function parseABI(connection::Web3Connection, ::FunctionABI, func)
    name = func["name"]
    args = join((arg-> arg["type"]).(func["inputs"]), ",")
    sig = "$name($args)"
    inputs = parseargs(func["inputs"])
    ABIFunction(
        func["constant"],
        utils.keccak(connection, sig)[1:4],
        inputs,
        name,
        haskey(func, "outputs") ? parseargs(func["outputs"]) : [],
        func["payable"],
        sig,
        func["stateMutability"],
        computetypes(name, inputs)
    )
end

function parseABI(connection::Web3Connection, ::ABIType{:event}, evt)
    name = evt["name"]
    args = join((arg-> arg["type"]).(evt["inputs"]), ",")
    sig = "$name($args)"
    ABIEvent(evt["name"], utils.keccak(connection, sig)[1:4], sig, get(evt, "anonymous", false), parseargs(evt["inputs"]))
end

computetype(decl::Decl{T, :int, SIZE}) where {T, SIZE} = soliditytojulia["int$SIZE"]
computetype(decl::Decl{T, :uint, SIZE}) where {T, SIZE} = soliditytojulia["uint$SIZE"]
computetype(decl::Decl{:tuple, SIZE}) where {T, SIZE} = "tuple[$SIZE]"
computetype(decl::Decl{:array, BASE, BITS, LENGTH}) where {BASE, BITS, LENGTH} = "array[LENGTH] of $BASE"
computetype(decl::Decl{T, :bool} where {T}) = Bool
computetype(decl::Decl{:string}) = String
computetype(decl::Decl{:bytes}) = Vector{UInt8}
computetype(decl::Decl{:dynamic, TYPE}) where TYPE = "array of $TYPE"

computetypes(func::ABIFunction) = computetypes(func.name, func.inputs)
function computetypes(name, decls::Array{T} where T <: Decl)
    computetype.(decls)
end

"""
    readABI(con::Web3Connection, contractname::String, stream::IO)

Read an ABI file for a contract
"""
function readABI(connection::Web3Connection, contractname::String, stream::IO)
    contract = Contract(contractname)
    d = JSON.parse(stream)
    close(stream)
    for json in (haskey(d, "abi") ? d["abi"] : d)
        obj = parseABI(connection, json)
        if isa(obj, ABIFunction)
            contract.functions[obj.name] = contract.functions[obj.hash] = obj
            if verbose
                println("$(bytes2hex(obj.hash)) $(obj.signature) $(repr(obj)) $(bytes2hex(utils.keccak(connection, obj.signature)))")
            end
        elseif isa(obj, ABIEvent)
            contract.events[obj.name] = contract.events[obj.hash] = obj
            if verbose
                println("$(bytes2hex(obj.hash)) $(obj.signature) $(repr(obj)) $(bytes2hex(utils.keccak(connection, obj.signature)))")
            end
        elseif verbose
            println(repr(obj))
        end
    end
    contracts[contractname] = contract
end

function basetypefor(typename)
    m = match(bitspattern, typename)
    (Symbol(m[1]), m[2] == "" ? 256 : parse(Int, m[2]))
end

function typefor(typename, arg)
    atype = arg["type"]
    if endswith(atype, "[]")
        (:dynamic, basetypefor(typename)..., :none)
    elseif endswith(atype, "]")
        m = match(fixedarraypattern, typename)
        (:array, basetypefor(typename)..., parse(Int, m[1]))
    elseif atype == "tuple"
        (:tuple, length(arg["components"]), :none, :none)
    elseif atype in ["string", "bytes"]
        (Symbol(atype), :none, :none, :none)
    else
        (:scalar, basetypefor(typename)..., :none)
    end
end

function parsearg(arg)
    typename = arg["type"]
    (enctype, atype, bits, len) = typefor(typename, arg)
    Decl{enctype, atype, bits, len}(
        arg["name"],
        typename,
        haskey(arg, "components") ? parseargs(arg["components"]) : [],
        get(arg, "indexed", false)
    )
end

parseargs(args) = parsearg.(args)

####################
# UTILS
####################

global verbose = false

struct ContractContext{contractid}
    connection
    contract
end

connection(con::ContractContext) = getfield(con, :connection)

contract(con::ContractContext) = getfield(con, :contract)

functions(con::ContractContext) = contract(con).functions

struct Val{Name} end

function gen(contract, con)
    funcs = collect(filter(p-> isa(p[1], String), contract.functions))
    funcdict = Dict(map(funcs) do ((name, func)) name => func end)
    methods = map(funcs) do ((name, func))
        argnames = map(a-> Symbol("a$a"), 1:length(func.argtypes))
        args = map(((name, type)::Tuple)-> :($name::$type), zip(argnames, func.argtypes))
        method = :($(Symbol(name)) = (
            send = ($(args...),; options...)-> send(context, $name, [$(argnames...)]; options...),
            call = ($(args...),; options...)-> call(context, $name, [$(argnames...)]; options...),
            estimategas = ($(args...),; options...)-> (),
            encodeabi = ($(args...),; options...)-> begin
                buf = IOBuffer()
                encodefunctioncall(buf, contract(context).functions[$(String(name))], [$(argnames...)])
                take!(buf)
            end
        ))
        if verbose println("\n$name = ", method, "\n") end
        method
    end
    type = ContractContext{Symbol(contract.id)}
    eval(:(Base.getproperty(context::$type, prop::Symbol) = ($(methods...),)[prop]))
end

"""
    send(id, data; options)

options are: gasprice, gas, value, nonce, chain, hardfork, common
"""
function send(context::ContractContext, name, args; options...)
    func = contract(context).functions[name]
    extractresult(func, basicsend(:eth_sendTransaction, context, name, args; options...))
end

function call(context::ContractContext, name, args; options...)
    func = contract(context).functions[name]
    extractresult(func, basicsend(:eth_call, context, name, args; options...))
end

function extractresult(func::ABIFunction, json::Dict)
    if haskey(json, "result")
        if verbose println("DECODING RESULT: ", json["result"]) end
        result = decode(IOBuffer(hex2bytes(json["result"][3:end])), func.outputs)
        if length(result) == 1
            result[1]
        else
            (result...,)
        end
    elseif haskey(json, "error")
        throw(json["error"])
    else
        throw("Unknown result: $(repr(json))")
    end
end

function basicsend(op, context::ContractContext, name, args; options...)
    #println("Call $name in contract $(contract(context).id)\nfrom account $(from)")
    transaction = merge!(Dict([
        :to => lowercase("0x" * contract(context).id)
        :data=> "0x" * bytes2hex(encodefunctioncall(contract(context).functions[name], args)) #encodearguments(args)
        :id => 1
    ]),
                         pairs(options))
    if verbose println(op, " $name DATA: ", transaction[:data]) end
    result = if op == :eth_sendTransaction
        rawjsonget(connection(context).url, op, transaction)
    else
        rawjsonget(connection(context).url, op, transaction, "latest")
    end
    if haskey(result, "error")
        err = result["error"]
        if verbose println("ERROR: $(err["message"])\n$(err["data"]["stack"])") end
    end
    if verbose println(op, " $name RESULT JSON: $(repr(result))") end
    result
end

encodearguments(args) = if args == () "" else :ENCODED_ARGUMENTS end

function ContractContext(url::String, contractid::String, jsonabifile::String)
    ContractContext(Web3Connection(url), contractid, jsonabifile)
end
ContractContext(con::Web3Connection, contractid::String, jsonabifile::String) = ContractContext(con, contractid, open(jsonabifile))

function ContractContext(con::Web3Connection, contractid::String, file::IO)
    contract = readABI(con, contractid, file)
    gen(contract, con)
    ContractContext{Symbol(contract.id)}(con, contract)
end

function setverbose(v)
    global verbose
    verbose = v
end

end # module


text:
# runtests.jl

using Test: @test, @testset, @test_throws, @inferred, detect_ambiguities
using Sound # sound, soundsc, record, pick_output
using SampledSignals: SampleBuf
using PortAudio: devices, PortAudioDevice

# include("aqua.jl") # fails in Interval, Unitful ...

@show devices()

if isempty(devices())
    @warn "No devices so no tests on this OS."
else

@testset "output" begin
    dev = @inferred Sound.get_default_output_device()
    @test dev isa PortAudioDevice

    index = findfirst(==(dev), devices())
    io_out = IOBuffer()
    io_in = IOBuffer("$index")
    tmp = @inferred Sound.pick_output( ; io_in, io_out)
    @test tmp === dev

    io_in = IOBuffer("\n")
    tmp = @inferred Sound.pick_output( ; io_in, io_out)
    @test tmp === dev

    io_in = IOBuffer("$index")
    sound(:pick, randn(500), 1000; io_in, io_out)
    sound(index, randn(500), 1000)
    @test_throws String sound(0, randn(500), 1000)
end

@testset "Sound" begin
    S = 8192 # sampling rate in Hz
    x = 0.8 * cos.(2π*(1:800)*440/S)
    y = 0.7 * sin.(2π*(1:800)*660/S)
    sound(x, S) # specify sampling rate
    @test_throws ErrorException sound(y) # unknown sampling rate
    sound([x y], S) # stereo
    soundsc(0.1 * reshape(x,1,:), S) # scale and row vector

    sb = SampleBuf(x, S)
    sound(sb)
    soundsc(sb)

    @test_throws String sound(ones(5,3), S) # array size

    @test isempty(detect_ambiguities(Sound))
end

@testset "record" begin
    dev = @inferred Sound.get_default_input_device()
    @test dev isa PortAudioDevice

    data, S = record(0.001)
    @test S isa Real
    @test data isa Vector
end

end # !isempty(devices())

include("hann.jl")
include("findpeaks.jl")
include("phase_vocoder.jl")

ext: jl


text:
# (31+42)::Float64 #> TypeError(:typeassert,"",Float64,73)
# convert(Int64, 7.01) #> ERROR: InexactError()
convert(Int64, 7.0) #> 7
int64(7.0) #> 7
int64(7.01) #> 7
# convert(Int64, "CV")
# ERROR: `convert` has no method matching convert(::Type{Int64}, ::ASCIIString)
# convert(Int64, "123") # same ERROR as above

promote(1, 2.5, 3//4) #> (1.0,2.5,0.75)
promote(1.5, im) #> (1.5 + 0.0im,0.0 + 1.0im)
promote(true, 'c', 1.0) #> (1.0,99.0,1.0)

promote_type(Int8, Uint16) #> Int64

max_stars_repo_path: Module 1/Chapter06/conversions.jl

max_issues_repo_path: Module 1/Chapter06/conversions.jl

max_forks_repo_path: Module 1/Chapter06/conversions.jl


text:
## wiener filter

@doc doc"""
    argmaxf_lnP(ϕ,                ds::DataSet; kwargs...)
    argmaxf_lnP(ϕ, θ::NamedTuple, ds::DataSet; kwargs...)
    argmaxf_lnP(Lϕ::LenseOp,      ds::DataSet; kwargs...)

Computes either the Wiener filter at fixed $\phi$, or a sample from this slice
along the posterior.

Keyword arguments:

* `which` — `:wf`, `:sample`, or `fluctuation` to compute 1) the Wiener filter,
  i.e. the best-fit of $\mathcal{P}(f\,|\,\phi,d)$, 2) a sample from
  $\mathcal{P}(f\,|\,\phi,d)$, or 3) a sample minus the Wiener filter, i.e. the
  fluctuation on top of the mean.
* `guess` — starting guess for `f` for the conjugate gradient solver
* `kwargs...` — all other arguments are passed to [`conjugate_gradient`](@ref)

"""
argmaxf_lnP(ϕ::Field,                ds; kwargs...) = argmaxf_lnP(cache(ds.L(ϕ),ds.d), ds();      kwargs...)
argmaxf_lnP(ϕ::Field, θ::NamedTuple, ds; kwargs...) = argmaxf_lnP(cache(ds.L(ϕ),ds.d), ds(;θ...); kwargs...)

function argmaxf_lnP(Lϕ::LenseOp, ds::DataSet; which=:wf, guess=nothing, kwargs...)

    check_hat_operators(ds)
    @unpack d, Cn, Cn̂, Cf, M, M̂, B, B̂, P = ds

    b = 0
    if (which in (:wf, :sample))
        b += Lϕ'*B'*P'*M'*(Cn\d)
    end
    if (which in (:fluctuation, :sample))
        b += Cf\simulate(Cf) + Lϕ'*B'*P'*M'*(Cn\simulate(Cn))
    end

    conjugate_gradient(
        pinv(Cf) + M̂'*B̂'*pinv(Cn̂)*B̂*M̂,
        pinv(Cf) + Lϕ'*B'*P'*M'*pinv(Cn)*M*P*B*Lϕ,
        b,
        guess==nothing ? 0*b : guess;
        kwargs...
    )

end

@doc doc"""
    Σ(ϕ, ds, ::Type{L}=LenseFlow) where {L}
    Σ(L::LenseOp, ds)

An operator for the data covariance, Cn + P*M*B*L*Cf*L'*B'*M'*P', which can
applied and inverted.
"""
Σ(ϕ, ds) = Σ(ds.L(ϕ),ds)
Σ(Lϕ::LenseOp, ds) = begin

    @unpack d,P,M,B,Cn,Cf,Cn̂,B̂ = ds

    SymmetricFuncOp(
        op   = x -> (Cn + P*M*B*Lϕ*Cf*Lϕ'*B'*M'*P')*x,
        op⁻¹ = x -> conjugate_gradient((Cn̂ .+ B̂*Cf*B̂'), Σ(Lϕ, ds), x, nsteps=100, tol=1e-1)
    )

end

@doc doc"""

    MAP_joint(ds::DataSet; kwargs...)

Compute the maximum a posteriori (i.e. "MAP") estimate of the joint posterior,
$\mathcal{P}(f,\phi,\theta\,|\,d)$, or compute a quasi-sample.

Keyword arguments:

* `ϕstart` — Starting point of the maximizer *(default:* $\phi=0$*)*
* `Nϕ` — Noise to use in the approximate hessian matrix. Can also give `Nϕ=:qe`
         to use the EB quadratic estimate noise *(default:* `:qe`*)*
* `quasi_sample` — `true` to iterate quasi-samples, or an integer to compute
                   a specific quasi-sample.
* `nsteps` — The number of iterations for the maximizer
* `Ncg` — Maximum number of conjugate gradient steps during the $f$ update
* `cgtol` — Conjugrate gradient tolerance (will stop at `cgtol` or `Ncg`, whichever is first)
* `αtol` — Absolute tolerance on $\alpha$ in the linesearch in the $\phi$ quasi-Newton-Rhapson step, $x^\prime = x - \alpha H^{-1} g$
* `αmax` — Maximum value for $\alpha$ in the linesearch
* `progress` — `false`, `:summary`, or `:verbose`, to control progress output

Returns a tuple `(f, ϕ, tr)` where `f` is the best-fit (or quasi-sample) field,
`ϕ` is the lensing potential, and `tr` contains info about the run.

"""
function MAP_joint(
    ds;
    ϕstart = nothing,
    Nϕ = nothing,
    quasi_sample = false,
    nsteps = 10,
    Ncg = 500,
    cgtol = 1e-1,
    αtol = 1e-5,
    αmax = 0.5,
    cache_function = nothing,
    callback = nothing,
    interruptable = false,
    progress = false)

    @assert progress in [false,:summary,:verbose]
    if !(isa(quasi_sample,Bool) || isa(quasi_sample,Int))
        throw(ArgumentError("quasi_sample should be true, false, or an Int."))
    end

    # since MAP estimate is done at fixed θ, we don't need to reparametrize to
    # ϕₘ = G(θ)*ϕ, so set G to constant here to avoid wasted computation
    @set! ds.G = IdentityOp
    @unpack d, D, Cϕ, Cf, Cf̃, Cn, Cn̂, L = ds

    f, f° = nothing, nothing
    ϕ = (ϕstart==nothing) ? zero(diag(Cϕ)) : ϕstart
    Lϕ = cache(L(ϕ),d)
    T = real(eltype(d))
    α = 0
    tr = []
    hist = nothing

    # compute approximate inverse ϕ Hessian used in gradient descent, possibly
    # from quadratic estimate
    if (Nϕ == :qe); Nϕ = quadratic_estimate(ds).Nϕ/2; end
    Hϕ⁻¹ = (Nϕ == nothing) ? Cϕ : pinv(pinv(Cϕ) + pinv(Nϕ))

    try
        @showprogress (progress==:summary ? 1 : Inf) "MAP_joint: " for i=1:nsteps

            # ==== f step ====

            # if we're doing a fixed quasi_sample, set the random seed here,
            # which controls the sample from the posterior we get from inside
            # `argmaxf_lnP`
            if isa(quasi_sample,Int); seed!(quasi_sample); end

            # recache Lϕ for new ϕ
            if i!=1; cache!(Lϕ,ϕ); end

            # run wiener filter
            (f, hist) = argmaxf_lnP(((i==1 && ϕstart==nothing) ? NoLensing() : Lϕ), ds,
                    which = (quasi_sample==false) ? :wf : :sample, # if doing a quasi-sample, we get a sample instead of the WF
                    guess = (i==1 ? nothing : f), # after first iteration, use the previous f as starting point
                    tol=cgtol, nsteps=Ncg, hist=(:i,:res), progress=(progress==:verbose))

            f° = Lϕ * D * f
            lnPcur = lnP(:mix,f°,ϕ,ds)

            # ==== show progress ====
            if (progress==:verbose)
                @printf("(step=%i, χ²=%.2f, Ncg=%i%s)\n", i, -2lnPcur, length(hist), (α==0 ? "" : @sprintf(", α=%.6f",α)))
            end
            push!(tr,@namedtuple(i,lnPcur,hist,ϕ,f,α))
            if callback != nothing
                callback(f, ϕ, tr)
            end

            # ==== ϕ step =====
            if (i!=nsteps)
                ϕnew = Hϕ⁻¹*(@ondemand(Zygote.gradient)(ϕ->lnP(:mix,f°,ϕ,ds), ϕ)[1])
                res = @ondemand(Optim.optimize)(α->(-lnP(:mix,f°,ϕ+α*ϕnew,ds)), T(0), T(αmax), abs_tol=αtol)
                α = res.minimizer
                ϕ = ϕ+α*ϕnew
            end

        end
    catch err
        if (err isa InterruptException) && interruptable
            println()
            @warn("Maximization interrupted. Returning current progress.")
        else
            rethrow(err)
        end
    end

    return @namedtuple(f, ϕ, tr)

end

@doc doc"""

    MAP_marg(ds; kwargs...)

Compute the maximum a posteriori (i.e. "MAP") estimate of the marginal posterior,
$\mathcal{P}(\phi,\theta\,|\,d)$.

"""
function MAP_marg(
    ds;
    ϕstart = nothing,
    L = LenseFlow,
    Nϕ = nothing,
    nsteps = 10,
    Ncg = 500,
    cgtol = 1e-1,
    α = 0.02,
    Nmc_det = 50,
    )

    @unpack Cf, Cϕ, Cf̃, Cn̂ = ds

    # compute approximate inverse ϕ Hessian used in gradient descent, possibly
    # from quadratic estimate
    if (Nϕ == :qe); Nϕ = quadratic_estimate(ds).Nϕ/2; end
    Hϕ⁻¹ = (Nϕ == nothing) ? Cϕ : pinv(pinv(Cϕ) + pinv(Nϕ))

    ϕ = (ϕstart != nothing) ? ϕstart : ϕ = zero(diag(Cϕ))
    tr = []

    for i=1:nsteps
        g, det_sims = δlnP_δϕ(ϕ, ds, progress=true, Nmc_det=Nmc_det, return_sims=true)
        ϕ += α * Hϕ⁻¹ * g
        push!(tr,@dictpack(i,g,det_sims,ϕ))
    end

    return ϕ, tr

end


text:
module Tokamak

using MacroTools, DataFlow
using DataFlow: TypeAssert, Lambda, OLambda, prewalk, postwalk, isconstant, value,
  inputnode, constant, iscall, spliceinputs, syntax

export @tk

include("utils.jl")
include("syntax.jl")
include("loops.jl")
include("types.jl")

include("codegen/cpu.jl")
include("codegen/dagger.jl")

end # module

max_stars_repo_name: JuliaTagBot/Tokamak

max_issues_repo_name: JuliaTagBot/Tokamak

max_forks_repo_name: FluxML/Tokamak


text:
function quaternion_map(ω, timestep)
    return Quaternion(sqrt(4 / timestep^2 - dot(ω, ω)), ω)
end

function quaternion_map_jacobian(ω::SVector{3}, timestep)
    msq = -sqrt(4 / timestep^2 - dot(ω, ω))
    return [ω' / msq; I]
end

function cayley(ω)
    Quaternion(1.0 / sqrt(1.0 + ω' * ω) * [1.0; ω])
end

function cayley_jacobian(ω)
    ω₁, ω₂, ω₃ = ω
    a = (sqrt(1.0 + ω₁^2 + ω₂^2 + ω₃^2)^-3)
    b = sqrt(1.0 + ω₁^2 + ω₂^2 + ω₃^2)^-1
    SMatrix{4,3}([
                 -ω₁*a -ω₂*a -ω₃*a;
                 (b - (ω₁^2)*a) (-ω₁ * ω₂ * a) (-ω₁ * ω₃ * a);
                 (-ω₁ * ω₂ * a) (b - (ω₂^2)*a) (-ω₂ * ω₃ * a);
                 (-ω₁ * ω₃ * a) (-ω₂ * ω₃ * a) (b - (ω₃^2)*a);
                 ])
end


text:
# Method:
# 1. van-e közös merőleges? nincs -> break
# 2. összes pont levetítése erre a síkra (egyik pont és a közös normális)
# 3. legkisebb és legnagyobb távolság megnézése
# 4. AABB
# 5. ha az AABB területe nagyon kicsi -> break
# 6. legnagyobb összefüggő terület
# 7. kontúr kiszedése: kell-e, hogy zárt görbe legyen? - szerintem kell -> 2 végpont összekötése
# 8. kör/egyenes illesztése
# 9. visszaellenőrzés?

abstract type AbstractTranslationalSurface <: RANSAC.FittedShape end

struct FittedTranslational <: AbstractTranslationalSurface
    coordframe
    contourindexes
    subsetnum::Int
end

struct ExtractedTranslational <: AbstractTranslationalSurface
    coordframe
    contour
    # center of gravity
    center
    # normal of the contour is parallel to the direction
    # towards the center of the contour?
    # == should flip the computed normals to direct outwards?
    # this is used in e.g. CSGBuilding
    # true means, that the computed normals must be turned to direct outside
    outwards::Int
    # should the computed normal be flipped to match the measured points
    # this is used in this package to ensure that in/outwards is correct
    # true means that computed normals must be turned to match the measured points
    flipnormal::Int
    # for visualization
    ft::FittedTranslational
end

function RANSAC.defaultshapeparameters(::Type{FittedTranslational})
    #ϵ = 0.3
    #α = deg2rad(5)
    # max deviaton of the normal being perpendicular to the translation direction
    #α_perpend = cosd(89)
    # maximum number of contours on a plane
    #TODO: delete this
    #max_group_num::Int = 3
    # maximum number of iterations of tryíng
    # to find < max_group_num number for contour patches
    #TODO: delete this
    #max_contour_it::Int = 5
    #thinning_par = 2.0
    # minimum % of the normals must be the same
    #min_normal_num = 0.9
    # extract translational surface even though normals are not ok
    #force_transl::Bool = false
    # thinning method: :slow/:fast/:deldir
    #thin_method::Symbol = :slow
    # how close must they be to consider them as the same point?
    #samep = Float64(eps(Float32))
    # check side parameter
    #checksidepar = 0.04
    # "disabled by default"
    #max_end_d = 10000.0

    tr = (ϵ = 0.3, α = deg2rad(5), α_perpend = cosd(89),
        diagthr = 0.1, max_group_num= 3, max_contour_it= 5, thinning_par = 2.0,
        min_normal_num = 0.9, force_transl = false, thin_method= :slow,
        samep = Float64(eps(Float32)), checksidepar = 0.04, max_end_d = 10000.0,)
    return (translational=tr,)
end

RANSAC.defaultshapeparameters(::Type{ExtractedTranslational}) = RANSAC.defaultshapeparameters(FittedTranslational)

Base.show(io::IO, x::FittedTranslational) = print(io, """FittedTranslational""")
Base.show(io::IO, x::ExtractedTranslational) = print(io, """ExtractedTranslational""")

Base.show(io::IO, ::MIME"text/plain", x::FittedTranslational) =
    print(io, """FittedTranslational""")

Base.show(io::IO, ::MIME"text/plain", x::ExtractedTranslational) =
    print(io, """ExtractedTranslational""")

RANSAC.strt(x::AbstractTranslationalSurface) = "translational"

RANSACVisualizer.getcolour(s::AbstractTranslationalSurface) = get(colorschemes[:seaborn_bright], 0.4)

function transldir(p, n, params)
    # 1. van-e közös merőleges? nincs -> break
    #@unpack α_perpend = params
    @extract params : myp=translational
    @extract myp : α_perpend
    is = [[1,2, 3], [2,3, 1], [3,1, 2]]
    ds = [abs(dot(n[i[1]], n[i[2]])) for i in is]
    sel_i = is[argmin(ds)]
    direction = normalize(cross(n[sel_i[1]], n[sel_i[2]]))
    isapprox(norm(direction), 1) || return (false, direction)
    abs(dot(direction, n[sel_i[3]])) < α_perpend && return (true, direction)
    return (false, direction)
end

"""
    project2sketchplane(pcr, indexes, transl_frame, params)

Project points defined by `indexes` (`pcr.vertices[indexes]`) to the plane defined by `transl_frame`.
All the points should be enabled.
Only those points are considered, whose normal is perpendicular to the plane.
"""
function project2sketchplane(pcr, indexes, transl_frame, params)
    #@unpack α_perpend = params
    @extract params : myp=translational
    @extract myp : α_perpend
    xv = transl_frame[1]
    yv = transl_frame[2]
    z = transl_frame[3]
    # enabled points and normals
    p = @view pcr.vertices[indexes]
    n = @view pcr.normals[indexes]

    projected = Array{SVector{2,Float64},1}(undef, 0)
    inds = Int[]
    for i in eachindex(p)
        # point is not part of the translational surface
        abs(dot(z, n[i])) > α_perpend && continue
        nv1 = dot(xv, p[i])
        nv2 = dot(yv, p[i])
        nvv = SVector{2, Float64}(nv1, nv2)
        push!(projected, nvv)
        push!(inds, indexes[i])
    end
    return projected, inds
end

"""
    project2sketchplane(points, transl_frame)

Just project the points to the plane.
"""
function project2sketchplane(points, transl_frame)
    xv = transl_frame[1]
    yv = transl_frame[2]
    projd = [SVector{2,Float64}(dot(xv, p), dot(yv, p)) for p in points]
    return projd
end

"""
    segmentpatches(points, ϵ_inrange)

Return connected patches of a pointcloud
"""
function segmentpatches(points, ϵ_inrange)
    btree = KDTree(points)
    uf = UnionFinder(size(points,1))

    for i in eachindex(points)
        inr = inrange(btree, points[i], ϵ_inrange)
        for j in eachindex(inr)
            # the point itself is always in range
            i == inr[j] && continue
            union!(uf, i, inr[j])
        end
    end
    return CompressedFinder(uf)
end

"""
    filtermultipoint!(points, indexes, params)

Filter out points that are close (params.samep) to each other.
Duplicates are removed inplace from the `points` and also their index from `indexes`.
"""
function filtermultipoint!(points, indexes, params)
    #@unpack samep = params
    @extract params : myp=translational
    @extract myp : samep
    btree = KDTree(points)
    trm = Int[]

    for i in eachindex(points)
        inr = inrange(btree, points[i], samep)
        for j in eachindex(inr)
            # i - current index
            # inr[j] - index of a point that is < eps at i
            # i itself is also in inr
            inr[j] <= i && continue
            push!(trm, inr[j])
        end
    end
    sort!(trm)
    unique!(trm)
    @logmsg IterLow1 "Deleted $(length(trm)) duplicates"
    deleteat!(indexes, trm)
    deleteat!(points, trm)
    return indexes
end

"""
    normaldirs(segments, points, normals, center, params)

Compute the normal directions.
This is based on the points, which are closer than `ϵ`.
Return a boolean first that indicates that the normals direct to the "same direction".
"""
function normaldirs(segments, points, normals, center, params)
    @assert size(points) == size(normals)
    function fff(msg)
        if ! (msg === nothing)
            @logmsg IterLow1 msg
        end
        return (false, false, false)
    end
    #@unpack ϵ, min_normal_num = params
    @extract params : myp=translational
    @extract myp : ϵ min_normal_num
    #calcs = [dist2segment(p, segments) for p in points]
    calcs = [dn2contour(p, segments) for p in points]
    compats = [abs(calcs[i][1]) < ϵ for i in eachindex(calcs)]
    compatsize = count(compats)
    compatsize == 0 && return fff("Compat size is 0.")
    # later working with points[compats]
    psize = size(points,1)
    flipnormal = Vector{Bool}(undef, psize)
    outwards = Vector{Bool}(undef, psize)
    for i in 1:psize
        # continue if not compatible
        compats[i] || continue
        # this is the fitted normal
        #contour_n = segmentnormal(segments, calcs[i][2])
        contour_n = calcs[i][2]
        #tocenter = normalize(center-midpoint(segments, calcs[i][3]))
        tocenter = normalize(center-segments[calcs[i][3]])

        # normal of the contour is parallel to the direction
        # towards the center of the contour?
        # == should flip the computed normals to direct outwards?
        # this is used in e.g. CSGBuilding
        outwards[i] = dot(contour_n, tocenter) < 0.0 ? false : true

        # should the computed normal be flipped to match the measured points
        # this is used in this package to ensure that in/outwards is correct
        flipnormal[i] = dot(contour_n, normals[i]) < 0.0 ? true : false
    end
    thisoutw = @view outwards[compats]
    outwr = count(thisoutw)/compatsize
    # can't agree on outwards
    (outwr > min_normal_num) || (outwr <= 1-min_normal_num) || return fff("Bad outwards.")

    # can't agree on flipnormals
    thisflip = @view flipnormal[compats]
    flipr = count(thisflip)/compatsize
    (flipr > min_normal_num) || (flipr <= 1-min_normal_num) || return fff("Bad flipsign.")
    # this means, that the computed normals must be turned to direct outside
    outwb = outwr > min_normal_num ? -1 : 1

    # this means that computed normals must be turned to match the measured points
    flipn = flipr > min_normal_num ? -1 : 1

    return (true, outwb, flipn)
end

"""
    checksides(points, multipl)

`multipl=0.02` for example.
"""
function checksides(points, params)
    #@unpack checksidepar = params
    @extract params : myp=translational
    @extract myp : checksidepar

    obb = RANSAC.findOBB_(points)
    sl1 = norm(obb[1]-obb[2])
    sl2 = norm(obb[1]-obb[3])
    sl3 = norm(obb[1]-obb[5])
    sls = [sl1, sl2, sl3]
    for i in 1:3
        for j in 1:3
            i == j && continue
            sls[i] < checksidepar*sls[j] && return (false, sls)
        end
    end
    return (true, sls)
end

function retnot(msg)
    if ! (msg === nothing)
        @logmsg IterLow1 msg
    end
    return nothing
end

function RANSAC.fit(::Type{FittedTranslational}, p, n, pcr, params)
    #@unpack α_perpend, diagthr, max_group_num = params
    #@unpack max_contour_it, thinning_par, ϵ, τ = params

    @extract params.iteration : τ
    @extract params : myp=translational
    @extract myp : α_perpend diagthr max_group_num
    @extract myp : max_contour_it thinning_par ϵ
    # Method:
    # 1. van-e közös merőleges? nincs -> break
    ok, dir = transldir(p, n, params)
    ok || return retnot("No translational direction.")
    # 2. összes pont levetítése erre a síkra (egyik pont és a közös normális)
    o = p[1]
    xv = n[1]
    zv = dir
    yv = normalize(cross(zv, xv))
    coordframe = [xv, yv, zv]

    #TODO: okosabb kéne ide
    # these are the enabled points from the first subset
    ien = pcr.isenabled
    subsnum = 1
    sbs = pcr.subsets[subsnum]
    # index in isenabled with the subset, to get those who are enabled in the subset
    # then use it to index into the subset
    used_i = @view sbs[ien[sbs]]
    projected, proj_ind = project2sketchplane(pcr, used_i, coordframe, params)
    size(projected, 1) < 2 && return retnot("No compatible points to transl. direction.")

    spatchs = segmentpatches(projected, ϵ)
    # hereby spatchs should contain maximum max_group_num of patches
    fitresults = Array{FittedTranslational,1}(undef, 0)
    @logmsg IterLow1 "Nof groups: $(spatchs.groups)"
    for i in 1:spatchs.groups
        # i if part of the i-th group
        cur_group = findall(x->x==i, spatchs.ids)
        # at least 3 points please...
        size(cur_group,1) < 3 && continue
        patch_indexes = proj_ind[cur_group]
        # only extract contours that contain at least the minimumsize number of points
        size(patch_indexes, 1) < τ && continue
        #original:
        #size(patch_indexes, 1) < τ/size(pcr.subsets,1) && continue
        #@logmsg IterLow1 "Nof contour points: $(length(patch_indexes))"

        # 4. OOBB
        # don't extract planes
        ppp = @view pcr.vertices[patch_indexes]
        goodside, sidel = checksides(ppp, params)
        if !goodside
            @logmsg IterLow1 "One side of OOBB is small: $sidel"
            continue
        end
        ft = FittedTranslational(coordframe, patch_indexes, subsnum)
        push!(fitresults, ft)
    end
    return fitresults
end

## scoring

function RANSAC.scorecandidate(pc, candidate::FittedTranslational, subsetID, params)
    inpoints = candidate.contourindexes
    subsID = candidate.subsetnum
    score = RANSAC.estimatescore(length(pc.subsets[subsID]), pc.size, length(inpoints))
    return (score, inpoints)
end

## refit

function compatiblesTranslational(shape, points, normals, params)
    @extract params : myp=translational
    @extract myp : α  ϵ

    # project to plane
    ps = project2sketchplane(points, shape.coordframe)
    ns = project2sketchplane(normals, shape.coordframe)

    calcs = (dn2shape_contour(p, shape) for p in ps)

    zpn = zip(calcs, ns)
    c = [(abs(cc[1]) < ϵ) && RANSAC.isparallel(cc[2], ni, α) for (cc,ni) in zpn]
    return c
end

"""
    refit(s::FittedTranslational, pc, params)

Refit translational.
"""
function RANSAC.refit(s::FittedTranslational, pc, params)
    @extract params : myp=translational
    @extract myp : ϵ force_transl thin_method
    @extract myp : thinning_par max_end_d

    cf = s.coordframe
    cidxs = s.contourindexes
    # 1. project points & normals
    old_p = @view pc.vertices[cidxs]
    o_pp = project2sketchplane(old_p, cf)

    filtermultipoint!(o_pp, cidxs, params)
    # this uses the filtered cidxs
    old_n = @view pc.normals[cidxs]
    o_np = project2sketchplane(old_n, cf)

    @logmsg IterLow1 "Thinning"
    # 2. thinning
    if thin_method === :fast
        thinned, _ = thinning(o_pp, thinning_par)
    elseif thin_method === :slow
        thinned, _ = thinning_slow(o_pp, thinning_par)
    else
        error("thin method: $(thin_method) is not implemented!")
    end
    closed = [SVector{2,Float64}(th) for th in thinned]
    if norm(closed[1]-closed[end]) > max_end_d
        return retnot("max_end_d: $(norm(closed[1]-closed[end]))")
    end
    c = centroid(closed)
    @logmsg IterLow1 "Normaldirs"
    # 3. normaldirs()
    isok, outw, flips = normaldirs(closed, o_pp, o_np, c, params)
    (isok || force_transl) || return retnot("Normals not ok in refit.")
    et = ExtractedTranslational(cf, closed, c, outw, flips, s)
    @logmsg IterLow1 "Compatibles"
    # 4. search for all enabled and compatibel points
    # TODO: use octree for that
    p = @view pc.vertices[pc.isenabled]
    n = @view pc.normals[pc.isenabled]
    cp = compatiblesTranslational(et, p, n, params)
    ip = ((1:pc.size)[pc.isenabled])[cp]
    return RANSAC.ExtractedShape(et, ip)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module Foo
include("SubFoo1.jl")
include(joinpath("subdir", "SubFoo2.jl"))
import Bar, Baz, Qux
this = "Foo1"
which = "path"
end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
module Automa

using Printf: @sprintf
import ScanByte

include("sdict.jl")
include("sset.jl")

# TODO: use StableDict and StableSet only where they are required
const Dict = StableDict
const Set = StableSet

include("byteset.jl")
include("re.jl")
include("precond.jl")
include("action.jl")
include("edge.jl")
include("nfa.jl")
include("dfa.jl")
include("machine.jl")
include("traverser.jl")
include("dot.jl")
include("memory.jl")
include("codegen.jl")
include("tokenizer.jl")
include("Stream.jl")

end # module


text:
using DataFrames
using Gadfly

# Input parameters
path = "/scratch/cns7ae/ciris/trackplot_final17/"
abundanceFile = "trackplot.csv"
plotname = "trackplot.pdf"
imax = 428.
jmax = 428.
kmax = 428.
isize = 1.3e-5*1e7 #nm
jsize = 1.3e-5*1e7 #cm
ksize = 1.3e-5*1e7 #cm
xtrack=(237.0/jmax)*(jsize)
ytrack=(249.0/kmax)*(ksize)

plottheme = Gadfly.Theme(
                        default_point_size=1.5pt,
                        line_width=2pt,
                        minor_label_font_size=13pt,
                        minor_label_color=colorant"black",
                        major_label_font_size=15pt,
                        major_label_color=colorant"black",
                        key_label_color=colorant"black",
                        key_title_color=colorant"black",
                        key_title_font_size=15pt,
                        key_label_font_size=10pt,
                        grid_color=colorant"black",
                        highlight_width=0pt)

println("Reading input file")
track = readtable(
                  path*abundanceFile,
                  names = [:z, :y, :x, :Species, :Action],
                  eltypes = [Float64, Float64, Float64, String, String],
                  header = false
                 )

for i in 1:nrow(track)
  track[:z][i] = ((track[:z][i]/imax)*isize)
  track[:y][i] = ((track[:y][i]/jmax)*jsize)
  track[:x][i] = ((track[:x][i]/kmax)*ksize)
end

track = track[isize.>track[:z].>0,:]
track = track[jsize.>track[:y].>0,:]
track = track[ksize.>track[:x].>0,:]

# trackview_df = DataFrame()
# trackview_df = vcat(trackview_df,
#                     DataFrame(
#                               x=track[:j],
#                               y=track[:i],
#                               Species=track[:Species],
#                               view="View 1"
#                               ),
#                     DataFrame(
#                               x=track[:k],
#                               y=track[:i],
#                               Species=track[:Species],
#                               view="View 2"
#                               ),
#                     DataFrame(
#                               x=track[:k],
#                               y=track[:j],
#                               Species=track[:Species],
#                               view="View 3"
#                               )
#                     )

dens_df = DataFrame()
dens_df = vcat(dens_df,
               DataFrame(
                         x=track[:x],
                         y=track[:y],
                         r=0.0,
                         Species=track[:Species]
                         )
              )

for i in 1:length(dens_df[:x])
  dens_df[:r][i] = sqrt((dens_df[:x][i]-xtrack)^2 + (dens_df[:y][i]-ytrack)^2)
end

#Sort the columns by r-value
sort!(dens_df, cols = [order(:r)])
#Find the radius containing 2σ of the damage
P = 95.45
tot = nrow(dens_df)
ind = round(Int,((P*tot)/100))
radius = dens_df[:r][ind]

function circletop(x,h,k,r)
  sqrt((r^2 - (x-h)^2)) + k
end

function circlebottom(x,h,k,r)
   -1*sqrt((r^2 - (x-h)^2)) + k
end

#Generate the circle that will be plotted over the top-down graph
xstart = xtrack-radius
xend = xtrack+radius
xpoint = xstart:1e-2:xend
ypoint = Array(Float32,size(xpoint,1))
ypoint2 = Array(Float32,size(xpoint,1))
for i in 1:size(xpoint,1)
  ytemp = 0.0
  xtemp = xpoint[i]
  ytemp = circletop(xtemp,xtrack,ytrack,radius)
  ypoint[i] = ytemp
  ypoint2[i] = circlebottom(xtemp,xtrack,ytrack,radius)
end

println("Now plotting data")

pij = plot(
           track[100.>track[:z].>1,:],
           x="y",
           y="z",
           Coord.cartesian(aspect_ratio=1),
           Geom.point,
           Guide.xlabel("y (nm)"),
           Guide.ylabel("z (nm)"),
           Guide.manual_color_key(
                                   "Species",
                                   ["Proton","Electron"],
                                   ["green","deepskyblue"],
                                  ),
           Scale.x_continuous(minvalue=50, maxvalue=100),
           Scale.y_continuous(minvalue=0, maxvalue=50),
           Scale.color_discrete_manual(colorant"magenta",colorant"green",colorant"deepskyblue")
           )
push!(pij,plottheme)

proton_track = DataFrame()
proton_track[:z] = linspace(1,100,100)
proton_track[:y] = 71.986
proton_track[:x] = 75.6308
proton_track[:Species] = "proton"

push!(pij,layer(proton_track[proton_track[:Species].=="proton",:],x="y",y="z",Geom.line,Theme(default_color=colorant"green",line_style=Gadfly.get_stroke_vector(:dash),line_width=2pt)))

pik = plot(
           track[100.>track[:z].>1,:],
           x="x",
           y="z",
           Geom.point,
           Guide.manual_color_key(
                                   "Species",
                                   ["Proton","Electron"],
                                   ["green","deepskyblue"],
                                  ),
           Guide.xlabel("x (nm)"),
           Guide.ylabel("z (nm)"),
           Coord.cartesian(aspect_ratio=1),
           Scale.x_continuous(minvalue=50, maxvalue=100),
           Scale.y_continuous(minvalue=0, maxvalue=50),
           Scale.color_discrete_manual(colorant"magenta",colorant"green",colorant"deepskyblue")
          #  Scale.color_discrete_manual(colorant"deepskyblue",colorant"magenta",colorant"green")
          )

push!(pik,plottheme)
push!(pik,layer(proton_track[proton_track[:Species].=="proton",:],x="x",y="z",Geom.line,Theme(default_color=colorant"green",line_style=Gadfly.get_stroke_vector(:dash),line_width=2pt)))

pjk_path = plot(
                track,
                Guide.xlabel("x (nm)"),
                Guide.ylabel("y (nm)"),
                Coord.cartesian(aspect_ratio=1.0),
                Scale.x_continuous(minvalue=50, maxvalue=100),
                Scale.y_continuous(minvalue=50, maxvalue=100),
                Scale.color_discrete_manual(colorant"magenta",colorant"green",colorant"deepskyblue"),
           Guide.manual_color_key(
                                   "Species",
                                   ["Proton","Electron"],
                                   ["green","deepskyblue"],
                                  ),
                layer(
                      x=[75.63],
                      y=[71.986],
                      Geom.point,
                      Theme(default_color=colorant"green",
                             default_point_size=3pt,
                             highlight_width=0pt)
                      ),
                layer(
                      x="x",
                      y="y",
                      Geom.point,
                      ),

                )
push!(pjk_path,plottheme)
#push!(pjk_path,layer(track[track[:Species].=="proton",:],x="x",y="y",Geom.point,Theme(default_color=colorant"green",default_point_size=5pt,highlight_width=0pt)))

pjk = plot(
#           dens_df[50.>dens_df[:x].>0,:],
           dens_df,
           x="x",
           y="y",
#           color="Species",
#           Geom.point,
#           Geom.contour,
           Guide.xlabel("x (nm)"),
           Guide.ylabel("y (nm)"),
           Coord.cartesian(fixed=true),
           Geom.histogram2d(xbincount=55, ybincount=55),
#           Scale.x_continuous(minvalue=0, maxvalue=50),
#           Scale.y_continuous(minvalue=0, maxvalue=50),
            Scale.color_log10(colormap=Scale.lab_gradient(
                                                          colorant"darkblue",
                                                          colorant"mediumblue",
                                                          colorant"deepskyblue",
                                                          colorant"turquoise",
                                                          colorant"lightgreen",
                                                          colorant"gold",
                                                          colorant"yellow"
                                                          ))
            )
push!(pjk,plottheme)

circle_df = DataFrame()
circle_df = vcat(circle_df,
                  DataFrame(
                            x=xpoint,
                            y=ypoint,
                            Action="radius= 18.2 nm"),
                  )

#push!(pjk,layer(circle_df,x="x",y="y",Geom.point,Theme(default_color=colorant"red")))
push!(pjk,layer(circle_df,x="x",y="y",Geom.path,Theme(default_color=colorant"red")))
push!(pjk,layer(x=xpoint,y=ypoint2,Geom.path,Theme(default_color=colorant"red")))

#set_default_plot_size(50cm, 7.5cm)
# p_tot = plot(trackview_df,
#              layer(
#                    trackview_df[trackview_df[:Species] .== "electron", :],
#                    Geom.subplot_grid(Geom.point),
#                    xgroup="view", x="x", y="y",
#                    color="Species",
#                    Theme(default_point_size=1.5pt),
#                    ),
#              layer(
#                    trackview_df[trackview_df[:Species] .== "proton", :],
# #                   Geom.subplot_grid(Geom.point),
#                    Geom.point,
#                    xgroup="view", x="x", y="y",
#                    color="Species",
#                    Theme(default_point_size=1.5pt),
#                    )
#
#             )
#

#push!(p_tot, layer(trackview_df[trackview_df[:Species] .== "low-energy electron", :], xgroup="view", x="x", y="y", color="Species",
#             Geom.subplot_grid(Geom.point)))

println("Now exporting plot")
draw(PDF(path*"f3a.pdf", 6inch, 5inch), pij)
draw(PDF(path*"f3b.pdf", 6inch, 5inch), pik)
draw(PDF(path*"f3c.pdf", 6inch, 5inch), pjk_path)
draw(PDF(path*"f3d.pdf", 6inch, 5inch), pjk)
println("Now ending script")


text:
function catch_blow_up_time(ode; time_stepper=RDPK3SpFSAL49(), abstol=1.0e-7, reltol=1.0e-7,
                            callbacks=AliveCallback(alive_interval = 100))
    integrator = init(ode, time_stepper, abstol=abstol, reltol=reltol,
                      save_everystep=false, callback=callbacks)
    try
        solve!(integrator)
    catch e
        @warn "Blow-up at time " e
    end
    return integrator.t
end

reinterpret_DGSEM_solution_as_DGMulti(u, semi) =
    reinterpret_DGSEM_solution_as_DGMulti(u, Trixi.mesh_equations_solver_cache(semi)...)

function reinterpret_DGSEM_solution_as_DGMulti(u, mesh::TreeMesh{NDIMS}, equations, dg::DGSEM, cache) where {NDIMS}
    nvars = nvariables(equations)
    uEltype = eltype(u)
    reshaped_u_raw = reinterpret(reshape, SVector{nvars, uEltype},
                                 reshape(u, nvars, nnodes(dg)^NDIMS, nelements(mesh, dg, cache)))
    reshaped_u = StructArray{SVector{nvars, uEltype}}(ntuple(_ -> zeros(uEltype, size(reshaped_u_raw)), nvars))
    reshaped_u .= reshaped_u_raw

    return reshaped_u
end

function integrate_entropy(sol::Trixi.TrixiODESolution)
    S = zeros(length(sol.t))
    for i in eachindex(sol.t)
        S[i] = integrate_entropy(sol.u[i], sol.prob.p)
    end
    return S
end

integrate_entropy(u, semi) = integrate_entropy(Trixi.wrap_array(u, semi), Trixi.mesh_equations_solver_cache(semi)...)

function integrate_entropy(u, mesh, equations, dg::DGMulti, cache)
    rd = dg.basis

    # Note: this evaluates at the inherent quadrature used by the DGMulti solver,
    # which may not be consistent with Lobatto DGSEM quadrature. However, we have
    # observed that using Lobatto quadrature can result in negative density/pressure
    # even when the Gauss collocation solver does not crash (since it uses different
    # points).
    @unpack Vq = rd
    @unpack wJq = mesh.md

    uq = similar(u, rd.Nq, size(u, 2))
    Trixi.apply_to_each_field(Trixi.mul_by!(Vq), uq, u)
    S = sum(wJq .* entropy.(uq, equations))
    return S
end

using LinearAlgebra: diagm

function integrate_entropy(u, mesh, equations, solver::DGSEM, cache)
    nvars = nvariables(equations)
    reshaped_u = reinterpret(reshape, SVector{nvars, Float64}, u)
    J = 1.0 / nelements(mesh, solver, cache)
    w1D = solver.basis.weights
    wq = vec([w1D[i] * w1D[j] for j in eachnode(solver), i in eachnode(solver)])
    wJq = diagm(wq) * J * ones(nnodes(solver)^2, nelements(mesh, solver, cache))
    S = sum(reshape(wJq, size(reshaped_u)) .* entropy.(reshaped_u, equations))
    return S
end

function integrate_entropy(u, mesh::TreeMesh{3}, equations, solver::DGSEM, cache)
    nvars = nvariables(equations)
    reshaped_u = reinterpret(reshape, SVector{nvars, Float64}, u)
    J = 1.0 / nelements(mesh, solver, cache)
    w1D = solver.basis.weights
    wq = vec([w1D[i] * w1D[j] * w1D[k] for k in eachnode(solver), j in eachnode(solver), i in eachnode(solver)])
    wJq = diagm(wq) * J * ones(nnodes(solver)^3, nelements(mesh, solver, cache))
    S = sum(reshape(wJq, size(reshaped_u)) .* entropy.(reshaped_u, equations))
    return S
end

function compute_entropy_over_time(sol::ODESolution)
    semi = sol.prob.p
    S = zeros(length(sol.u))
    for (i, u) in enumerate(sol.u)
      S[i] = compute_entropy_over_time(u, semi)
    end
    return S
end

function compute_entropy_over_time(u::AbstractVector, semi)
    reshaped_u = reinterpret_DGSEM_solution_as_DGMulti(u, semi)
    return compute_entropy_over_time(reshaped_u, semi)
end

function compute_entropy_over_time(u::AbstractMatrix, semi)
    if typeof(u) <: StructArray
        rho = StructArrays.component(u, 1)
    else
        rho = getindex.(u, 1)
    end
    @unpack equations = semi
    p = pressure.(u, equations)
    S = integrate_entropy(u, Trixi.mesh_equations_solver_cache(semi)...)
    return S
end

## Routines which depend on dimension

# converts a solution array `u` with size (nnodes, cells_per_dimension^2) to a logically
# Cartesian array. Assumes the mesh is Cartesian.
function convert_dg2D_array_to_Cartesian(x, y, u, polydeg, cells_per_dimension)
    h = 1 / cells_per_dimension
    cell_centroids = LinRange(-1 + h, 1 - h, cells_per_dimension)
    my_mean(x) = sum(x)/length(x)

    num_pts_per_element = polydeg + 1
    Nfft = cells_per_dimension * num_pts_per_element
    u_cartesian = zeros(eltype(u), Nfft, Nfft)
    for e in 1:size(u, 2)
        xc = my_mean(view(x, :, e))
        yc = my_mean(view(y, :, e))
        u_e = reshape(view(u, :, e), polydeg + 1, polydeg + 1)'

        tol = 1e2*eps()
        ex = findfirst(@. abs(cell_centroids - xc) < tol)
        ey = findfirst(@. abs(cell_centroids - yc) < tol)

        row_ids = (1:num_pts_per_element) .+ (ey - 1) * num_pts_per_element
        col_ids = (1:num_pts_per_element) .+ (ex - 1) * num_pts_per_element
        u_cartesian[row_ids, col_ids] .= u_e
    end
    return u_cartesian
end

# Given an ODESolution, returns Cartesian arrays representing the primitive solution fields.
convert_to_Cartesian_arrays(sol) =
    convert_to_Cartesian_arrays(sol.u[end], Trixi.mesh_equations_solver_cache(sol.prob.p)...)

convert_to_Cartesian_arrays(u, semi) =
    convert_to_Cartesian_arrays(u, Trixi.mesh_equations_solver_cache(semi)...)

function convert_to_Cartesian_arrays(u, mesh, equations::CompressibleEulerEquations2D, dg::DGMulti, cache)
    rd = dg.basis
    polydeg = rd.N

    # This scales the node positions to create an equispaced grid of nodes which avoids
    # boundary evaluations. The solution sampled at these points can be used to compute
    # the FFT and the turbulent energy spectra
    equispaced_fft_nodes = map(x -> x * (1 - 1 / (polydeg + 1)), StartUpDG.equi_nodes(Quad(), polydeg))
    interp_matrix_fft_nodes = StartUpDG.vandermonde(Quad(), polydeg, equispaced_fft_nodes...) / rd.VDM

    # interpolate to equispaced FFT nodes
    primitive_variables = StructArrays.components(cons2prim.(u, equations))
    rho, u, v, p = map(u -> interp_matrix_fft_nodes * u, primitive_variables)
    x, y = map(u -> interp_matrix_fft_nodes * u, mesh.md.xyz)

    cells_per_dimension = Int(sqrt(mesh.md.num_elements)) # assume uniform Cartesian mesh
    rho_cartesian = convert_dg2D_array_to_Cartesian(x, y, rho, polydeg, cells_per_dimension)
    u_cartesian   = convert_dg2D_array_to_Cartesian(x, y, u, polydeg, cells_per_dimension)
    v_cartesian   = convert_dg2D_array_to_Cartesian(x, y, v, polydeg, cells_per_dimension)
    p_cartesian   = convert_dg2D_array_to_Cartesian(x, y, p, polydeg, cells_per_dimension)
    return rho_cartesian, u_cartesian, v_cartesian, p_cartesian, x, y
end

function convert_to_Cartesian_arrays(u, mesh, equations::CompressibleEulerEquations2D, dg::DGSEM, cache)
    polydeg = Trixi.polydeg(dg)

    # This scales the node positions to create an equispaced grid of nodes which avoids
    # boundary evaluations. The solution sampled at these points can be used to compute
    # the FFT and the turbulent energy spectra
    h = 1 / (polydeg + 1)
    equispaced_fft_nodes_1D = LinRange(-1 + h, 1 - h, polydeg + 1)
    interp_matrix_fft_1D = Trixi.polynomial_interpolation_matrix(dg.basis.nodes, equispaced_fft_nodes_1D)
    interp_matrix_fft_nodes = kron(interp_matrix_fft_1D, interp_matrix_fft_1D)

    # interpolate solution to FFT points
    u_solution = reinterpret_DGSEM_solution_as_DGMulti(u, mesh, equations, dg, cache)
    primitive_variables = StructArrays.components(cons2prim.(u_solution, equations))
    rho, u, v, p = map(u -> interp_matrix_fft_nodes * u, primitive_variables)

    # reshape arrays to have the same format as DGMulti
    x_coordinates = view(cache.elements.node_coordinates, 1, :, :, :)
    y_coordinates = view(cache.elements.node_coordinates, 2, :, :, :)
    x = reshape(x_coordinates, nnodes(dg)^2, nelements(mesh, dg, cache))
    y = reshape(y_coordinates, nnodes(dg)^2, nelements(mesh, dg, cache))
    x, y = map(u -> interp_matrix_fft_nodes * u, (x, y))

    cells_per_dimension = Int(sqrt(nelements(mesh, dg, cache))) # assume uniform Cartesian mesh
    rho_cartesian = convert_dg2D_array_to_Cartesian(x, y, rho, polydeg, cells_per_dimension)
    u_cartesian   = convert_dg2D_array_to_Cartesian(x, y, u, polydeg, cells_per_dimension)
    v_cartesian   = convert_dg2D_array_to_Cartesian(x, y, v, polydeg, cells_per_dimension)
    p_cartesian   = convert_dg2D_array_to_Cartesian(x, y, p, polydeg, cells_per_dimension)
    return rho_cartesian, u_cartesian, v_cartesian, p_cartesian, x, y
end

function compute_2D_energy_spectrum(rho_cartesian, u_cartesian, v_cartesian, p_cartesian)
    uhat = fft(sqrt.(rho_cartesian) .* u_cartesian)
    vhat = fft(sqrt.(rho_cartesian) .* v_cartesian)
    Ek_full = @. 0.5 * (abs(uhat)^2 + abs(vhat)^2)

    Ek = Ek_full[1:end÷2+1, 1:end÷2+1] # remove duplicate modes

    wavenumbers = 1:size(Ek, 1)
    effective_wavenumbers = @. sqrt(wavenumbers^2 + wavenumbers'^2)

    N = length(wavenumbers)
    Ek_1D = zeros(N * (N - 1) ÷ 2 + N) # number of unique wavenumbers = triangular part + diagonal
    effective_wavenumber_1D = zeros(N * (N - 1) ÷ 2 + N)
    sk = 1
    for i in wavenumbers
        for j in wavenumbers[i:end] # use that the effective_wavenumber matrix is symmetric
            kk = sqrt(i^2 + j^2)
            ids = findall(@. kk - 0.5 < effective_wavenumbers < kk + 0.5) # find wavenumbers in LinRange
            effective_wavenumber_1D[sk] = kk
            Ek_1D[sk] = sum(Ek[ids])
            sk += 1
        end
    end
    return Ek_1D, effective_wavenumber_1D
end


text:
# Unit test functions
function unit_test_slinkedlist_empty()
    l = SLinkedList{Float64}()
    @test l.head.value === nothing
    @test l.tail.value === nothing
    @test length(l) == 0

    return nothing
end

function unit_test_slinkedlist_addlast()
    l = SLinkedList{Float64}()
    val1 = 1.0
    addlast!(l, val1)
    @test l.head.value == val1
    @test l.tail.value == val1
    @test length(l) == 1

    val2 = 2.0
    addlast!(l, val2)
    @test l.head.value == val1
    @test l.head.next.value == val2
    @test l.tail.value == val2
    @test length(l) == 2

    return nothing
end

function unit_test_slinkedlist_addfirst()
    l = SLinkedList{Float64}()
    val1 = 1.0
    addfirst!(l, val1)
    @test l.head.value == val1
    @test l.tail.value == val1
    @test length(l) == 1

    val2 = 2.0
    addfirst!(l, val2)
    @test l.head.value == val2
    @test l.head.next.value == val1
    @test l.tail.value == val1
    @test length(l) == 2

    return nothing
end

function unit_test_slinkedlist_add()
    l = SLinkedList{Float64}()
    val1 = 1.0
    add!(l, val1, 0)
    @test l.head.value == val1
    @test l.tail.value == val1
    @test length(l) == 1
    # 1

    val2 = 2.0
    add!(l, val2, 0)
    @test l.head.value == val2
    @test l.head.next.value == val1
    @test l.tail.value == val1
    @test length(l) == 2
    # 2 -> 1

    val3 = 3.0
    add!(l, val3, 5)
    @test l.head.value == val2
    @test l.head.next.value == val1
    @test l.head.next.next.value == val3
    @test l.tail.value == val3
    @test length(l) == 3
    # 2 -> 1 -> 3

    val4 = 4.0
    add!(l, val4, 2)
    @test l.head.value == val2
    @test l.head.next.value == val4
    @test l.head.next.next.value == val1
    @test l.head.next.next.next.value == val3
    @test l.tail.value == val3
    @test length(l) == 4
    # 2 -> 4 -> 1 -> 3

    val5 = 5.0
    add!(l, val5, 3)
    @test l.head.value == val2
    @test l.head.next.value == val4
    @test l.head.next.next.value == val5
    @test l.head.next.next.next.value == val1
    @test l.head.next.next.next.next.value == val3
    @test l.tail.value == val3
    @test length(l) == 5
    # 2 -> 4 -> 5 -> 1 -> 3

    return nothing
end

function unit_test_slinkedlist_add_interface()
    l = SLinkedList{Float64}()
    val1 = 1.0
    add!(l, val1, 0)
    @test l.head.value == val1
    @test l.tail.value == val1
    @test length(l) == 1
    # 1

    val2 = 2.0
    add!(l, val2)
    @test l.head.value == val1
    @test l.head.next.value == val2
    @test l.tail.value == val2
    @test length(l) == 2
    # 1 -> 2
end

function unit_test_slinkedlist_addlastfirst()
    l = SLinkedList{Float64}()
    val1 = 1.0
    addfirst!(l, val1)
    @test l.head.value == val1
    @test l.tail.value == val1
    @test length(l) == 1

    val2 = 2.0
    addlast!(l, val2)
    @test l.head.value == val1
    @test l.head.next.value == val2
    @test l.tail.value == val2
    @test length(l) == 2

    val3 = 3.0
    addfirst!(l, val3)
    @test l.head.value == val3
    @test l.head.next.value == val1
    @test l.head.next.next.value == val2
    @test l.tail.value == val2
    @test length(l) == 3

    return nothing
end

function unit_test_slinkedlist_removelast()
    l = SLinkedList{Float64}()
    val1 = 1.0
    val2 = 2.0
    val3 = 3.0
    addlast!(l, val1)
    addlast!(l, val2)
    addlast!(l, val3)
    @test l.head.value == val1
    @test l.head.next.value == val2
    @test l.head.next.next.value == val3
    @test l.tail.value == val3
    @test length(l) == 3
    @test removelast!(l) == val3
    @test l.head.value == val1
    @test l.head.next.value == val2
    @test l.tail.value == val2
    @test length(l) == 2
    @test removelast!(l) == val2
    @test l.head.value == val1
    @test l.tail.value == val1
    @test length(l) == 1
    @test removelast!(l) == val1
    @test l.head == l.tail
    @test length(l) == 0

    return nothing
end

function unit_test_slinkedlist_removefirst()
    l = SLinkedList{Float64}()
    val1 = 1.0
    val2 = 2.0
    val3 = 3.0
    addlast!(l, val1)
    addlast!(l, val2)
    addlast!(l, val3)
    @test l.head.value == val1
    @test l.head.next.value == val2
    @test l.head.next.next.value == val3
    @test l.tail.value == val3
    @test length(l) == 3
    @test removefirst!(l) == val1
    @test l.head.value == val2
    @test l.head.next.value == val3
    @test l.tail.value == val3
    @test length(l) == 2
    @test removefirst!(l) == val2
    @test l.head.value == val3
    @test l.tail.value == val3
    @test length(l) == 1
    @test removefirst!(l) == val3
    @test l.head == l.tail
    @test length(l) == 0

    return nothing
end

function unit_test_slinkedlist_remove()
    l = SLinkedList{Float64}()
    val1 = 1.0
    val2 = 2.0
    val3 = 3.0
    addlast!(l, val1)
    addlast!(l, val2)
    addlast!(l, val3)
    @test l.head.value == val1
    @test l.head.next.value == val2
    @test l.head.next.next.value == val3
    @test l.tail.value == val3
    @test length(l) == 3
    @test remove!(l, 2) == val2
    @test l.head.value == val1
    @test l.head.next.value == val3
    @test l.tail.value == val3
    @test length(l) == 2
    @test remove!(l, 0) == val1
    @test l.head.value == val3
    @test l.tail.value == val3
    @test length(l) == 1
    @test remove!(l, 2) == val3
    @test l.head == l.tail
    @test length(l) == 0

    return nothing
end

function unit_test_slinkedlist_remove_interface()
    l = SLinkedList{Float64}()
    val1 = 1.0
    val2 = 2.0
    addlast!(l, val1)
    addlast!(l, val2)
    @test l.head.value == val1
    @test l.head.next.value == val2
    @test l.tail.value == val2
    @test length(l) == 2
    @test remove!(l) == val1
    @test l.head.value == val2
    @test l.tail.value == val2
    @test length(l) == 1
    @test remove!(l) == val2
    @test l.head == l.tail
    @test length(l) == 0

    return nothing
end

function unit_test_slinkedlist_removelastfirst()
    l = SLinkedList{Float64}()
    val1 = 1.0
    val2 = 2.0
    val3 = 3.0
    addlast!(l, val1)
    addlast!(l, val2)
    addlast!(l, val3)
    @test l.head.value == val1
    @test l.head.next.value == val2
    @test l.head.next.next.value == val3
    @test l.tail.value == val3
    @test length(l) == 3
    @test removelast!(l) == val3
    @test l.head.value == val1
    @test l.head.next.value == val2
    @test l.tail.value == val2
    @test length(l) == 2
    @test removefirst!(l) == val1
    @test l.head.value == val2
    @test l.tail.value == val2
    @test length(l) == 1
    @test removelast!(l) == val2
    @test l.head == l.tail
    @test length(l) == 0

    return nothing
end

function unit_test_slinkedlist_peekfirst()
    val1 = 1.0
    val2 = 2.0
    l = SLinkedList{Float64}()
    addlast!(l, val1)
    @test peekfirst(l) == val1
    addlast!(l, val2)
    @test peekfirst(l) == val1

    return nothing
end

function unit_test_slinkedlist_peeklast()
    val1 = 1.0
    val2 = 2.0
    l = SLinkedList{Float64}()
    addlast!(l, val1)
    @test peeklast(l) == val1
    addlast!(l, val2)
    @test peeklast(l) == val2

    return nothing
end

function unit_test_dlinkedlist_empty()
    l = DLinkedList{Float64}()
    @test l.head.value === nothing
    @test l.tail.value === nothing
    @test length(l) == 0

    return nothing
end

function unit_test_dlinkedlist_addlast()
    l = DLinkedList{Float64}()
    val1 = 1.0
    addlast!(l, val1)
    @test l.head.value == val1
    @test l.head.next.prev.value == val1
    @test l.tail.value == val1
    @test l.tail.next.prev.value == val1
    @test length(l) == 1

    val2 = 2.0
    addlast!(l, val2)
    @test l.head.value == val1
    @test l.head.next.value == val2
    @test l.head.next.next.prev.value == val2
    @test l.head.next.prev.value == val1
    @test l.tail.value == val2
    @test l.tail.next.prev.value == val2
    @test length(l) == 2

    return nothing
end

function unit_test_dlinkedlist_addfirst()
    l = DLinkedList{Float64}()
    val1 = 1.0
    addfirst!(l, val1)
    @test l.head.value == val1
    @test l.head.next.prev.value == val1
    @test l.tail.value == val1
    @test l.tail.next.prev.value == val1

    val2 = 2.0
    addfirst!(l, val2)
    @test l.head.value == val2
    @test l.head.next.value == val1
    @test l.head.next.next.prev.value == val1
    @test l.head.next.prev.value == val2
    @test l.tail.value == val1
    @test l.tail.next.prev.value == val1

    return nothing
end

function unit_test_dlinkedlist_add()
    l = DLinkedList{Float64}()
    val1 = 1.0
    add!(l, val1, 0)
    @test l.head.value == val1
    @test l.head.next.prev.value == val1
    @test l.tail.value == val1
    @test l.tail.next.prev.value == val1
    # 1

    val2 = 2.0
    add!(l, val2, 0)
    @test l.head.value == val2
    @test l.head.next.value == val1
    @test l.head.next.next.prev.value == val1
    @test l.head.next.prev.value == val2
    @test l.tail.value == val1
    @test l.tail.next.prev.value == val1
    # 2 <-> 1

    val3 = 3.0
    add!(l, val3, 5)
    @test l.head.value == val2
    @test l.head.next.value == val1
    @test l.head.next.next.value == val3
    @test l.head.next.next.next.prev.value == val3
    @test l.head.next.next.prev.value == val1
    @test l.head.next.prev.value == val2
    @test l.tail.value == val3
    @test l.tail.next.prev.value == val3
    # 2 <-> 1 <-> 3

    val4 = 4.0
    add!(l, val4, 2)
    @test l.head.value == val2
    @test l.head.next.value == val4
    @test l.head.next.next.value == val1
    @test l.head.next.next.next.value == val3
    @test l.head.next.next.next.next.prev.value == val3
    @test l.head.next.next.next.prev.value == val1
    @test l.head.next.next.prev.value == val4
    @test l.head.next.prev.value == val2
    @test l.tail.value == val3
    @test l.tail.next.prev.value == val3
    # 2 <-> 4 <-> 1 <-> 3

    val5 = 5.0
    add!(l, val5, 3)
    @test l.head.value == val2
    @test l.head.next.value == val4
    @test l.head.next.next.value == val5
    @test l.head.next.next.next.value == val1
    @test l.head.next.next.next.next.value == val3
    @test l.head.next.next.next.next.next.prev.value == val3
    @test l.head.next.next.next.next.prev.value == val1
    @test l.head.next.next.next.prev.value == val5
    @test l.head.next.next.prev.value == val4
    @test l.head.next.prev.value == val2
    @test l.tail.value == val3
    @test l.tail.next.prev.value == val3
    # 2 <-> 4 <-> 5 <-> 1 <-> 3

    return nothing
end

function unit_test_dlinkedlist_add_interface()
    l = DLinkedList{Float64}()
    val1 = 1.0
    add!(l, val1, 0)
    @test l.head.value == val1
    @test l.head.next.prev.value == val1
    @test l.tail.value == val1
    @test l.tail.next.prev.value == val1
    # 1

    val2 = 2.0
    add!(l, val2)
    @test l.head.value == val1
    @test l.head.next.value == val2
    @test l.head.next.next.prev.value == val2
    @test l.head.next.prev.value == val1
    @test l.tail.value == val2
    @test l.tail.next.prev.value == val2
    # 1 <-> 2

    return nothing
end

function unit_test_dlinkedlist_addlastfirst()
    l = DLinkedList{Float64}()
    val1 = 1.0
    addlast!(l, val1)
    @test l.head.value == val1
    @test l.tail.value == val1
    @test l.head.next.prev.value == val1
    @test l.tail.next.prev.value == val1
    @test length(l) == 1

    val2 = 2.0
    addfirst!(l, val2)
    @test l.head.value == val2
    @test l.head.next.value == val1
    @test l.head.next.next.prev.value == val1
    @test l.head.next.prev.value == val2
    @test l.tail.value == val1
    @test l.tail.next.prev.value == val1
    @test length(l) == 2

    val3 = 3.0
    addlast!(l, val3)
    @test l.head.value == val2
    @test l.head.next.value == val1
    @test l.head.next.next.value == val3
    @test l.head.next.next.next.prev.value == val3
    @test l.head.next.next.prev.value == val1
    @test l.head.next.prev.value == val2
    @test l.tail.value == val3
    @test l.tail.next.prev.value == val3
    @test length(l) == 3

    return nothing
end

function unit_test_dlinkedlist_removelast()
    l = DLinkedList{Float64}()
    val1 = 1.0
    val2 = 2.0
    addlast!(l, val1)
    addlast!(l, val2)
    @test l.head.value == val1
    @test l.head.next.value == val2
    @test l.head.next.next.prev.value == val2
    @test l.head.next.prev.value == val1
    @test l.tail.value == val2
    @test l.tail.next.prev.value == val2
    @test length(l) == 2
    @test removelast!(l) == val2
    @test l.head.value == val1
    @test l.head.next.prev.value == val1
    @test l.tail.value == val1
    @test l.tail.next.prev.value == val1
    @test length(l) == 1
    @test removelast!(l) == val1
    @test l.head == l.tail
    @test length(l) == 0

    return nothing
end

function unit_test_dlinkedlist_removefirst()
    l = DLinkedList{Float64}()
    val1 = 1.0
    val2 = 2.0
    addlast!(l, val1)
    addlast!(l, val2)
    @test l.head.value == val1
    @test l.head.next.value == val2
    @test l.head.next.next.prev.value == val2
    @test l.head.next.prev.value == val1
    @test l.tail.value == val2
    @test l.tail.next.prev.value == val2
    @test length(l) == 2
    @test removefirst!(l) == val1
    @test l.head.value == val2
    @test l.head.next.prev.value == val2
    @test l.tail.value == val2
    @test l.tail.next.prev.value == val2
    @test length(l) == 1
    @test removefirst!(l) == val2
    @test l.head == l.tail
    @test length(l) == 0

    return nothing
end

function unit_test_dlinkedlist_remove()
    l = DLinkedList{Float64}()
    val1 = 1.0
    val2 = 2.0
    val3 = 3.0
    addlast!(l, val1)
    addlast!(l, val2)
    addlast!(l, val3)
    @test l.head.value == val1
    @test l.head.next.value == val2
    @test l.head.next.next.value == val3
    @test l.head.next.next.next.prev.value == val3
    @test l.head.next.next.prev.value == val2
    @test l.head.next.prev.value == val1
    @test l.tail.value == val3
    @test l.tail.next.prev.value == val3
    @test length(l) == 3
    @test remove!(l, 2) == val2
    @test l.head.value == val1
    @test l.head.next.value == val3
    @test l.head.next.next.prev.value == val3
    @test l.head.next.prev.value == val1
    @test l.tail.value == val3
    @test l.tail.next.prev.value == val3
    @test length(l) == 2
    @test remove!(l, 0) == val1
    @test l.head.value == val3
    @test l.head.next.prev.value == val3
    @test l.tail.value == val3
    @test l.tail.next.prev.value == val3
    @test length(l) == 1
    @test remove!(l, 2) == val3
    @test l.head == l.tail
    @test length(l) == 0

    return nothing
end

function unit_test_dlinkedlist_remove_interface()
    l = DLinkedList{Float64}()
    val1 = 1.0
    val2 = 2.0
    addlast!(l, val1)
    addlast!(l, val2)
    @test l.head.value == val1
    @test l.head.next.value == val2
    @test l.head.next.next.prev.value == val2
    @test l.head.next.prev.value == val1
    @test l.tail.value == val2
    @test l.tail.next.prev.value == val2
    @test length(l) == 2
    @test remove!(l) == val1
    @test l.head.value == val2
    @test l.head.next.prev.value == val2
    @test l.tail.value == val2
    @test l.tail.next.prev.value == val2
    @test length(l) == 1
    @test remove!(l) == val2
    @test l.head == l.tail
    @test length(l) == 0

    return nothing
end

function unit_test_dlinkedlist_removelastfirst()
    l = DLinkedList{Float64}()
    val1 = 1.0
    val2 = 2.0
    val3 = 3.0
    addlast!(l, val1)
    addlast!(l, val2)
    addlast!(l, val3)
    @test l.head.value == val1
    @test l.head.next.value == val2
    @test l.head.next.next.value == val3
    @test l.head.next.next.next.prev.value == val3
    @test l.head.next.next.prev.value == val2
    @test l.head.next.prev.value == val1
    @test length(l) == 3
    @test removelast!(l) == val3
    @test l.head.value == val1
    @test l.head.next.value == val2
    @test l.head.next.next.prev.value == val2
    @test l.head.next.prev.value == val1
    @test l.tail.value == val2
    @test l.tail.next.prev.value == val2
    @test length(l) == 2
    @test removefirst!(l) == val1
    @test l.head.value == val2
    @test l.head.next.prev.value == val2
    @test l.tail.value == val2
    @test l.tail.next.prev.value == val2
    @test length(l) == 1
    @test removelast!(l) == val2
    @test l.head == l.tail
    @test length(l) == 0

    return nothing
end

function unit_test_dlinkedlist_peekfirst()
    val1 = 1.0
    val2 = 2.0
    l = DLinkedList{Float64}()
    addlast!(l, val1)
    @test peekfirst(l) == val1
    addlast!(l, val2)
    @test peekfirst(l) == val1

    return nothing
end

function unit_test_dlinkedlist_peeklast()
    val1 = 1.0
    val2 = 2.0
    l = DLinkedList{Float64}()
    addlast!(l, val1)
    @test peeklast(l) == val1
    addlast!(l, val2)
    @test peeklast(l) == val2

    return nothing
end

function unit_test_linkedlist_length()
    val1 = 1.0
    val2 = 2.0
    l = SLinkedList{Float64}()
    addlast!(l, val1)
    @test length(l) == 1
    addlast!(l, val2)
    @test length(l) == 2

    d = DLinkedList{Float64}()
    addlast!(d, val1)
    @test length(d) == 1
    addlast!(d, val2)
    @test length(d) == 2
end

function unit_test_linkedlist_size()
    val1 = 1.0
    val2 = 2.0
    l = SLinkedList{Float64}()
    addlast!(l, val1)
    @test size(l) == 1
    addlast!(l, val2)
    @test size(l) == 2

    d = DLinkedList{Float64}()
    addlast!(d, val1)
    @test size(d) == 1
    addlast!(d, val2)
    @test size(d) == 2
end

function unit_test_linkedlist_equality()
    val1 = 1.0
    val2 = 2.0
    l1 = SLinkedList{Float64}()
    addlast!(l1, val1)
    addlast!(l1, val2)
    l2 = SLinkedList{Float64}()
    addlast!(l2, val1)
    addlast!(l2, val2)
    @test l1 == l2

    d1 = DLinkedList{Float64}()
    addlast!(d1, val1)
    addlast!(d1, val2)
    d2 = DLinkedList{Float64}()
    addlast!(d2, val1)
    addlast!(d2, val2)
    @test d1 == d2
end

function unit_test_linkedlist_eltype()

end

function unit_test_linkedlist_iterate()
    n = 100
    c = rand((1:n), n)
    l = SLinkedList{Int64}()
    for val in c
        addlast!(l, val)
    end
    for (i, v) in enumerate(l)
        @test v == c[i]
    end

    d = DLinkedList{Int64}()
    for val in c
        addlast!(d, val)
    end
    for (i, v) in enumerate(d)
        @test v == c[i]
    end

    return nothing
end

function unit_test_linkedlist_collect()
    n = 100
    c = rand((1:n), n)
    l = SLinkedList{Int64}()
    for val in c
        addlast!(l, val)
    end
    @test c == collect(l)

    d = DLinkedList{Int64}()
    for val in c
        addlast!(d, val)
    end
    @test c == collect(d)

    return nothing
end

# Run the unit tests
@testset "Singly linked list   " begin
    unit_test_slinkedlist_empty()
    unit_test_slinkedlist_addlast()
    unit_test_slinkedlist_addfirst()
    unit_test_slinkedlist_add()
    unit_test_slinkedlist_add_interface()
    unit_test_slinkedlist_addlastfirst()
    unit_test_slinkedlist_removelast()
    unit_test_slinkedlist_removefirst()
    unit_test_slinkedlist_remove()
    unit_test_slinkedlist_remove_interface()
    unit_test_slinkedlist_removelastfirst()
    unit_test_slinkedlist_peekfirst()
    unit_test_slinkedlist_peeklast()
end

@testset "Doubly linked list   " begin
    unit_test_dlinkedlist_empty()
    unit_test_dlinkedlist_addlast()
    unit_test_dlinkedlist_addfirst()
    unit_test_dlinkedlist_add()
    unit_test_dlinkedlist_add_interface()
    unit_test_dlinkedlist_addlastfirst()
    unit_test_dlinkedlist_removelast()
    unit_test_dlinkedlist_removefirst()
    unit_test_dlinkedlist_remove()
    unit_test_dlinkedlist_remove_interface()
    unit_test_dlinkedlist_removelastfirst()
    unit_test_dlinkedlist_peekfirst()
    unit_test_dlinkedlist_peeklast()
end

@testset "Linked list interface" begin
    unit_test_linkedlist_length()
    unit_test_linkedlist_size()
    unit_test_linkedlist_equality()
    unit_test_linkedlist_eltype()
    unit_test_linkedlist_iterate()
    unit_test_linkedlist_collect()
end


text:
"""
    DynamicAggregation

Functor object for using dynamic aggregation in an L-shaped algorithm. Create by supplying a [`DynamicAggregate`](@ref) object through `aggregate ` in the `LShapedSolver` factory function and then pass to a `StochasticPrograms.jl` model.

The following selection rules are available
- [`SelectUniform`](@ref)
- [`SelectDecaying`](@ref)
- [`SelectRandom`](@ref
- [`SelectClosest`](@ref)
- [`SelectClosestToReference`](@ref)

...
# Parameters
- `num_aggregates::Int`: Number of aggregates
- `rule::SelectionRule`: Rule that determines which aggregate an incoming cut should be placed in
- `lock_after::Function = (τ,n)->false`: Function that determines if the current aggregation scheme should be fixed, based on the current optimality gap `τ` and the number of iterations `n`
...
"""
struct DynamicAggregation{T <: AbstractFloat, S <: AbstractSelectionRule} <: AbstractAggregation
    aggregates::Vector{AggregatedOptimalityCut{T}}
    rule::S
    partitioning::Dict{Int,Int}
    lock::Function

    function DynamicAggregation(num_aggregates::Integer, rule::AbstractSelectionRule, lock::Function, ::Type{T}) where T <: AbstractFloat
        S = typeof(rule)
        aggregates = [zero(AggregatedOptimalityCut{T}) for _ = 1:num_aggregates]
        return new{T,S}(aggregates, rule, Dict{Int,Int}(), lock)
    end
end

function aggregate_cut!(lshaped::AbstractLShaped, ::DynamicAggregation, cut::HyperPlane)
    return add_cut!(lshaped, cut)
end

function aggregate_cut!(lshaped::AbstractLShaped, aggregation::DynamicAggregation{T}, cut::HyperPlane{OptimalityCut}) where T <: AbstractFloat
    if aggregation.lock(gap(lshaped), num_iterations(lshaped)) && haskey(aggregation.partitioning, cut.id)
        aggregation.aggregates[aggregation.partitioning[cut.id]] += cut
        return false
    end
    (idx, full) = select(aggregation.rule, aggregation.aggregates, cut)
    aggregation.aggregates[idx] += cut
    aggregation.partitioning[cut.id] = idx
    if full
        aggregate = aggregation.aggregates[idx]
        aggregation.aggregates[idx] = zero(AggregatedOptimalityCut{T})
        return add_cut!(lshaped, aggregate)
    end
    return false
end

function aggregate_cut!(cutqueue::CutQueue, aggregation::DynamicAggregation, ::MetaData, t::Integer, cut::HyperPlane, x::AbstractArray)
    put!(cutqueue, (t, cut))
    return nothing
end

function aggregate_cut!(cutqueue::CutQueue, aggregation::DynamicAggregation{T}, metadata::MetaData, t::Integer, cut::HyperPlane{OptimalityCut}, x::AbstractArray) where T <: AbstractFloat
    gap = fetch(metadata, t, :gap)
    if aggregation.lock(gap, t) && haskey(aggregation.partitioning, cut.id)
        aggregation.aggregates[aggregation.partitioning[cut.id]] += cut
        return nothing
    end
    (idx, full) = select(aggregation.rule, aggregation.aggregates, cut)
    aggregation.partitioning[cut.id] = idx
    aggregation.aggregates[idx] += cut
    if full
        aggregate = aggregation.aggregates[idx]
        put!(cutqueue, (t, aggregate))
        aggregation.aggregates[idx] = zero(AggregatedOptimalityCut{T})
    end
    return nothing
end

function num_thetas(nscenarios::Integer, ::DynamicAggregation)
    return nscenarios
end

function num_thetas(nscenarios::Integer, ::DynamicAggregation, ::AbstractScenarioProblems)
    return nscenarios
end

function flush!(lshaped::AbstractLShaped, aggregation::DynamicAggregation{T}) where T <: AbstractFloat
    added = false
    for (i,aggregate) in enumerate(aggregation.aggregates)
        if !iszero(aggregate)
            added |= add_cut!(lshaped, aggregate)
            aggregation.aggregates[i] = zero(AggregatedOptimalityCut{T})
        end
    end
    reset!(aggregation.rule)
    return added
end

function flush!(cutqueue::CutQueue, aggregation::DynamicAggregation{T}, ::MetaData, t::Integer, x::AbstractArray) where T <: AbstractFloat
    for (i,aggregate) in enumerate(aggregation.aggregates)
        if !iszero(aggregate)
            put!(cutqueue, (t, aggregate))
            aggregation.aggregates[i] = zero(AggregatedOptimalityCut{T})
        end
    end
    reset!(aggregation.rule)
    return nothing
end

# API
# ------------------------------------------------------------
"""
    DynamicAggregate(num_aggregates::Integer, rule::AbstractSelectionRule; lock_after::Function = (τ,n)->false)

Factory object for [`DynamicAggregation`](@ref). Pass to `aggregate ` in the `LShapedSolver` factory function. See ?DynamicAggregation for parameter descriptions.

"""
mutable struct DynamicAggregate <: AbstractAggregator
    num_aggregates::Int
    rule::AbstractSelectionRule
    lock::Function

    function DynamicAggregate(num_aggregates::Integer, rule::AbstractSelectionRule; lock_after = (τ,n)->false)
        return new(num_aggregates, rule, lock_after)
    end
end

struct SelectionRule <: AggregationParameter end

function MOI.get(aggregator::DynamicAggregate, ::SelectionRule)
    return aggregator.rule
end

function MOI.set(aggregator::DynamicAggregate, ::SelectionRule, rule::AbstractSelectionRule)
    return aggregator.rule = rule
end

function (aggregator::DynamicAggregate)(nscenarios::Integer, T::Type{<:AbstractFloat})
    if aggregator.rule isa SelectRandom && aggregator.rule.max == 1
        return NoAggregation()
    end
    return DynamicAggregation(aggregator.num_aggregates, aggregator.rule, aggregator.lock, T)
end

function remote_aggregator(aggregation::DynamicAggregation, ::AbstractScenarioProblems, ::Integer)
    return DynamicAggregate(length(aggregation.aggregates), aggregation.rule; lock_after = aggregation.lock)
end

function str(aggregator::DynamicAggregate)
    return "dynamic aggregation ruled by $(str(aggregator.rule))"
end


text:
# https://adventofcode.com/2020/day/22

example = "
Player 1:
9
2
6
3
1

Player 2:
5
8
4
7
10
"

function read_input(input)
    deck1, deck2 = split(input, "\n\n", keepempty=false)
    deck1 = parse.(Int, split(deck1, '\n', keepempty=false)[2:end])
    deck2 = parse.(Int, split(deck2, '\n', keepempty=false)[2:end])
    return deck1, deck2
end

@assert read_input(example) == ([9, 2, 6, 3, 1], [5, 8, 4, 7, 10])

function combat(deck1, deck2)
    player = 1
    while !isempty(deck1) && !isempty(deck2)
        x = popfirst!(deck1)
        y = popfirst!(deck2)
        if x > y
            append!(deck1, [x, y])
        else
            append!(deck2, [y, x])
        end
    end
    return isempty(deck1) ? deck2 : deck1
end

score(deck) = sum([i * x for (i, x) in enumerate(reverse(deck))])

function part1(input)
    deck1, deck2 = read_input(input)
    return score(combat(deck1, deck2))
end

@assert part1(example) == 306

function recursive_combat(deck1, deck2; verbose=false)
    previous_rounds = Set()
    deck1 = copy(deck1)
    deck2 = copy(deck2)

    verbose && println("Game\n")

    while !isempty(deck1) && !isempty(deck2)
        verbose && println("Round")
        verbose && println("Player's 1 deck: $(deck1)")
        verbose && println("Player's 2 deck: $(deck2)")

        game_state = hash((deck1, deck2))
        if game_state in previous_rounds
            verbose && println("We saw those cards!")
            return deck1, []
        else
            push!(previous_rounds, game_state)
        end

        x = popfirst!(deck1)
        y = popfirst!(deck2)

        verbose && println("Player 1 plays: $x")
        verbose && println("Player 2 plays: $y")

        if (length(deck1) >= x) && (length(deck2) >= y)
            # sub-game
            verbose && println("Playing a sub-game to determine the winner...")
            winner = findfirst(.!isempty.(recursive_combat(deck1[1:x], deck2[1:y])))
            verbose && println("...anyway, back\n")
        else
            winner = x > y ? 1 : 2
            verbose && println("Player $winner wins\n")
        end

        if winner == 1
            append!(deck1, [x, y])
        else
            append!(deck2, [y, x])
        end
    end
    return deck1, deck2
end

function part2(input)
    deck1, deck2 = read_input(input)
    deck1, deck2 = recursive_combat(deck1, deck2)
    winner = isempty(deck1) ? deck2 : deck1
    return score(winner)
end

@assert part2(example) == 291

test = read("data/day-22.txt", String)
println("Part 1: $(result1 = part1(test))")
println("Part 2: $(result2 = part2(test))")

@assert result1 == 32472
@assert result2 == 36463

max_stars_repo_name: twolodzko/advent-of-code

max_issues_repo_name: twolodzko/advent-of-code

max_forks_repo_name: twolodzko/advent-of-code