text:
using JuliaInterpreter
using JuliaInterpreter: enter_call_expr
using Test, InteractiveUtils, CodeTracking

pc = JuliaInterpreter.JuliaProgramCounter(2)
@test convert(Int, pc) == 2
@test convert(Int, pc+1) == 3
@test convert(Int, pc-1) == 1

module Isolated end

function summer(A)
    s = zero(eltype(A))
    for a in A
        s += a
    end
    return s
end

A = [0.12, -.99]
frame = JuliaInterpreter.enter_call(summer, A)
frame2 = JuliaInterpreter.enter_call(summer, A)
@test summer(A) == something(runframe(frame)) == something(runstack(frame2))

A = rand(1000)
@test @interpret(sum(A)) ≈ sum(A)  # note: the compiler can leave things in registers to increase accuracy, doesn't happen with interpreted
fapply() = (Core.apply_type)(Base.NamedTuple, (), Tuple{})
@test @interpret(fapply()) == fapply()
function fbc()
    bc = Broadcast.broadcasted(CartesianIndex, 6, [1, 2, 3])
    copy(bc)
end
@test @interpret(fbc()) == fbc()
@test @interpret(repr("hi")) == repr("hi")  # this tests kwargs and @generated functions

fkw(x::Int8; y=0, z="hello") = y
@test @interpret(fkw(Int8(1); y=22, z="world")) == fkw(Int8(1); y=22, z="world")

# generators that throw before returning the body expression
@test_throws ArgumentError("input tuple of length 3, requested 2") @interpret Base.fill_to_length((1,2,3), -1, Val(2))

# Throwing exceptions across frames
function f_exc_inner()
    error("inner")
end

f_exc_inner2() = f_exc_inner()

const caught = Ref(false)
function f_exc_outer1()
    try
        f_exc_inner()
    catch err    # with an explicit err capture
        caught[] = true
        rethrow(err)
    end
end

function f_exc_outer2()
    try
        f_exc_inner()
    catch        # implicit err capture
        caught[] = true
        rethrow()
    end
end

function f_exc_outer3(f)
    try
        f()
    catch err
        return err
    end
end

@test !caught[]
ret = @interpret f_exc_outer3(f_exc_outer1)
@test ret == ErrorException("inner")
@test caught[]

caught[] = false
ret = @interpret f_exc_outer3(f_exc_outer2)
@test ret == ErrorException("inner")
@test caught[]

caught[] = false
ret = @interpret f_exc_outer3(f_exc_inner2)
@test ret == ErrorException("inner")
@test !caught[]

stc = try f_exc_outer1() catch
    stacktrace(catch_backtrace())
end
sti = try @interpret(f_exc_outer1()) catch
    stacktrace(catch_backtrace())
end
@test_broken stc == sti

# issue #3
@test @interpret(joinpath("/home/julia/base", "sysimg.jl")) == "/home/julia/base/sysimg.jl"
@test @interpret(10.0^4) == 10.0^4
# issue #6
@test @interpret(Array.body.body.name) === Array.body.body.name
@test @interpret(Vararg.body.body.name) === Vararg.body.body.name
frame = JuliaInterpreter.prepare_thunk(Main, :(Vararg.body.body.name))
@test JuliaInterpreter.finish_and_return!(JuliaStackFrame[], frame, true) === Vararg.body.body.name
frame = JuliaInterpreter.prepare_thunk(Base, :(Union{AbstractChar,Tuple{Vararg{<:AbstractChar}},AbstractVector{<:AbstractChar},Set{<:AbstractChar}}))
@test JuliaInterpreter.finish_and_return!(JuliaStackFrame[], frame, true) isa Union

# issue #8
ex = quote
    if sizeof(JLOptions) === ccall(:jl_sizeof_jl_options, Int, ())
    else
        ccall(:jl_throw, Cvoid, (Any,), "Option structure mismatch")
    end
end
frame = JuliaInterpreter.prepare_thunk(Base, ex)
JuliaInterpreter.finish_and_return!(JuliaStackFrame[], frame, true)

# ccall with two Symbols
ex = quote
    @testset "Some tests" begin
       @test 2 > 1
    end
end
frame = JuliaInterpreter.prepare_thunk(Main, ex)
JuliaInterpreter.finish_and_return!(JuliaStackFrame[], frame, true)

@test @interpret Base.Math.DoubleFloat64(-0.5707963267948967, 4.9789962508669555e-17).hi ≈ -0.5707963267948967

# ccall with cfunction
fcfun(x::Int, y::Int) = 1
ex = quote   # in lowered code, cf is a Symbol
    cf = @eval @cfunction(fcfun, Int, (Int, Int))
    ccall(cf, Int, (Int, Int), 1, 2)
end
frame = JuliaInterpreter.prepare_thunk(Main, ex)
@test JuliaInterpreter.finish_and_return!(JuliaStackFrame[], frame, true) == 1
ex = quote
    let   # in lowered code, cf is a SlotNumber
        cf = @eval @cfunction(fcfun, Int, (Int, Int))
        ccall(cf, Int, (Int, Int), 1, 2)
    end
end
frame = JuliaInterpreter.prepare_thunk(Main, ex)
@test JuliaInterpreter.finish_and_return!(JuliaStackFrame[], frame, true) == 1

# From Julia's test/ambiguous.jl. This tests whether we renumber :enter statements correctly.
ambig(x, y) = 1
ambig(x::Integer, y) = 2
ambig(x, y::Integer) = 3
ambig(x::Int, y::Int) = 4
ambig(x::Number, y) = 5
ex = quote
    let
        cf = @eval @cfunction(ambig, Int, (UInt8, Int))
        @test_throws(MethodError, ccall(cf, Int, (UInt8, Int), 1, 2))
    end
end
frame = JuliaInterpreter.prepare_thunk(Main, ex)
JuliaInterpreter.finish_and_return!(JuliaStackFrame[], frame, true)

# Core.Compiler
ex = quote
    length(code_typed(fcfun, (Int, Int)))
end
frame = JuliaInterpreter.prepare_thunk(Main, ex)
@test JuliaInterpreter.finish_and_return!(JuliaStackFrame[], frame, true) == 1

# copyast
ex = quote
    struct CodegenParams
        cached::Cint

        track_allocations::Cint
        code_coverage::Cint
        static_alloc::Cint
        prefer_specsig::Cint

        module_setup::Any
        module_activation::Any
        raise_exception::Any
        emit_function::Any
        emitted_function::Any

        CodegenParams(;cached::Bool=true,
                       track_allocations::Bool=true, code_coverage::Bool=true,
                       static_alloc::Bool=true, prefer_specsig::Bool=false,
                       module_setup=nothing, module_activation=nothing, raise_exception=nothing,
                       emit_function=nothing, emitted_function=nothing) =
            new(Cint(cached),
                Cint(track_allocations), Cint(code_coverage),
                Cint(static_alloc), Cint(prefer_specsig),
                module_setup, module_activation, raise_exception,
                emit_function, emitted_function)
    end
end
frame = JuliaInterpreter.prepare_thunk(Isolated, ex)
JuliaInterpreter.finish_and_return!(JuliaStackFrame[], frame, true)
@test Isolated.CodegenParams(cached=false).cached === Cint(false)

# cglobal
val = @interpret(BigInt())
@test isa(val, BigInt) && val == 0
@test isa(@interpret(Base.GMP.version()), VersionNumber)

# "correct" line numbers
defline = @__LINE__() + 1
function f(x)
    x = 2x
    # comment
    # comment
    x = 2x
    # comment
    return x*x
end
frame = JuliaInterpreter.enter_call(f, 3)
@test whereis(frame, JuliaInterpreter.JuliaProgramCounter(1))[2] == defline + 1
@test whereis(frame, JuliaInterpreter.JuliaProgramCounter(3))[2] == defline + 4
@test whereis(frame, JuliaInterpreter.JuliaProgramCounter(5))[2] == defline + 6

# issue #28
let a = ['0'], b = ['a']
    @test @interpret(vcat(a, b)) == vcat(a, b)
end

# issue #51
if isdefined(Core.Compiler, :SNCA)
    ci = @code_lowered gcd(10, 20)
    cfg = Core.Compiler.compute_basic_blocks(ci.code)
    @test isa(@interpret(Core.Compiler.SNCA(cfg)), Vector{Int})
end

# llvmcall
function add1234(x::Tuple{Int32,Int32,Int32,Int32})
    Base.llvmcall("""%3 = extractvalue [4 x i32] %0, 0
      %4 = extractvalue [4 x i32] %0, 1
      %5 = extractvalue [4 x i32] %0, 2
      %6 = extractvalue [4 x i32] %0, 3
      %7 = extractvalue [4 x i32] %1, 0
      %8 = extractvalue [4 x i32] %1, 1
      %9 = extractvalue [4 x i32] %1, 2
      %10 = extractvalue [4 x i32] %1, 3
      %11 = add i32 %3, %7
      %12 = add i32 %4, %8
      %13 = add i32 %5, %9
      %14 = add i32 %6, %10
      %15 = insertvalue [4 x i32] undef, i32 %11, 0
      %16 = insertvalue [4 x i32] %15, i32 %12, 1
      %17 = insertvalue [4 x i32] %16, i32 %13, 2
      %18 = insertvalue [4 x i32] %17, i32 %14, 3
      ret [4 x i32] %18""",Tuple{Int32,Int32,Int32,Int32},
      Tuple{Tuple{Int32,Int32,Int32,Int32},Tuple{Int32,Int32,Int32,Int32}},
        (Int32(1),Int32(2),Int32(3),Int32(4)),
        x)
end
@test @interpret(add1234(map(Int32,(2,3,4,5)))) === map(Int32,(3,5,7,9))

# issue #74
let A = [1]
    wkd = WeakKeyDict()
    @interpret setindex!(wkd, 2, A)
    @test wkd[A] == 2
end

# issue #76
let TT = Union{UInt8, Int8}
    a = TT[0x0, 0x1]
    pa = pointer(a)
    @interpret unsafe_store!(pa, 0x1, 2)
end

# issue #92
let x = Core.TypedSlot(1, Any)
    f(x) = objectid(x)
    @test isa(@interpret(f(x)), UInt)
end

# Some expression can appear nontrivial but lower to nothing
@test isa(JuliaInterpreter.prepare_thunk(Main, :(@static if ccall(:jl_get_UNAME, Any, ()) == :NoOS 1+1 end)), Nothing)
@test isa(JuliaInterpreter.prepare_thunk(Main, :(Base.BaseDocs.@kw_str "using")), Nothing)

@testset "locals" begin
    f_locals(x::Int64, y::T, z::Vararg{Symbol}) where {T} = x
    frame = JuliaInterpreter.enter_call(f_locals, 1, 2.0, :a, :b)
    locals = JuliaInterpreter.locals(frame)
    @test JuliaInterpreter.Variable(1, :x, false) in locals
    @test JuliaInterpreter.Variable(2.0, :y, false) in locals
    @test JuliaInterpreter.Variable((:a, :b), :z, false) in locals
    @test JuliaInterpreter.Variable(Float64, :T, true) in locals

    function f_multi(x)
        c = x
        x = 2
        x = 3
        x = 4
        return x
    end
    frame = JuliaInterpreter.enter_call(f_multi, 1)
    stack = [frame]
    locals = JuliaInterpreter.locals(frame)
    @test length(locals) == 2
    @test JuliaInterpreter.Variable(1, :x, false) in locals
    JuliaInterpreter.step_expr!(stack, frame)
    JuliaInterpreter.step_expr!(stack, frame)
    locals = JuliaInterpreter.locals(frame)
    @test length(locals) == 3
    @test JuliaInterpreter.Variable(1, :c, false) in locals
    JuliaInterpreter.step_expr!(stack, frame)
    locals = JuliaInterpreter.locals(frame)
    @test length(locals) == 3
    @test JuliaInterpreter.Variable(2, :x, false) in locals
    JuliaInterpreter.step_expr!(stack, frame)
    locals = JuliaInterpreter.locals(frame)
    @test length(locals) == 3
    @test JuliaInterpreter.Variable(3, :x, false) in locals
end


text:
# Test routine to compare scaling for a number of QP Lasso problems (partially badly scaled)
workspace()
include("../../../src/Solver.jl")
include("../../Benchmarks/solverComparison/Compare.jl")

using OSQP, OSSDP, Base.Test, Compare, JuMP, Mosek

rng = Random.MersenneTwister(12345)
#rng = Random.MersenneTwister(211121)

nn = 20
timestamp = Dates.format(now(), "yyddmm_HH-MM")
dataFile = "./SC_" * timestamp * "LassoQP.jld"
problemType = "QP-Lasso"

sr1 = SolverResult(nn, problemType,"COSMO",timestamp,0,false)
sr2 = SolverResult(nn, problemType,"COSMO",timestamp,0,true)
sr3 = SolverResult(nn, problemType,"OSQP",timestamp,0,false)
sr4 = SolverResult(nn, problemType,"OSQP",timestamp,0,true)
resData = [sr1;sr2;sr3;sr4]
ws = 0
for iii =1:1:nn
  # generate problem data
  n = 8
  m = 50*n
  A = sprandn(rng,m,n,0.5)
  vtrue = 1/n*sprandn(rng,n,0.5)
  noise = 1/4*randn(rng,m,1)
  b = A*vtrue + noise
  λ = 0.2*norm(A'*b,Inf)

  pDims = [m;n;nnz(A)]

  # define lasso problem as QP
  Aa = [-A zeros(m,n) Matrix(1.0I,m,m);
         Matrix(1.0I,n,n) -Matrix(1.0I,n,n) zeros(n,m);
         -Matrix(1.0I,n,n) -Matrix(1.0I,n,n) zeros(n,m)]

  ba = [-b;zeros(2*n)]
  P = 2*Matrix(Diagonal([zeros(2*n);ones(m)])) # times two to cancel the 1/2 in the cost function
  q = [zeros(n);λ*ones(n);zeros(m)]

  # define cone membership
  K = OSSDPTypes.Cone(m,2*n,[],[])

  setOFF = OSSDPSettings(rho=0.1,sigma=1e-6,alpha=1.6,max_iter=1500,verbose=false,check_termination=1,scaling = 0,eps_abs = 1e-3,eps_rel=1e-3,adaptive_rho=true)
  setON = OSSDPSettings(rho=0.1,sigma=1e-6,alpha=1.6,max_iter=1500,verbose=false,check_termination=1,scaling = 10,eps_abs = 1e-3,eps_rel=1e-3,adaptive_rho=true)

  # Solve with OSSDP
  resOSSDP_unscaled,nothing = OSSDP.solve(P,q,Aa,ba,K,setOFF);
        print("\n.")

  resOSSDP_scaled,nothing = OSSDP.solve(P,q,Aa,ba,K,setON);
      print(".")

  # modify problem for OSQP (using inequality constraints)
  l = full([-b;-Inf*ones(n);zeros(n)][:])
  u = full([-b;zeros(n);Inf*ones(n)][:])
  Aa2 = [-A zeros(m,n) Matrix(1.0I,m,m);
       Matrix(1.0I,n,n) -Matrix(1.0I,n,n) zeros(n,m);
       Matrix(1.0I,n,n) Matrix(1.0I,n,n) zeros(n,m)]

  m1 = OSQP.Model()
  OSQP.setup!(m1; P=sparse(P), q=q, A=sparse(Aa2), l=l, u=u,scaling=0,check_termination=1,verbose=false,adaptive_rho = true,eps_abs = 1e-3,eps_rel=1e-3)
  resOSQP1_unscaled = OSQP.solve!(m1)
  m2 = OSQP.Model()
  OSQP.setup!(m2; P=sparse(P), q=q, A=sparse(Aa2), l=l, u=u,scaling=10,check_termination=1,verbose=false,adaptive_rho = true,eps_abs = 1e-3,eps_rel=1e-3)
  resOSQP1_scaled = OSQP.solve!(m2)

  resArray = [resOSSDP_unscaled;resOSSDP_scaled;resOSQP1_unscaled;resOSQP1_scaled]
  updateResults!(dataFile,resData,resArray,pDims,"QP-$(iii)",0.,true)
  printStatus(iii,nn,"QP-$(iii)",resData)
end

ext: jl


text:
@testset "KBS" begin

    A = Alphabet(['a', 'e', 'b', 'p'])
    KnuthBendix.set_inversion!(A, 'a', 'e')
    KnuthBendix.set_inversion!(A, 'b', 'p')

    a = Word([1,2])
    b = Word([2,1])
    c = Word([3,4])
    d = Word([4,3])
    ε = one(a)

    ba = Word([3,1])
    ab = Word([1,3])

    be = Word([3,2])
    eb = Word([2,3])

    pa = Word([4,1])
    ap = Word([1,4])

    pe = Word([4,2])
    ep = Word([2,4])

    lenlexord = LenLex(A)
    rs = RewritingSystem([a=>ε, b=>ε, c=>ε, d=>ε, ba=>ab], lenlexord)

    crs = Set(KnuthBendix.Rule.([a=>ε, b=>ε, c=>ε, d=>ε, ba=>ab, be=>eb, pa=>ap, pe=>ep]))

    @test Set(KnuthBendix.rules(KnuthBendix.knuthbendix(rs, implementation=:naive_kbs1))) == crs
    @test Set(KnuthBendix.rules(KnuthBendix.knuthbendix(rs, implementation=:naive_kbs2))) == crs
    @test Set(KnuthBendix.rules(KnuthBendix.knuthbendix(rs, implementation=:deletion))) == crs
    @test_throws String Set(KnuthBendix.rules(KnuthBendix.knuthbendix(rs, implementation=:automata))) == crs
end


text:
function jldopen(@nospecialize(f::Function), args...; kws...)
    jld = jldopen(args...; kws...)
    try
        return f(jld)
    finally
        close(jld)
    end
end

"""
    @save filename var1 [var2 ...]
    @save filename {compress=true} var1 name2=var2

Write one or more variables `var1,...` from the current scope to a JLD2 file
`filename`.

For interactive use you can save all variables in the current module's global
scope using `@save filename`. More permanent code should prefer the explicit
form to avoid saving unwanted variables.

# Example

To save the string `hello` and array `xs` to the JLD2 file example.jld2:

    hello = "world"
    xs = [1,2,3]
    @save "example.jld2" hello xs

For passing options to the saving command use {}

    @save "example.jld2" {compress=true} hello xs

For saving variables under a different name use regular assignment syntax

    @save "example.jld2" greeting=hello xarray = xs
"""
macro save(filename, vars...)
    fields = []
    options = []
    for var in vars
        # Capture options of the form {compress = true, mmaparrays=false}
        if @capture(var, {opts__})
            for opt in opts
                if @capture(opt, key_ = val_)
                    push!(options, :($key = $(esc(val))))
                else
                    return :(throw(ArgumentError("Invalid option syntax")))
                end
            end
        # Allow assignment syntax a = b
        elseif @capture(var, a_ = b_)
            push!(fields, :(write(f, $(string(a)), $(esc(b)), wsession)))
        # Allow single arg syntax a   → "a" = a
        elseif @capture(var, a_Symbol)
            push!(fields, :(write(f, $(string(a)), $(esc(a)), wsession)))
        else
            return :(throw(ArgumentError("Invalid field syntax")))
        end
    end
    if !isempty(fields)
        return quote
            let
                f = jldopen($(esc(filename)), "w"; $(Expr(:tuple,options...))...)
                wsession = JLDWriteSession()
                try
                    $(Expr(:block, fields...))
                catch e
                    rethrow(e)
                finally
                    close(f)
                end
            end
        end
    else
        # The next part is old code that handles saving the whole workspace
        quote
            let
                m = $(__module__)
                f = jldopen($(esc(filename)), "w"; $(Expr(:tuple,options...))...)
                wsession = JLDWriteSession()
                try
                    for vname in names(m; all=true)
                        s = string(vname)
                        if !occursin(r"^_+[0-9]*$", s) # skip IJulia history vars
                            (startswith(s, "##") || s == "ans") && continue
                            v = getfield(m, vname)
                            if !isa(v, Module)
                                try
                                    write(f, s, v, wsession)
                                catch e
                                    if isa(e, PointerException)
                                        @warn("skipping $vname because it contains a pointer")
                                    else
                                        rethrow(e)
                                    end
                                end
                            end
                        end
                    end
                finally
                    close(f)
                end
            end
        end
    end
end

"""
    @load filename var1 [var2 ...]

Load one or more variables `var1,...` from JLD2 file `filename` into the
current scope and return a vector of the loaded variable names.

For interactive use, the form `@load "somefile.jld2"` will load all variables
from `"somefile.jld2"` into the current scope. This form only supports literal
file names and should be avoided in more permanent code so that it's clear
where the variables come from.

# Example

To load the variables `hello` and `foo` from the file example.jld2, use

    @load "example.jld2" hello foo
"""
macro load(filename, vars...)
    if isempty(vars)
        if isa(filename, Expr)
            throw(ArgumentError("filename argument must be a string literal unless variable names are specified"))
        end
        # Load all variables in the top level of the file
        readexprs = Expr[]
        vars = Symbol[]
        f = jldopen(filename)
        try
            for n in keys(f)
                if !isgroup(f, lookup_offset(f.root_group, n))
                    push!(vars, Symbol(n))
                end
            end
        finally
            close(f)
        end
    end
    return quote
        ($([esc(x) for x in vars]...),) = jldopen($(esc(filename))) do f
            ($([:(read(f, $(string(x)))) for x in vars]...),)
        end
        $(Symbol[v for v in vars]) # convert to Array
    end
end

function loadtodict!(d::Dict, g::Union{JLDFile,Group}, prefix::String="")
    for k in keys(g)
        v = g[k]
        if v isa Group
            loadtodict!(d, v, prefix*k*"/")
        else
            d[prefix*k] = v
        end
    end
    return d
end

# Name used in JLD2 file to identify objects stored with `save_object`
const SINGLE_OBJECT_NAME = "single_stored_object"

"""
    save_object(filename, x)

Stores an object `x` in a new JLD2 file at `filename`. If a file exists at this
path, it will be overwritten.

Since the JLD2 format requires that all objects have a name, the object will be
stored as `single_sotred_object`. If you want to store more than one object, use
[`@save`](@ref) macro, [`jldopen`](@ref) or the FileIO API.

# Example

To save the string `hello` to the JLD2 file example.jld2:

    hello = "world"
    save_object("example.jld2", hello)
"""
function save_object(filename, x)
  jldopen(filename, "w") do file
    file[SINGLE_OBJECT_NAME] = x
  end
  return
end

"""
    load_object(filename)

Returns the only available object from the JLD2 file `filename` (The stored
object name is inconsequential). If the file contains more than one or no
objects, the function throws an `ArgumentError`.

For loading more than one object, use [`@load`](@ref) macro, [`jldopen`](@ref)
or the FileIO API.

# Example

To load the only object from the JLD2 file example.jld2:

    hello = "world"
    save_object("example.jld2", hello)
    hello_loaded = load_object("example.jld2")
"""
function load_object(filename)
  jldopen(filename, "r") do file
    all_keys = keys(file)
    length(all_keys) == 0 && throw(ArgumentError("File $filename does not contain any object"))
    length(all_keys) > 1 && throw(ArgumentError("File $filename contains more than one object. Use `load` or `@load` instead"))
    file[all_keys[1]] #Uses HDF5 functionality of treating the file like a dict
  end
end

"""
    jldsave(filename, compress=false; kwargs...)

Creates a JLD2 file at `filename` and stores the variables given as keyword arguments.

# Examples

    jldsave("example.jld2"; a=1, b=2, c)

is equivalent to

    jldopen("example.jld2, "w") do f
        f["a"] = 1
        f["b"] = 2
        f["c"] = c
    end

To choose the io type `IOStream` instead of the default `MmapIO` use
`jldsave(fn, IOStream; kwargs...)`.
"""
function jldsave(filename::AbstractString, compress=false, iotype::T=MmapIO;
                    kwargs...
                    ) where T<:Union{Type{IOStream},Type{MmapIO}}
    jldopen(filename, "w"; compress=compress, iotype=iotype) do f
        wsession = JLDWriteSession()
        for (k,v) in pairs(kwargs)
            write(f, string(k), v, wsession)
        end
    end
end

jldsave(filename::AbstractString, iotype::Union{Type{IOStream},Type{MmapIO}}; kwargs...) =
    jldsave(filename, false, iotype; kwargs...)


text:
"""
    function exponentiate(A, t::Number, x; kwargs...)
    function exponentiate(A, t::Number, x, algorithm)

Compute ``y = exp(t*A) x``, where `A` is a general linear map, i.e. a `AbstractMatrix` or
just a general function or callable object and `x` is of any Julia type with vector like
behavior.

### Arguments:

The linear map `A` can be an `AbstractMatrix` (dense or sparse) or a general function or
callable object that implements the action of the linear map on a vector. If `A` is an
`AbstractMatrix`, `x` is expected to be an `AbstractVector`, otherwise `x` can be of any
type that behaves as a vector and supports the required methods (see KrylovKit docs).

The time parameter `t` can be real or complex, and it is better to choose `t` e.g. imaginary
and `A` hermitian, then to absorb the imaginary unit in an antihermitian `A`. For the
former, the Lanczos scheme is used to built a Krylov subspace, in which an approximation to
the exponential action of the linear map is obtained. The argument `x` can be of any type
and should be in the domain of `A`.

### Return values:

The return value is always of the form `y, info = exponentiate(...)` with

  - `y`: the result of the computation, i.e. `y = exp(t*A)*x`

  - `info`: an object of type [`ConvergenceInfo`], which has the following fields

      + `info.converged::Int`: 0 or 1 if the solution `y` was approximated up to the
        requested tolerance `tol`.
      + `info.residual::Nothing`: value `nothing`, there is no concept of a residual in
        this case
      + `info.normres::Real`: a (rough) estimate of the error between the approximate and
        exact solution
      + `info.numops::Int`: number of times the linear map was applied, i.e. number of times
        `f` was called, or a vector was multiplied with `A`
      + `info.numiter::Int`: number of times the Krylov subspace was restarted (see below)

!!! warning "Check for convergence"

    By default (i.e. if `verbosity = 0`, see below), no warning is printed if the solution
    was not found with the requested precision, so be sure to check `info.converged == 1`.

### Keyword arguments:

Keyword arguments and their default values are given by:

  - `verbosity::Int = 0`: verbosity level, i.e. 0 (no messages), 1 (single message
    at the end), 2 (information after every iteration), 3 (information per Krylov step)
  - `krylovdim = 30`: the maximum dimension of the Krylov subspace that will be constructed.
    Note that the dimension of the vector space is not known or checked, e.g. `x₀` should
    not necessarily support the `Base.length` function. If you know the actual problem
    dimension is smaller than the default value, it is useful to reduce the value of
    `krylovdim`, though in principle this should be detected.
  - `tol = 1e-12`: the requested accuracy per unit time, i.e. if you want a certain
    precision `ϵ` on the final result, set `tol = ϵ/abs(t)`. If you work in e.g. single
    precision (`Float32`), you should definitely change the default value.
  - `maxiter::Int = 100`: the number of times the Krylov subspace can be rebuilt; see below
    for further details on the algorithms.
  - `issymmetric`: if the linear map is symmetric, only meaningful if `T<:Real`
  - `ishermitian`: if the linear map is hermitian
    The default value for the last two depends on the method. If an `AbstractMatrix` is
    used, `issymmetric` and `ishermitian` are checked for that matrix, otherwise the default
    values are `issymmetric = false` and `ishermitian = T <: Real && issymmetric`.
  - `eager::Bool = false`: if true, eagerly try to compute the result after every expansion
    of the Krylov subspace to test for convergence, otherwise wait until the Krylov subspace
    as dimension `krylovdim`. This can result in a faster return, for example if the total
    time for the evolution is quite small, but also has some overhead, as more computations
    are performed after every expansion step.

### Algorithm

This is actually a simple wrapper over more general method [`expintegrator`](@ref) for
for integrating a linear non-homogeneous ODE.
"""
function exponentiate end

exponentiate(A, t::Number, v; kwargs...) = expintegrator(A, t, v; kwargs...)
exponentiate(A, t::Number, v, alg::Union{Lanczos,Arnoldi}) = expintegrator(A, t, (v,), alg)


text:
module TestGrouping

using Base.Test
using DataArrays, DataFrames
using DataFramesMeta

d = DataFrame(n = 1:20, x = [3, 3, 3, 3, 1, 1, 1, 2, 1, 1, 2, 1, 1, 2, 2, 2, 3, 1, 1, 2])
g = groupby(d, :x)

@test  @where(d, :x .== 3) == where(d, x -> x[:x] .== 3)
@test  DataFrame(@where(g, length(:x) > 5)) == DataFrame(where(g, x -> length(x[:x]) > 5))
@test  DataFrame(@where(g, length(:x) > 5))[:n][1:3] == @data [5, 6, 7]

@test  DataFrame(orderby(g, x -> mean(x[:n]))) == DataFrame(@orderby(g, mean(:n)))
#@test  DataFrames.based_on(@orderby(g, mean(:n)), x -> x[1,:x])[:x1] == [3,1,2]

@test  (@transform(g, y = :n - median(:n)))[1,:y] == -5.0

end # module


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

module Repeat

export basismatrix

using Base: offset_if_vec
using StaticArrays:SVector,MVector,SMatrix,MMatrix
using DataFrames:DataFrame
using Colors
import LazySets

include("Lattice.jl")
include("HexagonalLattice.jl")
include("RectangularLattice.jl")
include("Array.jl")
include("Cluster.jl")
include("Multilens/Lenslets.jl")

end #module
export Repeat


text:
# Linear GP circuit chromosome functions
# When Parameters are used, numlevelsback specifies numregisters.
# A circuit is a sequence of instructions on the register array R
#   which has length numregisters + numinputs
# R is normally local to the function where is is used.
# The first numregisters components or R are calculation registers.
# The remaining numinputs components are read-only input registers which are
#   initialized with the components of the current context.
# Even though there is a nodearity keyword value, this file assumes nodearity==2
# An instruction is specified either by a 4-tuple of MyInts or by an integer of OutputType.
# The elements of the (2+nodearity)-tuple are (assuming nodearity==2):
#  1.  The index of the logical operation (such as AND, OR) in the funcs array.
#  2.  The index of the element of R where the output is stored
#  3.  The index of the element of R which is the first operand of the logical operation
#  4.  The index of the element of R which is the second operand of the logical operation
# The functions vect_to_int() and int_to_vect() convert between these representations.
# When parameters are used, numinputs and numoutputs have the same meaning for the cartesian and linear representatons.
# numinteriors means numinstructions for the linear representation
# numlevelsback means numregisters for the linear representation
# There are 3 alternative representations for the circuit_vects component of a  LinCicuit.
#   1.  A list of numinstructions vectors where each vector is  (2+nodearity)-tuple described above.
#   2.  A list of numinstructions Int64s where each of these is a 1-based encoding of the vectors
#   3.  A Int128 integer which is an encoding of the integers of step 3.
# See notes/3_24.txt for outline
# See test/testLinCircuit.jl for tests.
export LinCircuit, output_values
export execute_lcircuit, instruction_vect_to_instruction_int, instruction_int_to_instruction_vect, lincomplexity
export vect_to_int, int_to_vect, rand_ivect, execute_random_circuit, number_active, degeneracy, recover_phenotype
export num_instructions # This is the total number of possible instructions for a parameter setting
export rand_lcircuit, mutate_circuit!, mutate_instruction, mutate_circuit_all, mutate_all, print_circuit
export instruction_vects_to_instruction_ints, instruction_ints_to_instruction_vects
export instruction_ints_to_circuit_int, circuit_int_to_instruction_ints
export circuit_int_to_circuit, circuit_to_circuit_int
#export circuit_int
OutputType = Int64

#=
mutable struct LinCircuit
    # Here, numinteriors represents number of gates, numlevelsback represents number of computational registers
    # the First numoutputs registers are the output registers
    circuit_vects::Vector{Vector{MyInt}}
    params::Parameters
end
=#

function output_values( c::LinCircuit, funcs::Vector{Func}=Func[] )
  if length(funcs) == 0
    funcs = lin_funcs( c.params.numinputs )
  end
  #println("funcs: ",funcs)
  R = execute_lcircuit( c, funcs )
  R[1:c.params.numoutputs]
end

# Not used.  Assumes that R and funcs are in the execution environment
function vect_to_funct( lc::Vector{MyInt} )
  (lc)->R[lc[2]] = funcs[lc[1]].func(R[lc[3]],R[lc[4]])
end

# Executes the circuit, returns the register array R.
# numregisters is the number of calculation registers, not the total number of registers
# The outputs are  R[1:numoutputs]
function execute_lcircuit( circuit_ints::Vector{OutputType}, numregisters::Int64, numinputs::Int64, funcs::Vector{Func}; nodearity::Int64=2 )
  lcv = map( ci->instruction_int_to_instruction_vect( ci, numregisters, numinputs, funcs ), circuit_ints )
  execute_lcircuit( lcv, numregisters, numinputs, funcs )
end

function execute_lcircuit( circuit::LinCircuit, funcs::Vector{Func} )
  p = circuit.params
  R = fill(MyInt(0), p.numlevelsback + p.numinputs ) # numlevelsback is the number of computational registers
  R[(p.numlevelsback+1):end] = construct_context(p.numinputs)
  for lc in circuit.circuit_vects
    #println("lc: ",lc,"  R: ",R)
    R[lc[2]] = funcs[lc[1]].func(R[lc[3]],R[lc[4]])
  end
  R
end

# Executes the circuit, returns the register array R.
# The outputs are  R[1:numoutputs]
function execute_lcircuit( circuit_vects::Vector{Vector{MyInt}}, numregisters::Int64, numinputs::Int64, funcs::Vector{Func}; nodearity::Int64=2 )
  R = fill(MyInt(0), numregisters+numinputs )
  R[numregisters+1:end] = construct_context(numinputs)
  for lc in circuit_vects
    #println("lc: ",lc,"  func: ",funcs[lc[1]],"  R: ",R)
    R[lc[2]] = funcs[lc[1]].func(R[lc[3]],R[lc[4]])
    #println("lc: ",lc,"  func: ",funcs[lc[1]],"  R: ",R)
  end
  R
end

function lincomplexity( circuit_vects::Vector{Vector{MyInt}}, numregisters::Int64, numinputs::Int64, funcs::Vector{Func}; nodearity::Int64=2 )
  p = Parameters( numinputs, 1, length(circuit_vects), numregisters )
  lincomplexity( LinCircuit(circuit_vects,p), funcs, nodearity=nodearity )
end

# Executes the circuit, returns the array X of instruction outputs
# The outputs are  R[1:numoutputs]
function lincomplexity( circuit::LinCircuit, funcs::Vector{Func}; nodearity::Int64=2 )
  p = circuit.params
  R = fill(MyInt(0), p.numlevelsback+p.numinputs )
  X = fill(MyInt(0), length(circuit.circuit_vects) )
  R[p.numlevelsback+1:end] = construct_context(p.numinputs)
  i = 1
  for lc in circuit.circuit_vects
    #println("lc: ",lc,"  func: ",funcs[lc[1]],"  R: ",R)
    R[lc[2]] = X[i] = funcs[lc[1]].func(R[lc[3]],R[lc[4]])
    #println("lc: ",lc,"  func: ",funcs[lc[1]],"  R: ",R)
    i += 1
  end
  complexity5( X, p.numinputs )
end

# The number of possible instructions with these settings.
# This is distinct from numinstructions which is the number of instructions in a circuit.
function num_instructions( numregisters::Int64, numinputs::Int64, funcs::Vector{Func};
     nodearity::Int64=2 )
  length(funcs)*numregisters*(numinputs+numregisters)^nodearity
end

function num_instructions( p::Parameters, funcs::Vector{Func} )
  num_instructions( p.numlevelsback, p.numinputs, funcs, nodearity=p.nodearity )
end

function number_active( circuit::LinCircuit )  # Dummy function for now 10/16/21
  return 0
end

function degeneracy( circuit::LinCircuit )  # Dummy function for now 10/16/21
  return 0
end

function recover_phenotype( circuit::LinCircuit, max_steps_recover, maxtrials_recover, maxtries_recover )  # Dummy function for now 10/16/21
  return (0,0)
end

# Convert a vectors of integers to a unique integer.
# The components of the vector are 1-based rather than 0-based.
# Example 1:
# julia> vect_to_int( [3,4,5,6],[10,10,10,10])
#  2345
# Example 2:
# julia> vect_to_int([3,2,1,4],[7,2,5,8])
#  203
function vect_to_int( ints::Vector{Int64}, max::Vector{Int64} )
  result = Int128(ints[1])-1
  for i = 2:length(ints)
    result *= max[i]
    result += ints[i]-1
  end
  result
end

# julia> int_to_vect( 2345, [10,10,10,10])
# 4-element Array{Int64,1}:
#  Int64[ 3, 4, 5, 6 ]
# julia> int_to_vect(203,[7,2,5,8])
# 4-element Array{Int64,1}:
#  Int64[ 3, 2, 1, 4]
function int_to_vect( int::Integer, max::Vector{Int64} )
  ints = zeros(Int64,length(max))
  for i = length(max):-1:1
    ints[i] = int % max[i] + 1
    int = div(int,max[i])
  end
  ints
end

# Converts an instruction specified by a vector to the instruction specified by an integer.
function instruction_vect_to_instruction_int( inst_vect::Vector{MyInt}, numregisters::Int64, numinputs::Int64,
    funcs::Vector{Func}; nodearity::Int64=2 )
  result = Int64(inst_vect[1]-1)  # Function code
  multiplier = 1
  #println("multiplier: ",multiplier,"  result: ",result)
  multiplier = numregisters
  result *= multiplier
  result += inst_vect[2]-1
  #println("multiplier: ",multiplier,"  result: ",result)
  multiplier = numregisters+numinputs
  for j = 1:nodearity
    result *= multiplier
    result += inst_vect[2+j]-1
    #println("multiplier: ",multiplier,"  result: ",result)
  end
  result + 1    # Convert 1-based from 0-based
end

# Converts an instruction specified by an integer to the instruction specified by a vector.
function instruction_int_to_instruction_vect( inst_int::OutputType, numregisters::Int64, numinputs::Int64,
    funcs::Vector{Func}; nodearity::Int64=2 )
  #println("ii_to_iv inst_int: ",inst_int)
  result = fill(MyInt(0),2+nodearity)
  multiplier = numregisters+numinputs
  inst_int -= 1   # Convert from 1-based to 0-based so that mod and div work
  for j = nodearity:-1:1
    try
      result[j+2] = inst_int % multiplier + 1
    catch
      println("catch inst_int: ",inst_int,"  multiplier: ",multiplier)
      error(" error ")
    end
    inst_int = div(inst_int,multiplier)
    #println("j+2: ",j+2,"  multiplier: ", multiplier,"  inst_int: ",inst_int,"  result[j+2]: ",result[j+2])
  end
  multiplier = numregisters
  result[2] = inst_int % multiplier + 1
  inst_int = div(inst_int,multiplier)
  #println("i: ",2,"  multiplier: ", multiplier,"  inst_int: ",inst_int,"  result[2]: ",result[2])
  multiplier = length(funcs)
  result[1] = inst_int % multiplier + 1
  inst_int = div(inst_int,multiplier)
  #println("i: ",1,"  multiplier: ", multiplier,"  inst_int: ",inst_int,"  result[1]: ",result[1])
  result
end

function instruction_vects_to_instruction_ints( cv::Vector{Vector{MyInt}}, numregisters::Int64, numinputs::Int64, funcs::Vector{Func};
    nodearity::Int64=2 )
  map(x->instruction_vect_to_instruction_int(x, numregisters, numinputs, funcs, nodearity=nodearity ), cv )
end

function instruction_vects_to_instruction_ints( circ::LinCircuit, funcs::Vector{Func}=Func[] )
  if length(funcs) == 0
    funcs = lin_funcs( circ.params.numinputs )
  end
  instruction_vects_to_instruction_ints( circ.circuit_vects, circ.params.numlevelsback, circ.params.numinputs, funcs )
end

function instruction_ints_to_instruction_vects( c_ints::Vector{Int64}, numregisters::Int64, numinputs::Int64, funcs::Vector{Func};
    nodearity::Int64=2 )
  map(x->instruction_int_to_instruction_vect(x, numregisters, numinputs, funcs, nodearity=nodearity ), c_ints )
end

function instruction_ints_to_instruction_vects( c_ints::Vector{Int64}, p::Parameters, funcs::Vector{Func}=Func[];
    nodearity::Int64=2 )
  if length(funcs) == 0
    funcs = lin_funcs( p.numinputs )
  end
  instruction_ints_to_instruction_vects( c_ints, p.numlevelsback, p.numinputs, funcs )
end

#=
function circuit_vect_to_circuit_ints( c::LinCircuit, funcs::Vector{Func}=Func[] )
  if length(funcs) == 0
    funcs = lin_funcs( c.params.numinputs )
  end
  circuit_vect_to_circuit_ints( c.circuit_vects, c.params.numinteriors, c.params.numinputs, funcs )
end
=#

function instruction_ints_to_circuit_int( c_ints::Vector{Int64}, p::Parameters, funcs::Vector{Func} )
  numregisters = p.numlevelsback
  multiplier = num_instructions( p, funcs )
  #println("multiplier: ",multiplier)
  vect_to_int( c_ints, fill(multiplier,length(c_ints)))
end

function circuit_int_to_instruction_ints( c_int::Int128, p::Parameters, funcs::Vector{Func} )
  numinstructions = p.numinteriors
  multiplier = num_instructions( p, funcs )
  #println("multiplier: ",multiplier)
  int_to_vect( c_int, fill( multiplier, numinstructions ))
end

function circuit_to_circuit_int( circ::LinCircuit, funcs::Vector{Func} )
  instruction_ints_to_circuit_int( instruction_vects_to_instruction_ints( circ, funcs ), circ.params, funcs )
end

function circuit_int_to_circuit( c_int::Int128, p::Parameters, funcs::Vector{Func} )
  i_vects = instruction_ints_to_instruction_vects( circuit_int_to_instruction_ints( c_int, p, funcs), p, funcs )
  LinCircuit( i_vects, p )
end

# Random instruction
# Assumes gates are arity 2
function rand_ivect( numregisters::Int64, numinputs::Int64, funcs::Vector{Func}; nodearity::Int64=2 )
  MyInt[ rand(1:length(funcs)), rand(1:numregisters),  rand(1:(numregisters+numinputs)), rand(1:(numregisters+numinputs)) ]
end

# Random instruction
# Assumes gates are arity 2 which means that instructions have 2 components
function rand_ivect( p::Parameters, funcs::Vector{Func}; nodearity::Int64=2 )
  MyInt[ rand(1:length(funcs)), rand(1:p.numlevelsback),  rand(1:(p.numlevelsback+p.numinputs)),
    rand(1:(p.numlevelsback+p.numinputs)) ]
end

# Random circuit as a list of instruction integers
# Assumes gates are arity 2
function rand_lcircuit( n_instructions::Int64, numregisters::Int64, numinputs::Int64, funcs::Vector{Func} )
  p = Parameters( numinputs, 1, n_instructions, numregisters )
  LinCircuit( [ rand_ivect( numregisters, numinputs, funcs ) for _ = 1:n_instructions ], p )
end

# Random circuit as a list of instruction integers
# Assumes gates are arity 2
# p.numinteriors is equivalent to numinstructions
# p.numlevelsback is equivalent to numregisters
function rand_lcircuit( p::Parameters, funcs::Vector{Func}=lin_funcs(p.numinputs) )
  LinCircuit( [ rand_ivect( p.numlevelsback, p.numinputs, funcs ) for _ = 1:p.numinteriors], p )
end

function execute_random_circuits( ncircuits::Int64, prog_length::Int64, numregisters::Int64, numinputs::Int64, funcs::Vector{Func} )
  n_instructions = num_instructions( numregisters, numinputs, funcs )
  for i = 1:ncircuits
    circuit_ints = rand(0:(n_instructions-1),prog_length)
    execute_lcircuit( circuit_ints, numregisters, numinputs, funcs )
  end
end

# Mutates the instruction instruction_vect, returns the instruction vector of the mutated instruction
# mut_loc is the mutate location.
# Returns the mutated instruction as a vector of MyInts
# The mutated instruction is guaranteed to be different than instruction_vect.
function mutate_instruction( instruction_vect::Vector{MyInt}, numregisters::Int64, numinputs::Int64, funcs::Vector{Func}, mut_loc::Int64;
      nodearity::Int64=2 )
  lfuncs = length(funcs)
  result_vect = deepcopy(instruction_vect)
  num_mutate_locations = ((lfuncs>1) ? (lfuncs-1) : 0) + (numregisters-1) + nodearity*(numregisters+numinputs-1)
  @assert 1 <= mut_loc && mut_loc <= num_mutate_locations
  if mut_loc <= lfuncs-1
    result_vect[1] = (mut_loc<instruction_vect[1]) ? mut_loc : (mut_loc+1)
    #println("result_vect[1]: ",result_vect[1])
  elseif mut_loc <= (lfuncs-1) + (numregisters-1)
    loc = mut_loc - (lfuncs-1)
    result_vect[2] = (loc<instruction_vect[2]) ? loc : (loc+1)
    #println("result_vect[2]: ",result_vect[2])
  else
    for i = 1:nodearity
      if mut_loc <= (lfuncs-1) + (numregisters-1) + (numregisters+numinputs-1)*i
        loc = mut_loc - (lfuncs-1) - (numregisters-1) - (numregisters+numinputs-1)*(i-1)
        result_vect[2+i] = (loc<instruction_vect[2+i]) ? loc : (loc+1)
        #println("result_vect[",2+i,"]: ",result_vect[2+i])
        break
      end
    end
  end
  result_vect
end

function mutate_instruction( instruction_vect::Vector{MyInt}, p::Parameters, funcs::Vector{Func}, mut_loc::Int64;
      nodearity::Int64=2 )
  mutate_instruction( instruction_vect, p.numlevelsback, p.numinputs, funcs, mut_loc, nodearity=p.nodearity )
end

# Mutates the instruction instruction_vect at a random location, returns the instruction vector of the mutated instruction
# Returns the mutated instruction as a vector of MyInts
function mutate_instruction( instruction_vect::Vector{MyInt}, numregisters::Int64, numinputs::Int64, funcs::Vector{Func}; nodearity::Int64=2 )
  lfuncs = length(funcs)
  num_mutate_locations = ((lfuncs>1) ? (lfuncs-1) : 0) + (numregisters-1) + nodearity*(numregisters+numinputs-1)
  mut_loc = rand(1:num_mutate_locations)
  #println("num_mutate_locations: ",num_mutate_locations,"  mut_loc: ",mut_loc)
  mutate_instruction( instruction_vect, numregisters, numinputs, funcs, mut_loc, nodearity=nodearity )
end

function mutate_instruction( instruction_vect::Vector{MyInt}, p::Parameters, funcs::Vector{Func}; nodearity::Int64=2 )
  mutate_instruction( instruction_vect, p.numlevelsback, p.numinputs, funcs, nodearity=p.nodearity )
end

function mutate_instruction_all( instruction_vect::Vector{MyInt}, numregisters::Int64, numinputs::Int64, funcs::Vector{Func};
      nodearity::Int64=2 )
  lenfuncs = length(funcs)
  num_mutate_locations = ((lenfuncs>1) ? (lenfuncs-1) : 0) + (numregisters-1) + nodearity*(numregisters+numinputs-1)
  instructions = Vector{MyInt}[]
  for mut_loc = 1:num_mutate_locations
    new_instruction = mutate_instruction( instruction_vect, numregisters, numinputs, funcs, mut_loc, nodearity=nodearity )
    push!(instructions,new_instruction)
  end
  instructions
end

function mutate_circuit!( circuit_vect::Vector{Vector{MyInt}}, numregisters::Int64, numinputs::Int64, funcs::Vector{Func}; nodearity::Int64=2 )
  numinstructions = length(circuit_vect)
  ind = rand(1:numinstructions)
  #println("ind: ",ind)
  circuit_vect[ind] = mutate_instruction( circuit_vect[ind], numregisters::Int64, numinputs::Int64, funcs::Vector{Func}; nodearity=nodearity )
  circuit_vect
end

function mutate_circuit!( circuit_vect::Vector{Vector{MyInt}}, p::Parameters, funcs::Vector{Func}; nodearity::Int64=2 )
  mutate_circuit!( circuit_vect, p.numlevelsback, p.numinputs, funcs, nodearity=p.nodearity )
end

function mutate_circuit!( circuit::LinCircuit, funcs::Vector{Func}; nodearity::Int64=2 )
  #LinCircuit( params, mutate_circuit!( circuit.circuit_vects, circuit.params, funcs ) )
  LinCircuit( mutate_circuit!( circuit.circuit_vects, circuit.params, funcs ), circuit.params )
end

function mutate_circuit_all( circuit_vect::Vector{Vector{MyInt}}, numregisters::Int64, numinputs::Int64, funcs::Vector{Func}; nodearity::Int64=2 )
  lenfuncs = length(funcs)
  result = Vector{Vector{MyInt}}[]
  for ind = 1:length(circuit_vect)
    mutated_instructions=mutate_instruction_all(circuit_vect[ind],numregisters,numinputs,funcs,nodearity=nodearity)
    for mi in mutated_instructions
      mutated_circuit = deepcopy(circuit_vect)
      mutated_circuit[ind] = mi
      push!(result,mutated_circuit)
    end
  end
  result
end

function mutate_circuit_all( circuit_vect::Vector{Vector{MyInt}}, p::Parameters, funcs::Vector{Func}; nodearity::Int64=2 )
  mutate_circuit_all( circuit_vect, p.numlevelsback, p.numinputs, funcs, nodearity=p.nodearity )
end

function mutate_all( circuit::LinCircuit, lfuncs::Vector{Func}=Func[]; output_outputs::Bool=true, nodearity::Int64=2 )
  mca = mutate_circuit_all( circuit, lfuncs, nodearity=nodearity )
  map( x->output_values(LinCircuit(x,circuit.params)), mca )
end

function mutate_circuit_all( circuit::LinCircuit, funcs::Vector{Func}=Func[]; nodearity::Int64=2 )
  if length(funcs) == 0
    funcs = lin_funcs( circuit.params.numinputs )
  end
  mutate_circuit_all( circuit.circuit_vects, circuit.params, funcs )
end

function print_circuit( f::IO, circuit::LinCircuit, funcs::Vector{Func} )
  print(f,"(")
  i = 1
  for lc in circuit.circuit_vects
    print(f,"R$(lc[1])=$(funcs[lc[2]].name)(R$(lc[3]),R$(lc[4]))")
    if i < length(circuit.circuit_vects) print(f,",") end
    i += 1
  end
  println(f,")")
end

function print_circuit( circuit::LinCircuit, funcs::Vector{Func} )
  print_circuit( Base.stdout, circuit, funcs )
end

function instruction_distance( instruction_vect1::Vector{MyInt}, instruction_vect2::Vector{MyInt} )
  @assert length(instruction_vect1) == length(instruction_vect2)
  count_diffs = 0
  for i in 1:length(instruction_vect1)
    count_diffs += instruction_vect1[i] == instruction_vect2[i] ? 0 : 1
  end
  count_diffs
end

function circuit_distance( circuit_vect1::Vector{Vector{MyInt}},circuit_vect2::Vector{Vector{MyInt}} )
  @assert length(circuit_vect1) == length(circuit_vect2)
  sum_distances = 0
  for i in 1:length(circuit_vect1)
    sum_distances += instruction_distance( circuit_vect1[i], circuit_vect2[i] )
  end
  sum_distances
end


text:
function infidelity!(u,δt,Ut,Hd,Hc,H,U,X,D,V,u_last)
    # Calculate forward propagators
    calc_fprops!(U,X,D,V,u,δt,Hd,Hc,H,u_last)
    Uf = X[end]; N² = size(Ut,1)^2
    # Calculate infidelity:
    # fₑ = 1 - Φ where Φ = |⟨Ut,Uf⟩|²/N²
    return 1 - abs2(inner(Ut,Uf))/N²
end

function infidelityprime!(fp,u,δt,Ut,Hd,Hc,H,A,U,X,D,V,P,u_last)
    # Calculate forward and backward propagators
    calc_fprops!(U,X,D,V,u,δt,Hd,Hc,H,u_last)
    calc_bprops!(P,U,Ut)
    n = length(U); m = length(Hc); N² = size(Ut,1)^2
    # Calculate exact derivative of the infidelity function:
    # ∂Φ/∂uₖⱼ  = ⟨Pⱼ,∂Uⱼ/∂uₖⱼ*Xⱼ₋₁⟩⟨Xⱼ,Pⱼ⟩/N² + c.c.
    # ∂fₑ/∂uₖⱼ = -∂Φ/∂uₖⱼ
    #          = -2*Re(⟨Pⱼ,Jₖⱼ*Xⱼ₋₁⟩⟨Uf,Ut⟩)/N² where Jₖⱼ = ∂Uⱼ/∂uₖⱼ
    Jkj = similar(A)
    cisDj = Vector{ComplexF64}(length(D[1]))
    a = inner(X[end],Ut)
    for j = 1:n
        cisDj .= cis.(D[j])
        for k = 1:m
            Jmat!(Jkj,Hc[k],cisDj,D[j],V[j],δt,A)
            j==1 ? copy!(A,Jkj) : mul!(A,Jkj,X[j-1])
            fp[(k-1)*n+j] = -2*real(inner(P[j],A)*a)/N²
        end
    end
    return fp
end

function infidelityprimeapprox!(fp,u,δt,Ut,Hd,Hc,H,A,U,X,D,V,P,u_last)
    # Calculate forward and backward propagators
    calc_fprops!(U,X,D,V,u,δt,Hd,Hc,H,u_last)
    calc_bprops!(P,U,Ut)
    n = length(U); m = length(Hc); N² = size(Ut,1)^2
    # Calculate approximate derivative of fidelity error function:
    # ∂Φ/∂uₖⱼ  = ⟨Pⱼ,∂Uⱼ/∂uₖⱼ*Xⱼ₋₁⟩⟨Xⱼ,Pⱼ⟩/N² + c.c.
    # ∂fₑ/∂uₖⱼ = -∂Φ/∂uₖⱼ
    #          ≈ 2*Re(⟨Pⱼ,iδt*Hₖ*Xⱼ⟩⟨Uf,Ut⟩)/N²
    a = inner(X[end],Ut)
    for j = 1:n, k = 1:m
        mul!(A,Hc[k],X[j])
        fp[(k-1)*n+j] = -2δt*imag(inner(P[j],A)*a)/N²
    end
    return fp
end

function gen_opt_fun(Ut::Operator,Hd::Operator,Hc::Vector{<:Operator},t::Real,n::Integer)
    N = prod(dims(Hd))
    m = length(Hc)
    # Make sure we pass dense operators
    Ut_d = full(Ut)
    Hd_d = full(Hd)
    Hc_d = full.(Hc)
    # Generate cache for various objects
    H = promote_type(typeof(Hd_d),eltype(Hc_d))(N,N)
    A = Matrix{ComplexF64}(N,N)
    U = [similar(A) for i=1:n]
    X = deepcopy(U)
    P = deepcopy(U)
    # For the exact derivative we need to store eigenvectors and eigenvalues
    V = [similar(H) for i=1:n]
    D = [Vector{Float64}(N) for i=1:n] # eigenvalues are always real
    # Storage for last control ampitudes, NaN for first run
    u_last = fill(NaN64,m*n)
    # Create optimization function object
    f = (u) -> infidelity!(u,t/n,Ut_d,Hd_d,Hc_d,H,U,X,D,V,u_last)
    g! = (fp,u) -> infidelityprime!(fp,u,t/n,Ut_d,Hd_d,Hc_d,H,A,U,X,D,V,P,u_last)
    #g! = (fp,u) -> infidelityprimeapprox!(fp,u,t/n,Ut_d,Hd_d,Hc_d,H,A,U,X,D,V,P,u_last)
    # Return a OnceDifferentiable object with appropriate seed
    return OnceDifferentiable(f,g!,zeros(m*n)),X[end]
end

function grape(Ut::Operator,Hd::Operator,Hc::Vector{<:Operator},u_init,t::Real,n::Integer)
    m = length(Hc)
    n == size(u_init,1) || throw(ArgumentError("control amplitude matrix not consistent with number of timesteps"))
    m == size(u_init,2) || throw(ArgumentError("control amplitude matrix not consistent with number of control hamiltonians"))
    # Generate OnceDifferentiable object, and get reference to final U
    od, Uf = gen_opt_fun(Ut,Hd,Hc,t,n)
    # Optimization options
    #opt = Optim.Options(g_tol = 1E-9)
    # Run optimization
    res = optimize(od,vec(transpose(u_init)),ConjugateGradient())
    #res = optimize(od,vec(transpose(u_init)))
    # Reshape final control amplitude matrix
    uf = transpose(reshape(Optim.minimizer(res),m,n))
    return uf, Uf, res
end

function opt_2Q_QFT1(n)
    Hd = 0.5 * (σx⊗σx + σy⊗σy + σz⊗σz)
    Hc = [0.5*σx⊗σ0, 0.5*σy⊗σ0, 0.5*σ0⊗σx, 0.5*σ0⊗σy]
    t = 6.0
    u_init = zeros(n,4)
    #u_init = rand(n,4) .- 0.5
    Ut = Operator([0.5  0.5    0.5  0.5
                   0.5  0.5im -0.5 -0.5im
                   0.5 -0.5    0.5 -0.5
                   0.5 -0.5im -0.5  0.5im],
                  (2,2))
    grape(Ut,Hd,Hc,u_init,t,n)
end


text:
"""
    rpolyr(x, β, θ)

Generate a random integer `Y` such that `P(Y≤j) = g^{-1}(θ[j]-x'β)`.
`θ` has to be monotone `θ[1] ≤ ... ≤ θ[J-1]`.
"""
function rpolr(
    x::AbstractVector,
    β::AbstractVector,
    θ::AbstractVector,
    link::GLM.Link
    )
    # number of categories
    J = length(θ) + 1
    # check monotonicity of θj
    issorted(θ) || throw(ArgumentError("θ[j] should be nondecreasing."))
    # generate category according to cumulative probability
    iprod = dot(x, β)
    q = rand()
    for j in 1:J-1
        ηj = θ[j] - iprod
        cumprobj = GLM.linkinv(link, ηj)
        q ≤ cumprobj && (return j)
    end
    return J
end

"""
    rpolr(X, β, θ)

Generate a vector of random integers `Y` such that
`P(Y[i]≤j) = g^{-1}(θ[j]-X[:,j]'β)`.
`θ` has to be monotone `θ[1] ≤ ... ≤ θ[J-1]`.
"""
function rpolr(
    X::AbstractMatrix,
    β::AbstractVector,
    θ::AbstractVector,
    link::GLM.Link
    )
    @views Y = [rpolr(X[i, :], β, θ, link) for i in 1:size(X, 1)]
end


text:
module EMI

module Draw
include("emi_draw.jl")

export Point
export Curve, Line, CircleArc, EllipseArc
export Shape, Circle, Ellipse, Rectangle, Square, ClosedPolygon, NGon, NGonR, Loop
export BoundingBox
export Canvas, set_bbox!

end  # Draw

module Gmsh
using EMI.Draw
include("emi_gmsh.jl")

export gmsh_script, mesh

end  # Gmsh

end # EMI

max_stars_repo_name: MiroK/emi-cylinders

max_issues_repo_name: MiroK/emi-cylinders

max_forks_repo_name: MiroK/emi-cylinders


text:
using GZip

mutable struct GZBufferedStream <: IO
    io::GZipStream
    buf::Vector{UInt8}
    len::Int
    ptr::Int

    function GZBufferedStream(io::GZipStream)
        buf = Array(UInt8, io.buf_size)

        len = ccall((:gzread, GZip._zlib), Int32,
            (Ptr{Cvoid}, Ptr{Cvoid}, UInt32), io.gz_file, buf, io.buf_size)

        new(io, buf, len, 1)
    end
end

Base.close(io::GZBufferedStream) = close(io.io)

@inline function Base.read(io::GZBufferedStream, ::Type{UInt8})
    c = io.buf[io.ptr]
    io.ptr += 1

    if io.ptr == io.len+1 #No more data
        io.len = ccall((:gzread, GZip._zlib), Int32,
            (Ptr{Cvoid}, Ptr{Cvoid}, UInt32), io.io.gz_file, io.buf, io.io.buf_size)
        io.ptr = 1
    end
    c
end

@inline Base.eof(io::GZBufferedStream) = io.len == 0

#function bench()
    #io = GZBufferedStream(GZip.open("random.gz", "rb"))

    #thischar = 0x00
    #n = 0
    #while !eof(io)
        #thischar = read(io, UInt8)
        #n += 1
    #end
    #close(io)

    #println(n)
    #thischar
#end

#bench()

##random contains 10^8 rand(UInt8)s
#@time run(`gunzip -k -f random.gz`) #255 ms
#@time bench() #168 ms


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    VF2Matcher,
    isomorphisms, is_isomorphic,
    nodesubgraph_isomorphisms, nodesubgraph_is_isomorphic,
    edgesubgraph_isomorphisms, edgesubgraph_is_isomorphic,
    subgraph_monomorphisms, subgraph_is_monomorphic

"""
    VF2Matcher{T1<:UndirectedGraph,T2<:UndirectedGraph}(
        G::T1, H::T2, matchtype::Symbol; kwargs...
    )

Lazy iterator that generate all isomorphism mappings between `G` and `H`.

`matchtype` should be one of the followings
  - `:isomorphic`: G is isomorphic to H
  - `:subgraph_isomorphic`: a node-induced subgraph of G is isomorphic to H
  - `:monomorphic`: a subgraph of G is monomorphic to H

# Options

- `nodematcher::Function`: a function for semantic node attribute matching (default: (a, b) -> true)
- `edgematcher::Function`: a function for semantic edge attribute matching (default: (a, b) -> true)
- `mandatory::Dict{Int,Int}`: mandatory node mapping (if matchtype=:edgeinduced, edge mapping)
- `forbidden::Dict{Int,Int}`: forbidden node mapping (if matchtype=:edgeinduced, edge mapping)
- `timeout::Union{Int,Nothing}`: if specified, abort vf2 calculation when the time reached and return empty iterator (default: 10 seconds)

"""
mutable struct VF2Matcher{T1<:UndirectedGraph,T2<:UndirectedGraph}
    G::T1
    H::T2
    matchtype::Symbol
    nodematcher::Function
    edgematcher::Function
    mandatory::Dict{Int,Int}
    forbidden::Dict{Int,Int}
    expire::Union{UInt64,Nothing} # UInt64, nanoseconds

    stack::Vector{Tuple{Symbol,Int,Int}}
    currentpair::Union{Tuple{Int,Int},Nothing}
    g_core::Dict{Int,Int}
    h_core::Dict{Int,Int}
    g_term::Dict{Int,Int}
    h_term::Dict{Int,Int}
    timedout::Bool

    function VF2Matcher{T1,T2}(
                G::T1, H::T2, matchtype::Symbol;
                nodematcher=(a, b)->true, edgematcher=(a, b)->true,
                mandatory=Dict(), forbidden=Dict(), timeout=nothing
            ) where {T1<:UndirectedGraph,T2<:UndirectedGraph}
        expire = (timeout === nothing ?
            nothing : (time_ns() + timeout * 1_000_000_000)::UInt64)
        return new(
            G, H, matchtype, nodematcher, edgematcher,
            mandatory, forbidden, expire,
            [], nothing, Dict(), Dict(), Dict(), Dict(), false
        )
    end
end
VF2Matcher(
    G::UndirectedGraph, H::UndirectedGraph, matchtype::Symbol; kwargs...
) = VF2Matcher{typeof(G),typeof(H)}(G, H, matchtype; kwargs...)

function candidatepairs(iter::VF2Matcher)
    # Mandatory pair
    md = setdiff(keys(iter.mandatory), keys(iter.g_core))
    if !isempty(md)
        n = pop!(md)
        return [(n, iter.mandatory[n])]
    end

    pairs = Tuple{Int,Int}[]
    g_cand = setdiff(keys(iter.g_term), keys(iter.g_core))
    h_cand = setdiff(keys(iter.h_term), keys(iter.h_core))
    if isempty(g_cand) || isempty(h_cand)
        # New connected component
        g_cand = setdiff(nodeset(iter.G), keys(iter.g_core))
        h_cand = setdiff(nodeset(iter.H), keys(iter.h_core))
    end
    if !isempty(h_cand)
        # improved pivot for monomorphism match with many isolated nodes
        h_sorted = sort!(collect(h_cand))
        h_min = h_sorted[argmax([degree(iter.H, n) for n in h_sorted])]
        for g in g_cand
            # Forbidden pair
            if haskey(iter.forbidden, g) && iter.forbidden[g] == h_min
                continue
            end
            push!(pairs, (g, h_min))
        end
    end
    return pairs
end

function is_feasible(iter::VF2Matcher, g, h)
    # Note: assume simple graph (no self loops and multi edges)
    g_nbrs = adjacencies(iter.G, g)
    h_nbrs = adjacencies(iter.H, h)
    if iter.matchtype !== :monomorphic
        for n in intersect(g_nbrs, keys(iter.g_core))
            if !(iter.g_core[n] in h_nbrs)
                @debug "Infeasible: neighbor connectivity"
                return false
            end
        end
    end
    for n in intersect(h_nbrs, keys(iter.h_core))
        if !(iter.h_core[n] in g_nbrs)
            @debug "Infeasible: neighbor connectivity"
            return false
        end
    end

    if iter.matchtype !== :monomorphic
        g_term_count = length(setdiff(keys(iter.g_term), keys(iter.g_core)))
        h_term_count = length(setdiff(keys(iter.h_term), keys(iter.h_core)))
        if iter.matchtype === :exact && g_term_count != h_term_count
            @debug "Infeasible: terminal set size"
            return false
        elseif iter.matchtype === :subgraph_isomorphic && g_term_count < h_term_count
            @debug "Infeasible: terminal set size"
            return false
        end

        g_new_count = length(setdiff(nodeset(iter.G), keys(iter.g_term)))
        h_new_count = length(setdiff(nodeset(iter.H), keys(iter.h_term)))
        if iter.matchtype === :isomorphic && g_new_count != h_new_count
            @debug "Infeasible: yet unexplored nodes"
            return false
        elseif iter.matchtype === :subgraph_isomorphic && g_new_count < h_new_count
            @debug "Infeasible: yet unexplored nodes"
            return false
        end
    end
    @debug "Syntactically feasible"
    return true
end

function is_semantic_feasible(iter::VF2Matcher, g, h)
    if !iter.nodematcher(g, h)
        @debug "Infeasible: node attribute mismatch"
        return false
    end
    for (inc, adj) in neighbors(iter.G, g)
        haskey(iter.g_core, adj) || continue
        hinc = findedgekey(iter.H, h, iter.g_core[adj])
        iter.matchtype === :monomorphic && hinc === nothing && continue
        if !iter.edgematcher(inc, hinc)
            @debug "Infeasible: edge attribute mismatch"
            return false
        end
    end
    return true
end

function expand!(iter::VF2Matcher, g, h)
    iter.g_core[g] = h
    iter.h_core[h] = g
    depth = length(iter.g_core)
    if !haskey(iter.g_term, g)
        iter.g_term[g] = depth
    end
    if !haskey(iter.h_term, h)
        iter.h_term[h] = depth
    end
    g_nbrset = union([adjacencies(iter.G, n) for n in keys(iter.g_core)]...)
    for n in setdiff(g_nbrset, keys(iter.g_term))
        iter.g_term[n] = depth
    end
    h_nbrset = union([adjacencies(iter.H, n) for n in keys(iter.h_core)]...)
    for n in setdiff(h_nbrset, keys(iter.h_term))
        iter.h_term[n] = depth
    end
    return
end

function restore!(iter::VF2Matcher, g, h)
    depth = length(iter.g_core)
    if g !== nothing && h !== nothing
        delete!(iter.g_core, g)
        delete!(iter.h_core, h)
    end
    for (k, v) in iter.g_term
        v == depth && delete!(iter.g_term, k)
    end
    for (k, v) in iter.h_term
        v == depth && delete!(iter.h_term, k)
    end
end

function Base.iterate(iter::VF2Matcher, state=nothing)
    iter.timedout && return
    if iter.currentpair !== nothing
        g1, h1 = iter.currentpair
        push!(iter.stack, (:restore, g1, h1))
    end
    for (g, h) in candidatepairs(iter)
        @debug "candidates" g h
        is_feasible(iter, g, h) || continue
        is_semantic_feasible(iter, g, h) || continue
        push!(iter.stack, (:expand, g, h))
    end
    while !isempty(iter.stack)
        command, g1, h1 = pop!(iter.stack)
        if command === :restore
            @debug "restored" iter.g_core
            restore!(iter, g1, h1)
            continue
        end
        @debug "g_core" iter.g_core
        expand!(iter, g1, h1)
        iter.currentpair = (g1, h1)
        @debug "depth" length(iter.g_core)
        if length(iter.g_core) == nodecount(iter.H)
            @debug "done" iter.g_core
            return (copy(iter.g_core), state)
        elseif iter.expire !== nothing && time_ns() > iter.expire
            iter.timedout = true
            @debug "Timed out"
            return
        end
        push!(iter.stack, (:restore, g1, h1))
        for (g, h) in candidatepairs(iter)
            @debug "candidates" g h
            is_feasible(iter, g, h) || continue
            is_semantic_feasible(iter, g, h) || continue
            push!(iter.stack, (:expand, g, h))
        end
    end
    return  # no match found
end

# TODO: get rid of generics
Base.IteratorSize(::Type{VF2Matcher{T1,T2}}) where {T1<:UndirectedGraph,T2<:UndirectedGraph} = Base.SizeUnknown()
Base.IteratorEltype(::Type{VF2Matcher{T1,T2}}) where {T1<:UndirectedGraph,T2<:UndirectedGraph} = Base.EltypeUnknown()

# Graph isomorphism

"""
    isomorphisms(G::UndirectedGraph, H::UndirectedGraph; kwargs...) -> Iterator

Return an iterator that generate isomorphic mappings between `G` and `H`.
The returned iterator has `ig => ih` pairs that correspond to the indices of matching nodes in `G` and `H`, respectively.
"""
isomorphisms(G, H; kwargs...) = VF2Matcher(G, H, :isomorphic; kwargs...)

"""
    is_isomorphic(G::UndirectedGraph, H::UndirectedGraph; kwargs...) -> Bool

Return whether `G` and `H` are isomorphic.
"""
is_isomorphic(G, H; kwargs...) = !isempty(isomorphisms(G, H; kwargs...))

# Node-induced subgraph isomorphism

"""
    nodesubgraph_isomorphisms(G::UndirectedGraph, H::UndirectedGraph; kwargs...) -> Iterator

Return an iterator that generate isomorphic mappings between `H` and node-induced subgraphs of `G`.
The returned iterator has `ig => ih` pairs that correspond to the indices of matching nodes in `G` and `H`, respectively.
"""
nodesubgraph_isomorphisms(G, H; kwargs...
    ) = VF2Matcher(G, H, :subgraph_isomorphic; kwargs...)

"""
    nodesubgraph_is_isomorphic(G::UndirectedGraph, H::UndirectedGraph; kwargs...) -> Bool

Return whether a node-induced subgraph of `G` is isomorphic to `H`.
"""
nodesubgraph_is_isomorphic(G, H; kwargs...
    ) = !isempty(nodesubgraph_isomorphisms(G, H; kwargs...))

# Edge-induced subgraph isomorphism

"""
    edgesubgraph_isomorphisms(
        G::UndirectedGraph, H::UndirectedGraph;
        nodematcher=(g,h)->true, edgematcher=(g,h)->true,
        kwargs...) -> Iterator

Return an iterator that generate isomorphic mappings between `H` and edge-induced subgraphs of `G`.
The returned iterator has `ig => ih` pairs that correspond to the indices of matching edges in `G` and `H`, respectively.

`nodematcher` and `edgematcher` control the features needed to be counted as a match.

See [`Graph.edgesubgraph`](@ref) to construct the subgraphs that result from the match.
"""
function edgesubgraph_isomorphisms(
        G, H; nodematcher=(g,h)->true, edgematcher=(g,h)->true, kwargs...)
    lg = linegraph(G)
    lh = linegraph(H)
    lgedge = lgedgematcher(lg, lh, nodematcher)
    lgnode = lgnodematcher(lg, lh, nodematcher, edgematcher)
    matcher = VF2Matcher(
        lg, lh, :subgraph_isomorphic;
        nodematcher=lgnode, edgematcher=lgedge, kwargs...)
    return Iterators.filter(matcher) do mapping
        return !delta_y_mismatch(G, H, mapping)
    end
end

"""
    edgesubgraph_is_isomorphic(
        G::UndirectedGraph, H::UndirectedGraph; kwargs...) -> Bool

Return whether an edge-induced subgraph of `G` is isomorphic to `H`.
"""
edgesubgraph_is_isomorphic(G, H; kwargs...
    ) = !isempty(edgesubgraph_isomorphisms(G, H; kwargs...))

# Subgraph monomorphism

"""
    subgraph_monomorphisms(
        G::UndirectedGraph, H::UndirectedGraph; kwargs...) -> Iterator

Generate monomorphism mappings between `H` and subgraphs of `G`.
The returned iterator has `ig => ih` pairs that correspond to the indices of matching nodes in `G` and `H`, respectively.
"""
subgraph_monomorphisms(G, H; kwargs...) = VF2Matcher(G, H, :monomorphic; kwargs...)

"""
    subgraph_is_monomorphic(
        G::UndirectedGraph, H::UndirectedGraph; kwargs...) -> Bool

Return whether a subgraph of `G` is monomorphic to `H`.
"""
subgraph_is_monomorphic(G, H; kwargs...
    ) = !isempty(subgraph_monomorphisms(G, H; kwargs...))


text:
module WTA
using SpikingNeuralNets, SpikingNeuralNets.SNNs, SpikingNeuralNets.GraphUtils, SpikingNeuralNets.SNNPlot
using LightGraphs, SimpleWeightedGraphs

function wta_layout(g::AbstractGraph)
    N = round(Int, (nv(g) - 2) / 2) # number of inputs/outputs
    xloc = Float64[fill(1, N)..., fill(2, 2)..., fill(3, N)...]
    yloc = Float64[((1:N) .- N/2)..., -1, 1, ((1:N) .- N/2)...]
    return xloc, yloc
end

function two_inhibitor_snn(n::Integer, γ::Real)
    graph = SimpleWeightedDiGraph(2n+2)
    for i in 1:n
        yi = n+2+i
        add_edge!(graph, i, yi, 3γ)      # inputs -> outputs
        add_edge!(graph, yi, yi, 2γ)
        add_edge!(graph, yi, n+1, γ)
        add_edge!(graph, yi, n+2, γ)
        add_edge!(graph, n+1, yi, -γ)
        add_edge!(graph, n+2, yi, -γ)
    end

    bias = [fill(0, n)..., -γ/2, -3γ/2, fill(-3γ, n)...]
    return SNN{Float64, Int}(graph; bias=bias)
end

function run(input::Vector{T}, γ::Real, args...; iterations=100, kwargs...) where {T<:Integer}
    N = length(input)
    snn = two_inhibitor_snn(N, γ)
    step!(snn, iterations, args...; input=input, transfer=sigmoid, f=fire, kwargs...)
    display(rasterplot(spiketrain(snn, iterations)))
end
end


text:
"""
    These functions read particles in a given range over multiple files, based
    on the relevant peano-hilbert keys.
    This code is based on read_particles_in_box.pro by Dr. habil. Klaus Dolag.
"""

using Dates
using Base.Threads
#using ThreadsX

"""
    Functionality to read the header of the key files
"""

struct KeyHeader
    nkeys_file::Vector{Int32}
    domain_corners::Vector{Float64}
    domain_fac::Float64
    bits::Int32
    nkeys_total::Vector{UInt32}
    nkeys_total_highword::Vector{UInt32}
end

function read_keyheader(filename::String)

    f = open(filename)
    blocksize = read(f, Int32)

    if blocksize == 8
        swap = 0
        snap_format = 2
    elseif blocksize == 256
        swap = 0
        snap_format = 1
    else
        blocksize = bswap(blocksize)
        if blocksize == 8
            swap = 1
            snap_format = 2
        elseif blocksize == 256
            swap = 1
            snap_format = 1
        else
            error("incorrect file format encountered when reading header of $filename")
        end
    end

    if snap_format == 2
        seek(f, 16)
        skip_line = read(f, Int32)
    end

    nkeys = read!(f, Array{Int32,1}(undef,6))
    domain_corners = read!(f, Array{Float64,1}(undef,3))
    domain_fac = read(f, Float64)

    read!(f, Array{Float64,1}(undef,2))
    read(f, Float64)
    read(f, Float64)
    read(f, Int32)

    bits = read(f, Int32)
    nkeys_total = read!(f, Array{UInt32,1}(undef,6))

    read(f, Int32)
    read(f, Int32)
    read(f, Float64)
    read(f, Float64)
    read(f, Float64)
    read(f, Float64)
    read(f, Int32)
    read(f, Int32)

    nkeys_total_highword = read!(f, Array{UInt32,1}(undef,6))

    close(f)

    return KeyHeader( nkeys, domain_corners, domain_fac, bits, nkeys_total, nkeys_total_highword )
end

"""
    Helper function for Peano Hilbert stuff
"""
@inline function get_int_pos(pos::AbstractFloat, domain_corner::Float64, domain_fac::Float64 )
    return Int64( floor(( pos - domain_corner ) * domain_fac) )
end

"""
    Lookup tables for peano hilbert keys
"""
const global rottable3 = [36  28  25  27  10  10  25  27;
                          29  11  24  24  37  11  26  26;
                           8   8  25  27  30  38  25  27;
                           9  39  24  24   9  31  26  26;
                          40  24  44  32  40   6  44   6;
                          25   7  33   7  41  41  45  45;
                           4  42   4  46  26  42  34  46;
                          43  43  47  47   5  27   5  35;
                          33  35  36  28  33  35   2   2;
                          32  32  29   3  34  34  37   3;
                          33  35   0   0  33  35  30  38;
                          32  32   1  39  34  34   1  31;
                          24  42  32  46  14  42  14  46;
                          43  43  47  47  25  15  33  15;
                          40  12  44  12  40  26  44  34;
                          13  27  13  35  41  41  45  45;
                          28  41  28  22  38  43  38  22;
                          42  40  23  23  29  39  29  39;
                          41  36  20  36  43  30  20  30;
                          37  31  37  31  42  40  21  21;
                          28  18  28  45  38  18  38  47;
                          19  19  46  44  29  39  29  39;
                          16  36  45  36  16  30  47  30;
                          37  31  37  31  17  17  46  44;
                          12   4   1   3  34  34   1   3;
                           5  35   0   0  13  35   2   2;
                          32  32   1   3   6  14   1   3;
                          33  15   0   0  33   7   2   2;
                          16   0  20   8  16  30  20  30;
                           1  31   9  31  17  17  21  21;
                          28  18  28  22   2  18  10  22;
                          19  19  23  23  29   3  29  11;
                           9  11  12   4   9  11  26  26;
                           8   8   5  27  10  10  13  27;
                           9  11  24  24   9  11   6  14;
                           8   8  25  15  10  10  25   7;
                           0  18   8  22  38  18  38  22;
                          19  19  23  23   1  39   9  39;
                          16  36  20  36  16   2  20  10;
                          37   3  37  11  17  17  21  21;
                           4  17   4  46  14  19  14  46;
                          18  16  47  47   5  15   5  15;
                          17  12  44  12  19  6   44   6;
                          13   7  13   7  18  16  45  45;
                           4  42   4  21  14  42  14  23;
                          43  43  22  20   5  15   5  15;
                          40  12  21  12  40   6  23   6;
                          13   7  13   7  41  41  22  20 ]

const global subpix3 = [0  7  1  6  3  4  2  5 ;
                        7  4  6  5  0  3  1  2 ;
                        4  3  5  2  7  0  6  1 ;
                        3  0  2  1  4  7  5  6 ;
                        1  0  6  7  2  3  5  4 ;
                        0  3  7  4  1  2  6  5 ;
                        3  2  4  5  0  1  7  6 ;
                        2  1  5  6  3  0  4  7 ;
                        6  1  7  0  5  2  4  3 ;
                        1  2  0  3  6  5  7  4 ;
                        2  5  3  4  1  6  0  7 ;
                        5  6  4  7  2  1  3  0 ;
                        7  6  0  1  4  5  3  2 ;
                        6  5  1  2  7  4  0  3 ;
                        5  4  2  3  6  7  1  0 ;
                        4  7  3  0  5  6  2  1 ;
                        6  7  5  4  1  0  2  3 ;
                        7  0  4  3  6  1  5  2 ;
                        0  1  3  2  7  6  4  5 ;
                        1  6  2  5  0  7  3  4 ;
                        2  3  1  0  5  4  6  7 ;
                        3  4  0  7  2  5  1  6 ;
                        4  5  7  6  3  2  0  1 ;
                        5  2  6  1  4  3  7  0 ;
                        7  0  6  1  4  3  5  2 ;
                        0  3  1  2  7  4  6  5 ;
                        3  4  2  5  0  7  1  6 ;
                        4  7  5  6  3  0  2  1 ;
                        6  7  1  0  5  4  2  3 ;
                        7  4  0  3  6  5  1  2 ;
                        4  5  3  2  7  6  0  1 ;
                        5  6  2  1  4  7  3  0 ;
                        1  6  0  7  2  5  3  4 ;
                        6  5  7  4  1  2  0  3 ;
                        5  2  4  3  6  1  7  0 ;
                        2  1  3  0  5  6  4  7 ;
                        0  1  7  6  3  2  4  5 ;
                        1  2  6  5  0  3  7  4 ;
                        2  3  5  4  1  0  6  7 ;
                        3  0  4  7  2  1  5  6 ;
                        1  0  2  3  6  7  5  4 ;
                        0  7  3  4  1  6  2  5 ;
                        7  6  4  5  0  1  3  2 ;
                        6  1  5  2  7  0  4  3 ;
                        5  4  6  7  2  3  1  0 ;
                        4  3  7  0  5  2  6  1 ;
                        3  2  0  1  4  5  7  6 ;
                        2  5  1  6  3  4  0  7 ]

# ; This function computes a Peano-Hilbert key for an integer triplet (x,y,z),
# ; with x,y,z in the range between 0 and 2^bits-1.
function peano_hilbert_key(bits::Integer, x::Integer, y::Integer, z::Integer)

    rotation = 1
    key = 0
    mask = 1 << (bits-1)

    for imask = 1:bits

        pix = (
                (( x & mask) > 0 ? 4 : 0)  +
                (( y & mask) > 0 ? 2 : 0)  +
                (( z & mask) > 0 ? 1 : 0)
              ) + 1

        key = key << 3

        key = key | subpix3[rotation, pix]

        rotation = rottable3[rotation, pix] + 1

        mask = mask >> 1

    end

    return key
end

function read_key_index(file_key_index::String)

    finfo = stat(file_key_index)

    f = open(file_key_index)
    n = read(f, Int32)

    if finfo.size == 4+n*8*2+n*4
        int_type = UInt64
    else
        int_type = UInt32
    end

    low_list = read!(f, Array{int_type,1}(undef,n))
    high_list = read!(f, Array{int_type,1}(undef,n))
    file_list = read!(f, Array{UInt32,1}(undef,n))

    close(f)

    return low_list, high_list, file_list
end

function get_index_bounds(ids::Vector{Int}, low_bounds::Vector{UInt}, high_bounds::Vector{UInt})

    nids = length(ids)
    nbounds = length(low_bounds)

    ind_all = zeros(Int64, nids)

    icountall = 1

    icountids    = 1
    icountbounds = 1

    lend = false

    while !lend

        if low_bounds[icountbounds] <= ids[icountids] <= high_bounds[icountbounds]

            ind_all[icountall] = icountbounds
            icountall += 1

            lend2 = false
            while !lend2
                if ids[icountids] > high_bounds[icountbounds]
                    lend2 = true
                else # ids[icountids] > high_bounds[icountbounds]
                    icountids += 1

                    if icountids >= nids
                        lend2 = true
                    end # icountids >= nids
                end # ids[icountids] > high_bounds[icountbounds]
            end # while !lend2

            icountbounds += 1

        else # low_bounds[icountbounds] <= ids[icountids] <= high_bounds[icountbounds]

            if ids[icountids] < low_bounds[icountbounds]

                icountids += 1

                if icountids <= nids

                    lend2 = false
                    while !lend2

                        if ids[icountids] >= low_bounds[icountbounds]
                            lend2 = true
                        else # ids[icountids] >= low_bounds[icountbounds]
                            icountids += 1

                            if icountids >= nids
                                lend2 = true
                            end
                        end # if ids[icountids] >= low_bounds[icountbounds]

                    end # while !lend2

                else # if icountids < nids

                    if ids[icountids] > high_bounds[icountbounds]

                        icountbounds += 1

                        if icountbounds < nbounds

                            lend2 = false

                            while !lend2

                                if ids[icountids] <= high_bounds[icountbounds]
                                    lend2 = true
                                else # if ids[icountids] <= high_bounds[icountbounds]
                                    icountbounds += 1
                                    if icountbounds >= nbounds
                                        lend2 = true
                                    end
                                end # if ids[icountids] <= high_bounds[icountbounds]
                            end # while !lend2
                        end # if icountbounds < nbounds

                    end # if ids[icountids] > high_bounds[icountbounds]

                end # if icountids < nids

            else # if ids[icountids] < low_bounds[icountbounds]
                icountbounds += 1

                if icountbounds < nbounds
                    lend2 = false

                    while !lend2
                        if ids[icountids] <= high_bounds[icountbounds]
                            lend2 = true
                        else # ids[icountids] <= high_bounds[icountbounds]
                            icountbounds += 1
                            if icountbounds >= nbounds
                                lend2 = true
                            end # icountbounds >= nbounds
                        end # ids[icountids] <= high_bounds[icountbounds]
                    end # while !lend2

                end # icountbounds < nbounds

            end # if ids[icountids] < low_bounds[icountbounds]

        end # if low_bounds[icountbounds] <= ids[icountids] <= high_bounds[icountbounds]

        if icountids >= nids
            lend = true
        end
        if icountbounds >= nbounds
            lend = true
        end

        #
        # println("icountids     = ", icountids)
        # println("icountbounds  = ", icountbounds)

    end # while !lend

    if icountall > 1
        ind_out = ind_all[1:icountall]
        if ind_out[end] == 0
            return ind_out[1:end-1]
        else
            return ind_out
        end
    else
        return -1
    end
end

"""
This should be fixed! It's A LOT faster than the simpler version!
"""
function find_files_for_keys_old(filebase::String, nfiles::Integer, keylist::Vector{UInt})

    file_key_index = filebase * ".key.index"

    # if index file does not exist all key files need to be read
    if !isfile(file_key_index)
        return collect(0:nfiles-1)
    end

    # get the data from the index file
    low_list, high_list, file_list = read_key_index(file_key_index)

    # sort the keys
    key_sort = sortperm(keylist)

    index_bounds = get_index_bounds(keylist[key_sort], low_list, high_list)

    file_sort = sort(file_list[index_bounds])

    files = file_list[file_sort]

    return Int64.(unique!(files))
end

function find_files_for_keys(filebase::String, nfiles::Integer, keylist::Vector{Int})

    file_key_index = filebase * ".key.index"

    # if index file does not exist all key files need to be read
    if !isfile(file_key_index)
        return collect(0:nfiles-1)
    end

    # get the data from the index file
    low_list, high_list, file_list = read_key_index(file_key_index)

    mask = falses(length(low_list))

    for key in keylist
        @. mask = mask | ( (key >= low_list ) & ( key <= high_list ))
    end

    return Int64.(unique!(sort!(file_list[mask])))
end

function get_index_list_old(idarr1, idarr2)

    narr1 = length(idarr1)
    narr2 = length(idarr2)

    ind_all   = zeros(Int64, narr1)
    not_arr2t = zeros(Int64, narr1)
    not_arr1t = zeros(Int64, narr2)

    icountall     = 1
    icountarr1    = 1
    icountarr2    = 1
    icountnotarr1 = 1
    icountnotarr2 = 1

    lend = false

    iiarr1 = sortperm(idarr1[:,1])
    iiarr2 = sortperm(idarr2[:,1])

    while !lend

        if idarr2[iiarr2[icountarr2]] == idarr1[iiarr1[icountarr1]]

            ind_all[icountall] = iiarr2[icountarr2]
            icountall  += 1
            icountarr1 += 1
            icountarr2 += 1
        else  # idarr2[iiarr2[icountnotarr2]] == idarr1[iiarr1[icountnotarr1]]
            if idarr2[iiarr2[icountarr2]] < idarr1[iiarr1[icountarr1]]
                not_arr1t[icountnotarr1] = iiarr2[icountarr2]
                icountarr2    += 1
                icountnotarr1 += 1
            else # idarr2[iiarr2[icountnotarr2]] < idarr1[iiarr1[icountnotarr1]]
                not_arr2t[icountnotarr2] = iiarr1[icountarr1]
                icountarr1    += 1
                icountnotarr2 += 1
            end # idarr2[iiarr2[icountnotarr2]] < idarr1[iiarr1[icountnotarr1]]
        end # idarr2[iiarr2[icountnotarr2]] == idarr1[iiarr1[icountnotarr1]]
        if (icountarr2 >= narr2 ) || (icountarr1 >= narr1)
            lend = true
        end
    end # while !lend

    rest = narr1 - icountarr1
    if rest > 0
        not_arr2t[icountnotarr2:icountnotarr2+rest] = iiarr1[icountarr1:icountarr1+rest]
        icountnotarr2 += rest
    end # rest > 0

    if icountall > 1
        ind_out = ind_all[1:icountall]
    # else
    #     error("Not enough data found!")
    end

    return ind_out
end

function get_keylist(h_key::KeyHeader, x0, x1)

    ix0 = zeros(Int, 3)
    ix1 = zeros(Int, 3)
    dix = zeros(Int, 3)

    nkeys = 1

    @inbounds for i = 1:3
        @inbounds ix0[i] = get_int_pos( x0[i], h_key.domain_corners[i], h_key.domain_fac )
        @inbounds ix1[i] = get_int_pos( x1[i], h_key.domain_corners[i], h_key.domain_fac )
        @inbounds dix[i] = ix1[i] - ix0[i] + 1
        @inbounds nkeys *= dix[i]
    end

    keylist = zeros(Int, nkeys)

    i = 1
    for ix = ix0[1]:ix1[1], iy = ix0[2]:ix1[2], iz = ix0[3]:ix1[3]
        @inbounds keylist[i] = peano_hilbert_key(h_key.bits, ix, iy, iz)
        i += 1
    end

    return keylist
end

@inline function get_index_list(keylist::Vector{Int}, keys_in_file)

    dict = Dict((n, i) for (i, n) in enumerate(keys_in_file))
    result = Vector{Int}(undef, length(keylist))
    len = 0

    for k in keylist
        i = get(dict, k, nothing)
        if i !== nothing
            len += 1
            @inbounds result[len] = i
        end
    end
    return resize!(result, len)
end

@inline function get_index_list_serial(keylist::Vector{Int}, keys_in_file)

    dict = Dict((n, i) for (i, n) in enumerate(keys_in_file))
    result = Vector{Int}(undef, length(keylist))
    len = 0

    for k in keylist
        i = get(dict, k, nothing)
        if i !== nothing
            len += 1
            @inbounds result[len] = i
        end
    end
    return resize!(result, len)
end

@inline function get_index_list_parallel(keylist::Vector{Int}, keys_in_file)

    dict = Dict((n, i) for (i, n) in enumerate(keys_in_file))
    result = Vector{Int}(undef, length(keylist))
    #len = 0
    len = Threads.Atomic{Int}(0)

    @threads for k = 1:length(keylist)
        i = get(dict, keylist[k], nothing)
        if i !== nothing
            #len += 1
            Threads.atomic_add!(len, 1)
            @inbounds result[len[]] = i
        end
    end
    return resize!(result, len[])
end

function join_blocks(offset_key, part_per_key)

    use_block = trues(length(offset_key))

    icount = 1

    for i = 2:length(offset_key)-1

        if offset_key[i] == offset_key[icount] + part_per_key[icount]
            part_per_key[icount] += part_per_key[i]
            use_block[i] = false
        else
            icount += 1
        end # if

    end # for

    return use_block, part_per_key
end

"""
    Read particles in box main
"""

"""
    read_particles_in_box(filename::String, blocks::Vector{String},
                          corner_lowerleft, corner_upperright;
                          parttype::Integer=0, verbose::Bool=true)

Reads all particles within a box defined by a lower left and upper right corner
for a given particle type. Returns a dictionary with all requested blocks.
"""
function read_particles_in_box(filename::String, blocks::Vector{String},
                               corner_lowerleft, corner_upperright;
                               parttype::Integer=0, verbose::Bool=true)

    if verbose
        println()
        @info "Running on $(nthreads()) threads"
    end

    filebase = filename
    file_key_index = filebase * ".key.index"

    # if the snapshot does not exist it may be split into multiple files
    if !isfile(filebase)

        if verbose
            @info "File: $filebase not found, looking for sub-files."
        end

        # try reading the first of the distributed snapshots
        filename = filebase * ".0"

        # throw error if file does not exist
        if !isfile(filename)
            error("File: $filename not found!")
        end

        h = head_to_obj(filename)

        if verbose
            @info "$(h.num_files) sub-files found."
        end

        nfiles = h.num_files
    else
        nfiles = 1
    end

    no_mass_block = false

    # check if requested blocks are present
    blocks_in_file = print_blocks(filename, verbose=false)
    for blockname in blocks
        if !block_present(filename, blockname, blocks_in_file)
            if blockname == "MASS"

                no_mass_block = true
                deleteat!(blocks, "MASS")
            else
                error("Block $blockname not present!")
            end
        end
    end

    # read info blocks once here
    snap_info = read_info(filename)
    key_info  = read_info(filename * ".key")

    # prepare dictionary for particle storage, this looks super ugly...
    d = Dict(blocks[i] => Array{snap_info[getfield.(snap_info, :block_name) .== blocks[i]][1].data_type,2}(undef, 0,
            snap_info[getfield.(snap_info, :block_name) .== blocks[i]][1].n_dim) for i = 1:length(blocks))

    # allocate mass array, if it's not in a block
    if no_mass_block
        d["MASS"] = Array{Float32,2}(undef, 0, 1)
    end

    if verbose
        @info "All requested blocks present and dictionary allocated!"
        @info "Checking for .key files..."
    end

    # check if key files are present
    file_key = filename * ".key"
    if !block_present(file_key, "KEY")
        error("No .key file present!")
    end

    if verbose
        @info ".key files found!"
        @info "Calculating peano-hilbert keys..."
        t1 = Dates.now()
    end

    # first read the header
    h_key = read_keyheader(file_key)

    # get a list of the required peano-hilbert keys
    keylist = get_keylist(h_key, corner_lowerleft, corner_upperright)

    if verbose
        t2 = Dates.now()
        @info "$(length(keylist)) Peano-Hilbert keys found. Took: $(t2 - t1)"
        @info "Looking for relevant files..."
        t1 = Dates.now()
    end

    files = find_files_for_keys(filebase, nfiles, keylist)

    if verbose
        t2 = Dates.now()
        @info "$(length(files)) files found. Took: $(t2 - t1)"
    end

    for i = 1:length(files)
        if nfiles > 1
            filename = filebase * ".$(files[i])"
        else
            filename = filebase
        end

        # read header of the file
        h = head_to_obj(filename)

        if h.npart[parttype+1] == 0
            error("No particles of type $parttype in file!")
        end

        filename_keyfile = filename * ".key"

        # read key file data
        h_key = read_keyheader(filename_keyfile)
        keys_in_file = read_block_by_name(filename_keyfile, "KEY",
                                          info = key_info[getfield.(key_info, :block_name) .== "KEY"][1],
                                          parttype = parttype)

        if verbose
            @info "Calculating index list..."
            t1 = Dates.now()
        end

        index_list = get_index_list(keylist, keys_in_file)

        if verbose
            t2 = Dates.now()
            @info "Index list done. Took: $(t2 - t1)"
            @info "Reading $(length(index_list)) key segments..."
        end

        part_per_key = read_block_by_name(filename_keyfile, "NKEY",
                                          info = key_info[getfield.(key_info, :block_name) .== "NKEY"][1],
                                          parttype = parttype)

        n_to_read = sum(part_per_key[index_list])

        if verbose
            @info "Reading $n_to_read particles..."
        end

        offset_key = read_block_by_name(filename_keyfile, "OKEY",
                                          info = key_info[getfield.(key_info, :block_name) .== "OKEY"][1],
                                          parttype = parttype)

        sorted_offset = sortperm(offset_key[index_list])

        # save sorted arrays
        offset_key   = offset_key[index_list[sorted_offset]]
        part_per_key = part_per_key[index_list[sorted_offset]]

        part_per_key_b = part_per_key

        # check if blocks can be joined
        use_block, part_per_key = join_blocks(offset_key, part_per_key)

        if verbose
            @info "Reduced independent blocks from $(length(offset_key)) to $(length(use_block[use_block]))"
        end

        offset_key   = offset_key[use_block]
        part_per_key = part_per_key[use_block]

        pos = get_block_positions(filename)

        # double-check if blocks are present
        blocks_in_file = String.(keys(pos))
        for blockname in blocks
            if !block_present(filename, blockname, blocks_in_file)
                error("Block $blockname not present in file $filename !")
            end
        end

        if verbose
            @info "Reading Blocks..."
            t1 = Dates.now()
        end

        # read blocks in parallel
        @threads for j = 1:length(blocks)

            block_info = snap_info[getfield.(snap_info, :block_name) .== blocks[j]][1]

            # add offset of particle types that should not be read
            offset = 0
            for i=1:length(h.npart)
                if (block_info.is_present[i] > 0) & (h.npart[i] > 0) & ( i < parttype + 1)
                    offset += h.npart[i]
                end
            end
            d[blocks[j]] = read_block_with_offset(filename, d[blocks[j]], pos[blocks[j]],
                                              block_info, offset, offset_key,
                                              n_to_read, part_per_key)

        end # loop over blocks

        if no_mass_block
            d["MASS"] = [d["MASS"]; h.massarr[parttype+1] .* ones(Float32, n_to_read, 1)]
        end

        if verbose
            t2 = Dates.now()
            @info "Blocks read. Took: $(t2 - t1)"
        end

    end # for i = 1:length(files)

    return d
end

"""
    read_particles_in_box(filename::String, blocks::String,
                          corner_lowerleft, corner_upperright;
                          parttype::Integer=0, verbose::Bool=true)

Like `read_particles_in_box` but for a single block. Returns the block as an array.

See also: [`read_particles_in_box`](@ref)
"""
function read_particles_in_box(filename::String, blocks::String,
                               corner_lowerleft, corner_upperright;
                               parttype::Integer=0, verbose::Bool=true)

    d = read_particles_in_box(filename, [blocks], x0, x1, parttype=parttype, verbose=verbose)

    return d[blocks]
end

"""
    read_particles_in_box(filename::String, blocks::Vector{String},
                          center_pos, radius;
                          parttype::Integer=0, verbose::Bool=true)

Reads all particles within a box encapsulating a volume defined by center position
and radius for a given particle type. Returns a dictionary with all requested blocks.

See also: [`read_particles_in_box`](@ref)
"""
function read_particles_in_volume(filename::String, blocks::Vector{String},
                                  center_pos, radius;
                                  parttype::Integer=0, verbose::Bool=true)

    # calculate lower left and upper right corner
    x0 = center_pos .- radius
    x1 = center_pos .+ radius

    return read_particles_in_box(filename, blocks, x0, x1, parttype=parttype, verbose=verbose)
end

"""
    read_particles_in_box(filename::String, blocks::String,
                          center_pos, radius;
                          parttype::Integer=0, verbose::Bool=true)

Like `read_particles_in_volume` but for a single block. Returns the block as an array.

See also: [`read_particles_in_volume`](@ref)
"""
function read_particles_in_volume(filename::String, blocks::String,
                                  center_pos, radius;
                                  parttype::Integer=0, verbose::Bool=true)

    d = read_particles_in_volume(filename, [blocks], center_pos, parttype,
                                 parttype=parttype, verbose=verbose)

    return d[blocks]
end

"""
    Test stuff -- Sandbox
"""

"""

using GadJet
using BenchmarkTools

filebase = "/HydroSims/Magneticum/Box2/hr_bao/snapdir_140/snap_140"
file_key = filebase * ".0.key"

corner_lowerleft  = [-1000.0, -1000.0, -1000.0]
corner_upperright = [ 1000.0,  1000.0,  1000.0]

typeof(corner_lowerleft)

# first read the header
h_key = read_keyheader(file_key)

# get a list of the required peano-hilbert keys
@benchmark keylist = get_keylist(h_key, corner_lowerleft, corner_upperright)
keylist = get_keylist(h_key, corner_lowerleft, corner_upperright)

files_ref = find_files_for_keys(filebase, 512, keylist)
println(files_ref)

@benchmark files = find_files_for_keys_old(filebase, 512, keylist)
files = find_files_for_keys_old(filebase, 512, keylist)
println(files)

filename = filebase * ".324"
filename_keyfile = filename * ".key"

# read key file data
h_key = read_keyheader(filename_keyfile)
key_info = read_info(filename_keyfile)
keys_in_file = read_block_by_name(filename_keyfile, "KEY",
                                  info = key_info[getfield.(key_info, :block_name) .== "KEY"][1],
                                  parttype = 0)

@benchmark get_index_list_serial(keylist, keys_in_file)
@benchmark get_index_list_parallel(keylist, keys_in_file)

idx_serial   = get_index_list_serial(keylist, keys_in_file)
idx_parallel = get_index_list_parallel(keylist, keys_in_file)

sort!(idx_serial)

idx_serial == idx_parallel

@benchmark get_index_list_old(keylist, keys_in_file)

typeof(1) <: Integer

typeof(Array{Float32,2}(undef, 10, 1)) <: AbstractArray

typeof(Array{Integer,2})

typeof(UInt32) <: Int

subtypes(Integer)

supertype(Signed)

typeof(Vector{AbstractFloat}(undef,2)) <: Array{AbstractFloat,1}

typeof(Vector{AbstractFloat}(undef,2))

"""

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


text:
using LBFGSB
using Test

# define a function
function f(x)
    n = length(x)
    y = 0.25 * (x[1] - 1)^2
    for i = 2:n
        y += (x[i] - x[i-1]^2)^2
    end
    4y
end

# and its gradient function that maps a vector x to a vector z
function g!(z, x)
    n = length(x)
    t₁ = x[2] - x[1]^2
    z[1] = 2 * (x[1] - 1) - 1.6e1 * x[1] * t₁
    for i = 2:n-1
        t₂ = t₁
        t₁ = x[i+1] - x[i]^2
        z[i] = 8 * t₂ - 1.6e1 * x[i] * t₁
    end
    z[n] = 8 * t₁
end

@testset "wrapper" begin
    # the first argument is the dimension of the largest problem to be solved
    # the second argument is the maximum number of limited memory corrections
    optimizer = L_BFGS_B(1024, 17)
    n = 25  # the dimension of the problem
    x = fill(Cdouble(3e0), n)  # the initial guess
    # set up bounds
    bounds = zeros(3, n)
    for i = 1:n
        bounds[1,i] = 2  # represents the type of bounds imposed on the variables:
                         #  0->unbounded, 1->only lower bound, 2-> both lower and upper bounds, 3->only upper bound
        bounds[2,i] = isodd(i) ? 1e0 : -1e2  #  the lower bound on x, of length n.
        bounds[3,i] = 1e2  #  the upper bound on x, of length n.
    end
    # call the optimizer
    # - m: the maximum number of variable metric corrections used to define the limited memory matrix
    # - factr: the iteration will stop when (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch,
    #     where epsmch is the machine precision, which is automatically generated by the code. Typical values for factr:
    #     1e12 for low accuracy, 1e7 for moderate accuracy, 1e1 for extremely high accuracy
    # - pgtol: the iteration will stop when max{|proj g_i | i = 1, ..., n} <= pgtol where pg_i is the ith component of the projected gradient
    # - iprint: controls the frequency of output. iprint < 0 means no output:
    #     iprint = 0 print only one line at the last iteration
    #     0 < iprint < 99 print also f and |proj g| every iprint iterations
    #     iprint = 99 print details of every iteration except n-vectors
    #     iprint = 100 print also the changes of active set and final x
    #     iprint > 100 print details of every iteration including x and g
    # - maxfun: the maximum number of function evaluations
    # - maxiter: the maximum number of iterations
    fout, xout = optimizer(f, g!, x, bounds, m=5, factr=1e7, pgtol=1e-5, iprint=-1, maxfun=15000, maxiter=15000)
    @test fout ≈ 1.083490083518441e-9

    #testing the original wrapper:
    function g(x)
        n = length(x)
        z = zeros(n)
        t₁ = x[2] - x[1]^2
        z[1] = 2 * (x[1] - 1) - 1.6e1 * x[1] * t₁
        for i = 2:n-1
            t₂ = t₁
            t₁ = x[i+1] - x[i]^2
            z[i] = 8 * t₂ - 1.6e1 * x[i] * t₁
        end
        z[n] = 8 * t₁
        z
    end
    func(x) = (f(x), g(x))
    fout2, xout2 = optimizer(func, x, bounds, m=5, factr=1e7, pgtol=1e-5, iprint=-1, maxfun=15000, maxiter=15000)
    @test fout2 ≈ 1.083490083518441e-9
end

@testset "lbfgsb" begin
    #  0->unbounded, 1->only lower bound, 2-> both lower and upper bounds, 3->only upper bound
    lb=[-Inf,-Inf,-2.,6f0]
    ub=[Inf,3.,2,Inf]
    btyps=[0,3,2,1]

    @test LBFGSB.typ_bnd.(lb,ub) == btyps

    @test_throws ErrorException LBFGSB.typ_bnd(6,-1.)

    n,m=4,6
    o,b=LBFGSB._opt_bounds(n,m,lb,ub)

    @test size(b)==(3,n)
    @test b[2,:] == lb
    @test b[3,:] == ub
    @test b[1,:] == btyps

    od=L_BFGS_B(n,m)
    @test all(getproperty(od,i)==getproperty(o,i) for i in fieldnames(L_BFGS_B))

    fout,xout=lbfgsb(x->0.5sum(abs,x),(g,x)->g.=x,100*rand(3),lb=[0.5,-Inf,-0.1],ub=[Inf,0.2,Inf])
    @test fout ≈ 0.25
    @test xout ≈ [0.5,0,0]
 end


text:
using FractionalSystems, ControlSystems
using Test

@testset "Test FractionalSystems.jl" begin
    include("fotf.jl")
    include("foss.jl")
    include("FOLE.jl")
end


text:
# Autogenerated wrapper script for ensmallen_jll for any
JLLWrappers.@generate_wrapper_header("ensmallen")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module Advection
  export AdvectionLaw

  import ..Atum
  using ..Atum: constants
  using StaticArrays: SVector

  struct AdvectionLaw{FT, D, C} <: Atum.AbstractBalanceLaw{FT, D, 1, C}
    function AdvectionLaw{FT, D}(u⃗ = ntuple(d->FT(1), D)) where {FT, D}
      new{FT, D, (;u⃗)}()
    end
  end

  Atum.flux(law::AdvectionLaw, q, aux) = SVector(constants(law).u⃗) * q'
  Atum.wavespeed(law::AdvectionLaw, n⃗, q, aux) = abs(n⃗' * SVector(constants(law).u⃗))
  Atum.entropy(law::AdvectionLaw, q, aux) = q' * q / 2
  Atum.entropyvariables(law::AdvectionLaw, q, aux) = q
end