text:
export
    DubinsPathType, SegmentType, DubinsPath,
    LSL, LSR, RSL, RSR, RLR, LRL,
    EDUBOK, EDUBCOCONFIGS, EDUBPARAM, EDUBBADRHO, EDUBNOPATH, EDUBBADINPUT

@enum DubinsPathType LSL LSR RSL RSR RLR LRL
@enum SegmentType L_SEG S_SEG R_SEG

DIRDATA = Dict{Int,Vector{SegmentType}}(
                                        Int(LSL) => [L_SEG, S_SEG, L_SEG],
                                        Int(LSR) => [L_SEG, S_SEG, R_SEG],
                                        Int(RSL) => [R_SEG, S_SEG, L_SEG],
                                        Int(RSR) => [R_SEG, S_SEG, R_SEG],
                                        Int(RLR) => [R_SEG, L_SEG, R_SEG],
                                        Int(LRL) => [L_SEG, R_SEG, L_SEG]
                                       )

"""
The data structure that holds the full dubins path.

Its data fields are as follows:

* the initial configuration, qi,
* the params vector that contains the length of each segment, params,
* the turn-radius, ρ, and,
* the Dubins path type given by the @enum DubinsPathType
"""
mutable struct DubinsPath
    qi::Vector{Float64}            # the initial configuration
    params::Vector{Float64}        # the lengths of the three segments
    ρ::Float64                     # turn radius
    path_type::DubinsPathType   # the path type
end

"""
Empty constructor for the DubinsPath type
"""
DubinsPath() = DubinsPath(zeros(3), zeros(3), 1., LSL)

"""
This data structure holds the information to compute the Dubins path
in the transformed coordinates where the initial (x,y) is translated to the
origin, the final the coordinate axis is rotated to make the x-axis aligned with
the line joining the two points. The variable names follow the convention used
in the paper "Classification of the Dubins set" by Andrei M. Shkel and Vladimir Lumelsky
"""
mutable struct DubinsIntermediateResults
    α::Float64                  # transformed α
    β::Float64                  # transformed β
    d::Float64                  # transformed d
    sa::Float64                 # sin(α)
    sb::Float64                 # sin(β)
    ca::Float64                 # cos(α)
    cb::Float64                 # cos(β)
    c_ab::Float64               # cos(α-β)
    d_sq::Float64               # d²
end

"""
Empty constructor for the DubinsIntermediateResults data type
"""
function DubinsIntermediateResults(q0::Vector{Float64}, q1::Vector{Float64}, ρ::Float64)

    ir = DubinsIntermediateResults(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.)

    dx = q1[1] - q0[1]
    dy = q1[2] - q0[2]
    D = sqrt(dx*dx + dy*dy)
    d = D / ρ
    Θ = 0

    # test required to prevent domain errors if dx=0 and dy=0
    (d > 0) && (Θ = mod2pi(atan(dy, dx)))
    α = mod2pi(q0[3] - Θ)
    β = mod2pi(q1[3] - Θ)
    ir.α = α
    ir.β = β
    ir.d = d
    ir.sa = sin(α)
    ir.sb = sin(β)
    ir.ca = cos(α)
    ir.cb = cos(β)
    ir.c_ab = cos(α-β)
    ir.d_sq = d*d

    return ir
end

const EDUBOK = 0                # no error
const EDUBCOCONFIGS = 1         # colocated configurations
const EDUBPARAM = 2             # path parameterization error
const EDUBBADRHO = 3            # the rho value is invalid
const EDUBNOPATH = 4            # no connection between configurations with this word
const EDUBBADINPUT = 5          # uninitialized inputs to functions
TOL = 1e-10                     # tolerance


text:
include("../utils.jl")
include("group_utils.jl")

@testset "Semidirect product group" begin
    M1 = TranslationGroup(2)
    A = TranslationAction(M1, M1)
    G = SemidirectProductGroup(M1, M1, A)
    @test G === GroupManifold(
        TranslationGroup(2) × TranslationGroup(2),
        Manifolds.SemidirectProductOperation(A),
    )
    @test repr(G) == "SemidirectProductGroup($(M1), $(M1), $(A))"
    @test repr(G.op) == "SemidirectProductOperation($(A))"
    M = base_manifold(G)
    @test M === TranslationGroup(2) × TranslationGroup(2)

    ts1 = Vector{Float64}.([1:2, 2:3, 3:4])
    ts2 = Vector{Float64}.([1:2, 2:3, 3:4]) .* 10
    tuple_pts = [zip(ts1, ts2)...]
    tuple_v = ([-1.0, 1.0], [2.0, 3.0])

    shape_se = Manifolds.ShapeSpecification(Manifolds.ArrayReshaper(), M.manifolds...)
    pts = [Manifolds.prod_point(shape_se, tp...) for tp in tuple_pts]
    v_pts = [Manifolds.prod_point(shape_se, tuple_v...)]

    X = log(G, pts[1], pts[1])
    Y = zero_tangent_vector(G, pts[1])
    Z = Manifolds.allocate_result(G, zero_tangent_vector, pts[1])
    Z = zero_tangent_vector!(M, Z, pts[1])
    @test norm(G, pts[1], X) ≈ 0
    @test norm(G, pts[1], Y) ≈ 0
    @test norm(G, pts[1], Z) ≈ 0

    e = Identity(G, (zeros(2), zeros(2)))
    @test inv(G, e) === e
    @test identity(G, e) === e
    @test identity!(G, e, e) === e

    @test compose(G, e, pts[1]) == pts[1]
    @test compose(G, pts[1], e) == pts[1]
    @test compose(G, e, e) === e
end

ext: jl


text:
test_str = "35
20
15
25
47
40
62
55
65
95
102
117
150
182
127
219
299
277
309
576"

test_list = split(test_str, '\n')

list = readlines("day9_input.txt")

function check(sublist)
    sublen = length(sublist)
    last = sublist[end]
    for j in 1:sublen-2
        for k in j+1:sublen-1
            if sublist[j] + sublist[k] == last
                return true
            end
        end
    end
    return false
end

function invalid(preamble, list)
    int_list = map(x-> parse(Int, x), list)
    for i in preamble+1:length(int_list)
        if !check(int_list[i-preamble:i])
            return i, int_list[i]
        end
    end
    return nothing, nothing
end

@show invalid(5, test_list)
@show invalid(25, list)

function contiguous(preamble, list)
    i_inv, num_inv = invalid(preamble, list)
    int_list = map(x-> parse(Int, x), list)
    len = length(list)
    for i in 1:len
        j = i
        sublistsum = int_list[i]
        while sublistsum <= num_inv
            if sublistsum == num_inv
                return minimum(int_list[i:j]) + maximum(int_list[i:j])
            end
            j += 1
            sublistsum += int_list[j]
        end
    end
end

@show contiguous(5, test_list)
@show contiguous(25, list)

nothing
#=
```julia
julia> @btime invalid(25, list);
  166.235 μs (526 allocations: 155.34 KiB)

julia> @btime contiguous(25, list);
  503.429 μs (531 allocations: 163.75 KiB)
```
 =#
# I restart counting with the second counter j every time I step up the
# first counter i, but that is indeed a waste. I could just subtract the
# first value in the sublist


text:
using jInv.LinearSolvers
using Test
using ConvDiffMIPDECO
using jInv.ForwardShare
using jInv.Mesh
using LinearAlgebra
using SparseArrays
using MUMPSjInv

domain = [0 1. 0 1]

N = 2*[16;32;64;128]
v = randn(2)
v = v/norm(v)
V = X -> [v[1]*ones(size(X,1)) v[2]*ones(size(X,1))]
bc = (:neu,:dir,:neu,:dir)

u  = X -> cos.(pi*X[:,1]).*cos.(pi*X[:,2])

ux = X -> -pi*sin.(pi*X[:,1]).*cos.(pi*X[:,2])
uy = X -> -pi*cos.(pi*X[:,1]).*sin.(pi*X[:,2])
sig = .001;

f  = X->  +2*sig*pi^2*u(X)+v[1]*ux(X)+v[2]*uy(X)

err = zeros(length(N))

bicg = (A,b; M=identity,tol=1e-10,maxIter=500,out=1)-> bicgstb(A,b,M1=identity,tol=tol,maxIter=maxIter,out=out)
times = zeros(length(N),3)
for j=1:2
	for k=1:length(N)
		Ainvs = (getJuliaSolver(), getMUMPSsolver(), ConvDiffMIPDECO.getBICGSTB(out=1))

	nk      = [N[k];N[k]]
	Mk  = getRegularMesh(domain,nk)
	xck = getCellCenteredGrid(Mk)

	pFork   = getConvDiffParam(Mk,V,sig=sig,gd=u,gn=ux,bc=bc,Ainv=Ainvs[j])

	dobsk,pFork = getData(f(xck), pFork)

	utruek = u(xck)
	times[k,j]= @elapsed begin
		err[k] = norm(pFork.Fields - utruek,Inf)/norm(utruek,Inf)
	end
	println("err[$k] = $(err[k])")
end
@test all(abs.(diff(log2.(err))).>0.8)
end

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
export f_∇f

function f_∇f(objective::ObjectiveType,sim::Simulation ; solver  = LU() , args...)
	linsys = LinearSystem(sim)
	A ,b = linsys.A , linsys.b

    x , x_adj = linsolve(A,b,adjoint_src(objective), solver ; args...)

    ∇ = normalize(2 * (2pi/sim.λ₀)^2 * real(x .* x_adj),Inf)

    ∇ϵx = sim(extractreshape,∇,x̂)
    ∇ϵy = sim(extractreshape,∇,ŷ)
    ∇ϵz = sim(extractreshape,∇,ẑ)

    sim.E = x
	return objective_value(objective,x) , Dict([(x̂, ∇ϵx), (ŷ, ∇ϵy),(ẑ,∇ϵz)])

end


text:
@compat abstract type CubicSpline{T <: Real} end

immutable ImmutableCubicSpline{T<:Real} <: CubicSpline{T}
    x::Vector{T}
    fx::Vector{T}
    fdp::Vector{T}
    npts::Ref{Cint}
end

type MutableCubicSpline{T<:Real} <: CubicSpline{T}
    x::Vector{T}
    fx::Vector{T}
    fdp::Vector{T}
    npts::Ref{Cint}
    dirty::Bool
end

const spline_lib = Pkg.dir()*"\CubicSplineNR\lib\spline1"

"""
    _dodp!(...)

    calls Fortran function from Numerical Recipes, provided by Tony
    Smith

    It computes the second derivatives needed for calculating the
    cubic spline interpolation.
"""
function _dodp!(x, fx, fdp, npts)
    ccall((:__procedures_MOD_dodp, spline_lib),
       Void,
       (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}),
       x, fx, fdp, npts)
end

function CubicSpline(x::Vector{Float64}, fx::Vector{Float64}; mutable::Bool=false)
    N = length(x)
    npts = Ref{Cint}(N)
    fdp = zeros(N)

    _dodp!(x, fx, fdp, npts)

    if mutable
        return MutableCubicSpline(copy(x), copy(fx), fdp, npts, false)
    else
        return ImmutableCubicSpline(x, fx, fdp, npts)
    end
end

function update!(cs::MutableCubicSpline, fx::Vector{Float64})
    cs.fdp .= 0.0
    cs.fx .= fx .+ 0.0

    dodp!(cs.x, cs.fx, cs.fdp, cs.npts)

    cs.dirty = false
end

function update!(cs::MutableCubicSpline, x::Vector{Float64}, fx::Vector{Float64})
    cs.fdp .= 0.0
    cs.x .= x .+ 0.0
    cs.fx .= fx .+ 0.0

    dodp!(cs.x, cs.fx, cs.fdp, cs.npts)

    cs.dirty = false
end

# function set!(cs::MutableCubicSpline, i::Int, value::Float64)
#     cs.fx[i] = value
#     cs.dirty = true
# end

function _interp!(point::Float64,
                 x::Vector, fx::Vector, fdp::Vector, npts,
                 level::Int, deriv::Int, second::Int,
                 y::Vector{Float64}, yp::Vector{Float64}, ydp::Vector{Float64})
    ny = Ref{Cint}(level)
    nyp = Ref{Cint}(deriv)
    nydp = Ref{Cint}(second)

    point = Ref{Cdouble}(point)

    ccall((:__procedures_MOD_interp, "spline1"), Void,
(Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
    point, x, fx, fdp, y, yp, ydp, ny, nyp, nydp, npts)
end

function interp_level{CS <: CubicSpline{Float64}}(cs::CS, x::Float64)
    point = Ref{Cdouble}(x)

    y = zeros(1)
    yp = zeros(1)
    ydp = zeros(1)

    _interp!(x, cs.x, cs.fx, cs.fdp, cs.npts, 1, 0, 0, y, yp, ydp)
#     ccall((:__procedures_MOD_interp, "spline1"), Void,
# (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
#     point, cs.x, cs.fx, cs.fdp, y, yp, ydp, ny, nyp, nydp, cs.npts)
     y[1]
end

#function interp_level{T<:Float64, CS <: CubicSpline{T}}(cs::CS, x::ForwardDiff.Dual{T})
function interp_level{T<:Float64, CS <: CubicSpline{Float64}}(cs::CS, x::ForwardDiff.Dual{T})

    ny = Ref{Cint}(1)
    nyp = Ref{Cint}(1)
    nydp = Ref{Cint}(0)

    point = Ref{Cdouble}(ForwardDiff.value(x))

    y = zeros(1)
    yp = zeros(1)
    ydp = zeros(1)

    ccall((:__procedures_MOD_interp, "spline1"), Void,
(Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
    point, cs.x, cs.fx, cs.fdp, y, yp, ydp, ny, nyp, nydp, cs.npts)
    ForwardDiff.Dual(y[1]::T, yp[1]::T * ForwardDiff.partials(x))
end

function(itp::ImmutableCubicSpline)(x::Real) #where T <: Real
    interp_level(itp, x)
end

function(itp::MutableCubicSpline)(x::Real; update::Bool=true) #where T <: Real
    if update && itp.dirty
        update!(itp, itp.fx)
    end
    interp_level(itp, x)
end


text:
"""
    QemuRunner

A `QemuRunner` represents an "execution context", an object that bundles all
necessary information to run commands within the container that contains
our crossbuild environment.  Use `run()` to actually run commands within the
`QemuRunner`, and `runshell()` as a quick way to get an interactive shell
within the crossbuild environment.
"""
mutable struct QemuRunner <: Runner
    qemu_cmd::Cmd
    append_line::String
    sandbox_cmd::Cmd
    env::Dict{String, String}
    platform::Platform
end

qemu_url = "https://github.com/Keno/QemuBuilder/releases/download/hvf/qemu.x86_64-apple-darwin14.tar.gz"
qemu_hash = "92fa439350970f673a7324b61ed5bc4894d6665543175f01135550cb4a458cde"
kernel_url = "https://github.com/Keno/LinuxBuilder/releases/download/v4.15-rc2/linux.x86_64-linux-gnu.tar.gz"
kernel_hash = "0d02413529e635d4af6d2122f6aba22d261f43616bb286da3855325988f9ac3b"

"""
    update_qemu(;verbose::Bool = false)

Update our QEMU and Linux kernel installations, downloading and installing them
into the `qemu_cache` directory that defaults to `deps/qemu`.
"""
function update_qemu(;verbose::Bool = false)
    download_verify_unpack(
        qemu_url,
        qemu_hash,
        qemu_cache;
        tarball_path=downloads_dir(basename(qemu_url)),
        verbose=verbose,
        force=true
    )
    download_verify(
        kernel_url,
        kernel_hash,
        downloads_dir(basename(kernel_url));
        verbose=verbose,
        force=true
    )
    if !isfile(kernel_path())
        unpack(downloads_dir(basename(kernel_url)), qemu_cache; verbose=verbose)
    end
end

qemu_path() = joinpath(qemu_cache, "usr/local/bin/qemu-system-x86_64")
kernel_path() = joinpath(qemu_cache, "bzImage")

function platform_def_mappings(platform)
    tp = triplet(platform)
    mapping = Pair{String,String}[
        shard_path_squashfs(tp) => joinpath("/opt", tp)
    ]

    # We might also need the x86_64-linux-gnu platform for bootstrapping,
    # so make sure that's always included
    if platform != Linux(:x86_64)
        ltp = triplet(Linux(:x86_64))
        push!(mapping, shard_path_squashfs(ltp) => joinpath("/opt", ltp))
    end

    # If we're trying to run macOS and we have an SDK directory, mount that!
    if platform == MacOS()
        sdk_version = "MacOSX10.10.sdk"
        sdk_shard_path = joinpath(shards_cache, sdk_version)
        push!(mapping, sdk_shard_path => joinpath("/opt", tp, sdk_version))
    end

    # Reverse mapping order, because `sandbox` reads them backwards
    reverse!(mapping)
    return mapping
end

platform_accelerator() = Compat.Sys.islinux() ? "kvm" : "hvf"

function QemuRunner(workspace_root::String; cwd = nothing,
                      platform::Platform = platform_key(),
                      extra_env=Dict{String, String}(),
                      verbose::Bool = false,
                      mappings::Vector = platform_def_mappings(platform),
                      workspaces::Vector = [])
    global use_ccache

    # Ensure the rootfs for this platform is downloaded and up to date
    update_rootfs(
        platform != Linux(:x86_64) ?
        [triplet(platform), triplet(Linux(:x86_64))] :
        triplet(platform); verbose=verbose, squashfs=true, mount=false)
    update_qemu(;verbose=verbose)

    qemu_cmd = ```
        $(qemu_path()) -kernel $(kernel_path())
        -M accel=$(platform_accelerator()) -nographic
        -drive if=virtio,file=$(rootfs_path_squashfs()),format=raw
        -nodefaults -rtc base=utc,driftfix=slew
        -device virtio-serial -chardev stdio,id=charconsole0
        -device virtconsole,chardev=charconsole0,id=console0
        -device virtserialport,chardev=charcomm0,id=comm0
        -device virtio-net-pci,netdev=networking -netdev user,id=networking
    ```
    append_line = "quiet console=hvc0 root=/dev/vda rootflags=ro rootfstype=squashfs noinitrd init=/sandbox"

    sandbox_cmd = ``
    if verbose
        # We put --verbose at the beginning here so that we can see argument parsing
        sandbox_cmd = `$sandbox_cmd --verbose`
    end

    if cwd != nothing
        sandbox_cmd = `$sandbox_cmd --cd $cwd`
    end

    # We always include the workspace as one of our read-write workspace mountings, and always as the first
    insert!(workspaces, 1, workspace_root => "/workspace")

    # If we're enabling ccache, then mount in a read-writeable volume at /root/.ccache
    if use_ccache
        if !isdir(ccache_dir())
            mkpath(ccache_dir())
        end
        push!(workspaces, ccache_dir() => "/root/.ccache")
    end

    devchr = 'b'
    devnum = 0
    # Mount in read-writable workspaces over plan9
    for (outside, inside) in workspaces
        qemu_cmd = `$qemu_cmd -fsdev local,security_model=none,id=fsdev$(devnum),path=$(outside) -device virtio-9p-pci,id=fs$(devnum),fsdev=fsdev$(devnum),mount_tag=workspace$(devnum)`
        sandbox_cmd = `$sandbox_cmd --workspace 9p/workspace$(devnum):$(inside)`
        devnum += 1
    end

    # Mount in read-only mappings over plan9 if they're raw directories, and directly if they're .squashfs files
    for (outside, inside) in mappings
        # Squashfs files get mounted in directly, actual directories have to be shared over plan9:
        if endswith(outside, ".squashfs")
            qemu_cmd = `$qemu_cmd -drive if=virtio,file=$(outside),format=raw`
            sandbox_cmd = `$sandbox_cmd --map /dev/vd$(devchr):$(inside)`
            devchr = Char(Int(devchr) + 1)
        else
            qemu_cmd = `$qemu_cmd -fsdev local,security_model=none,id=fsdev$(devnum),path=$(outside),readonly -device virtio-9p-pci,id=fs$(devnum),fsdev=fsdev$(devnum),mount_tag=mapping$(devnum)`
            sandbox_cmd = `$sandbox_cmd --map 9p/mapping$(devnum):$(inside)`
            devnum += 1
        end
    end

    QemuRunner(qemu_cmd, append_line, sandbox_cmd, merge(target_envs(triplet(platform)), extra_env), platform)
end

function show(io::IO, x::QemuRunner)
    p = x.platform
    # Displays as, e.g., Linux x86_64 (glibc) QemuRunner
    write(io, "$(typeof(p).name.name)", " ", arch(p), " ",
          Compat.Sys.islinux(p) ? "($(p.libc)) " : "",
          "QemuRunner")
end

function qemu_gen_cmd(qr::QemuRunner, cmd::Cmd, comm_socket_path::String)
    @async begin
        # This is what we'll send to `sandbox`
        sandbox_args = String.(`$(qr.sandbox_cmd) -- $(cmd)`.exec)
        sandbox_env = ["$k=$v" for (k, v) in qr.env]

        # Wait until QEMU starts up the communications channel
        start_time = time()
        while !issocket(comm_socket_path)
            sleep(0.01)
            if time() - start_time > 2
                Compat.@info("Unable to establish communication with QEMU, does $(comm_socket_path) exist?")
            end
        end

        # Once it has, open it up,
        commsock = connect(comm_socket_path)

        # We're going to spit out:
        #  [4 bytes] number of sandbox args
        #    for each sandbox arg:
        #      [4 bytes] length of arg
        #      [N bytes] arg literal string
        #  [4 bytes] number of environment mappings
        #    for each environment mapping
        #      [4 bytes] length of mapping
        #      [N bytes] mapping literal string
        write(commsock, UInt32(length(sandbox_args)))
        for idx in 1:length(sandbox_args)
            write(commsock, UInt32(length(sandbox_args[idx])))
            write(commsock, sandbox_args[idx])
        end

        write(commsock, UInt32(length(sandbox_env)))
        for idx in 1:length(sandbox_env)
            write(commsock, UInt32(length(sandbox_env[idx])))
            write(commsock, sandbox_env[idx])
        end

        # Wait for acknowledgement
        read(commsock, UInt8)

        close(commsock)
    end

    long_cmd = ```
        $(qr.qemu_cmd)
        -device virtio-serial -chardev socket,path=$(comm_socket_path),server,nowait,id=charcomm0
        -append $(qr.append_line)
    ```

    return `$long_cmd`
end

function Base.run(qr::QemuRunner, cmd, logpath::AbstractString; verbose::Bool = false, tee_stream=stdout)
    return temp_prefix() do prefix
        comm_socket_path = joinpath(prefix.path, "qemu_comm.socket")
        # Launch QEMU
        cmd = qemu_gen_cmd(qr, cmd, comm_socket_path)

        oc = OutputCollector(cmd; verbose=verbose, tee_stream=tee_stream)
        did_succeed = wait(oc)

        if !isempty(logpath)
            # Write out the logfile, regardless of whether it was successful
            mkpath(dirname(logpath))
            open(logpath, "w") do f
                # First write out the actual command, then the command output
                println(f, cmd)
                print(f, merge(oc))
            end
        end

        # Return whether we succeeded or not
        return did_succeed
    end
end

function run_interactive(qr::QemuRunner, cmd::Cmd; stdin = nothing, stdout = nothing, stderr = nothing)
    temp_prefix() do prefix
        comm_socket_path = joinpath(prefix.path, "qemu_comm.socket")
        cmd = qemu_gen_cmd(qr, cmd, comm_socket_path)
        if stdin isa Base.AbstractCmd || stdin isa Base.TTY
            cmd = pipeline(cmd, stdin=stdin)
        end
        if stdout isa Base.AbstractCmd || stdout isa Base.TTY
            cmd = pipeline(cmd, stdout=stdout)
        end
        if stderr isa Base.AbstractCmd || stderr isa Base.TTY
            cmd = pipeline(cmd, stderr=stderr)
        end

        if stdout isa IOBuffer
            if !(stdin isa IOBuffer)
                stdin = devnull
            end
            out, process = open(cmd, "r", stdin)
            @schedule begin
                while !eof(out)
                    write(stdout, read(out))
                end
            end
            wait(process)
        else
            run(cmd)
        end
    end

    # Qemu appears to mess with our terminal
    Base.reseteof(Compat.stdin)
end

function runshell(qr::QemuRunner, args...; kwargs...)
    println("Setting parent shell interrupt to ^] (use ^C as usual, ^] to interrupt Julia)")
    @static if Compat.Sys.isapple() run(`stty intr ^\]`) end
    run_interactive(qr, `/bin/bash`, args...; kwargs...)
    @static if Compat.Sys.isapple() run(`stty intr ^\C`) end
end

function runshell(::Type{QemuRunner}, platform::Platform = platform_key(); verbose::Bool = false)
    qr = QemuRunner(
        pwd();
        cwd="/workspace/",
        platform=platform,
        verbose=verbose
    )
    return runshell(qr)
end


text:
using Random
using MLJBase
using Test

mutable struct Foo <: MLJType
    rng::AbstractRNG
    x::Int
end

mutable struct Bar <: MLJType
    rng::AbstractRNG
    x::Int
end

f1 = Foo(MersenneTwister(7), 1)
f2 = Foo(MersenneTwister(8), 1)
@test f1.rng != f2.rng
@test f1 == f2
f1.x = 2
@test f1 != f2
b = Bar(MersenneTwister(7), 1)
@test f1 != b

true


text:
using NCon
using Test
using TensorOperations

# Some random tensors of various sizes to be used in the tests.
A = rand(3,2,5)
B = rand(5,6)
C = rand(2,2,5,4,5)
D1, D2, D3 = rand(10,8), rand(10,8), rand(10,8)
E = rand(3,2,3,2)
I = rand(-100:100, 3,2)
J = rand(-100:100, 5,2,3,5)

# Test various in valid calls. check_indices should catch all of these.
# More index lists than tensors.
@test_throws(ArgumentError,
             ncon((A, B), ([-1,-2,1], [1,-3], [-4]); check_indices=true))
# Fewer index lists than tensors.
@test_throws(ArgumentError,
             ncon((A, B), ([-1,-2,1]); check_indices=true))
# Contraction index appearing thrice.
@test_throws(ArgumentError,
             ncon((A, B), ([-1,1,1], [1,-3]); check_indices=true))
# Free index appearing twice.
@test_throws(ArgumentError,
             ncon((A, B), ([-1,-1,1], [1,-3]); check_indices=true))
# Zero in index lists.
@test_throws(ArgumentError,
             ncon((A, B), ([-1,-2,0], [0,-3]); check_indices=true))
# Incompatible contraction indices (different dimension).
@test_throws(ArgumentError,
             ncon((A, B), ([-1,-2,1], [-3,1]); check_indices=true))
# Negative number in order.
@test_throws(ArgumentError,
             ncon((A, E), ([2,1,-3], [-1,-2,2,1]); order=[1,2,-5],
                  check_indices=true))
# Not all contraction indices in order.
@test_throws(ArgumentError,
             ncon((A, E), ([2,1,-3], [-1,-2,2,1]); order=[1],
                  check_indices=true))
# An extraneous positive number in order.
@test_throws(ArgumentError,
             ncon((A, E), ([2,1,-3], [-1,-2,2,1]); order=[1,2,3],
                  check_indices=true))
# Positive number in forder.
@test_throws(ArgumentError,
             ncon((A, E), ([2,1,-3], [-1,-2,2,1]); forder=[-1,-2,1],
                  check_indices=true))
# Not all free indices in forder.
@test_throws(ArgumentError,
             ncon((A, E), ([2,1,-3], [-1,-2,2,1]); forder=[-1],
                  check_indices=true))
# An extraneous negative number in forder.
@test_throws(ArgumentError,
             ncon((A, E), ([2,1,-3], [-1,-2,2,1]); forder=[-1,-2,-5],
                  check_indices=true))

# Do various contractions using ncon and the @tensor macro, as well as a pure
# index permutation and a noop, and check that the results match.
# Note that the different calls use tuples/Arrays for the arguments in
# different combinations. They should all be equally valid.
con = ncon((A, B), ([-1,-2,1], [1,-3]); check_indices=true)
@tensor reco[a,b,c] := A[a,b,i] * B[i,c]
@test isapprox(con, reco)

con = ncon((A, B, C), ((-1,1,2), (2,-2), (-4,1,3,-3,3)); check_indices=true)
@tensor reco[a,b,c,d] := A[a,k,j] * B[j,b] * C[d,k,i,c,i]
@test isapprox(con, reco)

con = ncon((A, C), ((-1,2,1), (-2,2,-3,-4,1)); check_indices=true)
@tensor reco[a,b,c,d] := A[a,j,i] * C[b,j,c,d,i]
@test isapprox(con, reco)

@test permutedims(C, [3,1,2,4,5]) == ncon(C, [-11,-22,-33,-44,-55];
                                          forder=[-33,-11,-22,-44,-55],
                                          check_indices=true)

@test D1 == ncon(D1, [-1,-2]; check_indices=true)

con = ncon((D1, D2), ([-1,-2], [-3,-4]); check_indices=true)
@tensor reco[a,b,c,d] := D1[a,b] * D2[c,d]
@test isapprox(con, reco)

con = ncon((D1, D2, D3), ([-1,-2], [-6,-5], [-3,-4]);
           forder=[-2,-4,-5,-1,-6,-3], check_indices=true)
@tensor reco[-2,-4,-5,-1,-6,-3] := D1[-1,-2] * D2[-6,-5] * D3[-3,-4]
@test isapprox(con, reco)

con = ncon(E, [1,2,1,2]; check_indices=true)
@tensor reco[] := E[i,j,i,j]
@test isapprox(con, reco)

con = ncon((I, J), ([1,2], [-2,2,1,-1]); check_indices=true)
reco = tensorcontract(I, [:i,:j], J, [:b,:j,:i,:a], [:a,:b])
@test isapprox(con, reco)

con = ncon(J, [1,-1,-2,1]; check_indices=true)
@tensor reco[a,b] := J[i,a,b,i]
@test isapprox(con, reco)

# test whether it works for non-blastypes
A = Float16.(A)
B = Float16.(B)
C = Float16.(C)
con = ncon((A, B, C), ((-1,1,2), (2,-2), (-4,1,3,-3,3)); check_indices=true)
@tensor reco[a,b,c,d] := A[a,k,j] * B[j,b] * C[d,k,i,c,i]
@test isapprox(con, reco)


text:
#=
Here I show the exact correspondence
my spiking implementation and code by
Litwin-Kumar / Auguste Schulz

Here I just compare weight matrix,definition of assemblies, input function

=#

# include inbuilt modules
using Distributions
using SparseArrays
using Dates
using LinearAlgebra
using Random
using Plots ; theme(:dark)
using Random ; Random.seed!(0)
using ProgressMeter
using SpikingRNNs; const global S = SpikingRNNs
using SparseArrays

# --------------------- include functions  ------------------------------------------------
# include runs the respective julia code, i.e. defined functions are then in the workspace

#function gennextsequencepretraining(stim::Array{Float64,2}, firstimg, novelimg; Nimg = 4, Ntrain = 20, Nreps = 20, Nseq = 5, Nblocks = 1, stimstart = 1000, lenstim = 1000, lenpause = 3000, strength = 8)
function genstimparadigmpretraining(stim::Array{Float64,2}; Nimg = 4, Nass = 20, Ntrain = 20, stimstart = 1000, lenstim = 1000, lenpause = 3000, strength = 8)
	# generate a stimulus sequence
	# where each assembly is
	# Nimg Number of images per sequence
	# Nreps Number of seqeunce repititions in one block
	# Nseq Number of sequences
	# Nblocks Number of block repititions

	if Ntrain != 0
		Nstim = Nass * Ntrain
		tempstim = zeros(Nstim, 4)

		#tempstim[1,1] = firstimg
		# repeat sequence several times, i.e. assemblie numbers 123412341234...
		#tempstim[:,1] = repeat(firstimg:(firstimg+Nimg-1), outer = Nreps)
		# ensure that images are shuffled but never the same image after the other
		images = 1:Nass
		assemblysequence = shuffle(copy(images))
		for rep = 2:Ntrain
			shuffleimg = shuffle(images)
			#println("begin $(shuffleimg[1])")
			while shuffleimg[1] == assemblysequence[end]
				println("shuffle again")
				shuffleimg = shuffle(images)
				#println("new begin $(shuffleimg[1])")
			end
		assemblysequence = vcat(assemblysequence, shuffleimg)
		end
		tempstim[:,1] .= assemblysequence
		if stim[end,3] == 0
			tempstim[1,2] = stimstart
		else
			tempstim[1,2] = stim[end,3]+lenpause
		end
		tempstim[1,3] = tempstim[1,2] + lenstim
		tempstim[1,4] = strength

		for i = 2:Nstim#Nimg*Nreps #number in image times
			# # lenstim ms stimulation and lenpause ms wait time
			tempstim[i,2] = tempstim[i-1,2] + lenstim + lenpause
			tempstim[i,3] = tempstim[i,2]+lenstim
			tempstim[i,4] = strength
		end
		if stim[end,3] == 0
			stim = tempstim
		else
			stim = vcat(stim, tempstim)
		end
		#tempstim = nothing
	end
	return stim
end # function genstim

function gennextsequenceNonovel(stim::Array{Float64,2}, firstimg, novelimg; Nimg = 4, Nreps = 20, Nseq = 5, Nblocks = 1, stimstart = 1000, lenstim = 1000, lenpause = 3000, strength = 8)
	# generate a stimulus sequence
	# where each assembly is
	# Nimg Number of images per sequence
	# Nreps Number of seqeunce repititions in one block
	# Nseq Number of sequences
	# Nblocks Number of block repititions
	Nstim = Nimg*Nreps
	tempstim = zeros(Nstim, 4)
    #tempstim[1,1] = firstimg
	tempstim[:,1] = repeat(firstimg:(firstimg+Nimg-1), outer = Nreps)

	if stim[end,3] == 0
		tempstim[1,2] = stimstart
	else
    	tempstim[1,2] = stim[end,3]+lenpause
	end
    tempstim[1,3] = tempstim[1,2] + lenstim
    tempstim[1,4] = strength

    for i = 2:Nimg*Nreps #number in image times
        # lenstim ms stimulation and lenpause ms wait time
        tempstim[i,2] = tempstim[i-1,2] + lenstim + lenpause
        tempstim[i,3] = tempstim[i,2]+lenstim
        tempstim[i,4] = strength
    end
	#tempstim[end-Nimg,1] = novelimg
	if stim[end,3] == 0
		stim = tempstim
	else
		stim = vcat(stim, tempstim)
	end
	#tempstim = nothing
    return stim
end # function genstim

function genstimparadigmNonovel(stimulus; Nimg = 4, Nreps = 20, Nseq = 5, Nblocks = 1, stimstart = 1000, lenstim = 1000, lenpause = 3000, strength = 8)

	firstimages = collect(1:Nimg:Nimg*Nseq)
	novelimg = firstimages .+ Nimg*Nseq
	novelrep = collect(Nimg*Nseq:Nimg*Nseq+Nimg).+1

	# for img in firstimages
	# 	stimulus = gennextsequence!(stimulus, img)
	# end
	blockonset = []
	storefirstimages = copy(firstimages)
	for b = 1:Nblocks
		if b == 1
			for img in firstimages
				append!(blockonset, size(stimulus,1) + 1)
				stimulus = gennextsequenceNonovel(stimulus, img, img + Nseq * Nimg, Nimg = Nimg, Nreps = Nreps, Nseq = Nseq, Nblocks = Nblocks, stimstart = stimstart, lenstim = lenstim, lenpause = lenpause, strength = strength)
			end
		else
			shuffleimg = shuffle(firstimages)
			while shuffleimg[1] == storefirstimages[end]
				shuffleimg = shuffle(firstimages)
			end
			storefirstimages = vcat(storefirstimages, shuffleimg)
			for img in shuffleimg # since we shuffle ensure that it is still the correct novel image
				append!(blockonset, size(stimulus,1) + 1)
				stimulus = gennextsequenceNonovel(stimulus, img, img + Nseq * Nimg,Nimg = Nimg, Nreps = Nreps, Nseq = Nseq, Nblocks = Nblocks, stimstart = stimstart, lenstim = lenstim, lenpause = lenpause, strength = strength)
			end
		end
	end
	# ensure mapping from novel img to right value despite shuffling
	for i = 1:length(novelimg)
		stimulus[stimulus[:,1] .== novelimg[i],1] .= novelrep[i]
	end

	# ensure that blockonset 1 when we start with an emptz array
	if blockonset[1] == 2
		blockonset[1] = 1
	end
	return stimulus, convert(Array{Int64,1}, blockonset)

end

function genstimparadigmnovelcont(stimulus; Nimg = 4, Nreps = 20, Nseq = 5, Nblocks = 1, stimstart = 1000, lenstim = 1000, lenpause = 3000, strength = 8)
  """	generate the stimulation paradigm ensure that for each block there is a different novel image and
  that two consecutive blocks are never the same

  generate seq. violation stimulation paradigm
  ABCABCABC....ABN1ABC
  DEFDEFDEF....DEN2DEF
  XYZXYZXYZ...XYN3XYZ
  DEFDEFDEF....DEN4DEF
  ABCABCABC....ABN1ABC
  """

    firstimages = collect(1:Nimg:Nimg*Nseq)
    novelimg = Nimg*Nseq + 1 # start counting from the first assembly after the core assemblies
    novelrep = collect(Nimg*Nseq:Nimg*Nseq+Nimg).+1

    storefirstimages = copy(firstimages)
    blockonset = []

    for b = 1:Nblocks
      #println("block $(b) -------------------")
      if b == 1
        for img in firstimages
          append!(blockonset, size(stimulus,1) + 1)
          stimulus = gennextsequencenovel(stimulus, img, novelimg, Nimg = Nimg, Nreps = Nreps, Nseq = Nseq, Nblocks = Nblocks, stimstart = stimstart, lenstim = lenstim, lenpause = lenpause, strength = strength)
          novelimg += 1
        end
      else
        # ensure that the first images of two consecutive blocks differ
        shuffleimg = shuffle(firstimages)
        #println("begin $(shuffleimg[1])")
        while shuffleimg[1] == storefirstimages[end] && Nseq > 1 # added Nseq larger 1 cause otgerwise always the same
          println("shuffle again")
          shuffleimg = shuffle(firstimages)
          #println("new begin $(shuffleimg[1])")
        end
        storefirstimages = vcat(storefirstimages, shuffleimg)
        #println(storefirstimages)
        for img in shuffleimg  # since we shuffle ensure that it is still the correct novel image
          append!(blockonset, size(stimulus,1) + 1)
          stimulus = gennextsequencenovel(stimulus, img, novelimg ,Nimg = Nimg, Nreps = Nreps, Nseq = Nseq, Nblocks = Nblocks, stimstart = stimstart, lenstim = lenstim, lenpause = lenpause, strength = strength)
          novelimg += 1
        end
      end
    end

    # ensure that blockonset 1 when we start with an emptz array
    if blockonset[1] == 2
      blockonset[1] = 1
    end
    return stimulus, convert(Array{Int64,1}, blockonset)

  end

function gennextsequencenovel(stim::Array{Float64,2}, firstimg, novelimg; Nimg = 4, Nreps = 20, Nseq = 5, Nblocks = 1, stimstart = 1000, lenstim = 1000, lenpause = 3000, strength = 8)
	"""generate a stimulus sequence with novel stim
	Nimg Number of images per sequence
	Nreps Number of seqeunce repititions in one block
	Nseq Number of sequences
	Nblocks Number of block repititions"""
	Nstim = Nimg*Nreps
	tempstim = zeros(Nstim, 4)
    #tempstim[1,1] = firstimg
	# repeat sequence several times, i.e. assemblie numbers 123412341234...
	tempstim[:,1] = repeat(firstimg:(firstimg+Nimg-1), outer = Nreps)
	if stim[end,3] == 0
		tempstim[1,2] = stimstart
	else
    	tempstim[1,2] = stim[end,3]+lenpause
	end
    tempstim[1,3] = tempstim[1,2] + lenstim
    tempstim[1,4] = strength

    for i = 2:Nimg*Nreps #number in image times
        # # lenstim ms stimulation and lenpause ms wait time
        tempstim[i,2] = tempstim[i-1,2] + lenstim + lenpause
        tempstim[i,3] = tempstim[i,2]+lenstim
        tempstim[i,4] = strength
    end
	tempstim[end-Nimg,1] = novelimg
	if stim[end,3] == 0
		stim = tempstim
	else
		stim = vcat(stim, tempstim)
	end
	#tempstim = nothing
    return stim
end # function genstim

function genstimparadigmssa(stimulus; Nimg = 4, Nreps = 20, Nseq = 5, Nblocks = 1, stimstart = 1000, lenstim = 1000, lenpause = 3000, strength = 8)
    """
    generate SSA stimulation paradigm
    A A A A A A A B A A A ... Block 1
    B B B B B 				 Block 2 ...
    Nreps 20 times
    Blocks
    """
    Nimg = 2
    Nreps = 5
    Nseq = 5
    Nblocks = 10
    stimulus, blockonset = genstimparadigmnovelcont(stimulus, Nimg = 2, Nreps = 5, Nseq = 5, Nblocks = 10, stimstart = stimstart, lenstim = 100, lenpause = 300, strength = strength )

    #last_onset = blockonset[1]
    freq = 1
    dev = 2
    for bb =1:length(blockonset)
      #idxfreq = collect(blockonset[bb]:(blockonset[bb]+Nimg*Nreps))
      freq_start = blockonset[bb]
      freq_end = blockonset[bb]+Nimg*Nreps
      dev_end = blockonset[bb]+Nimg*Nreps
      dev_start = blockonset[bb]+Nreps
      if bb == length(blockonset)
        stimulus[freq_start:end,1] .= freq
        #println(idxfreq)
        #idxdev = collect(blockonset[bb]+Nreps:Nreps:(blockonset[bb]+Nimg*Nreps))
        stimulus[dev_start:Nreps:end,1] .= dev
        #printl
      else
        stimulus[freq_start:freq_end,1] .= freq
        #println(idxfreq)
        #idxdev = collect(blockonset[bb]+Nreps:Nreps:(blockonset[bb]+Nimg*Nreps))
        stimulus[dev_start:Nreps:dev_end,1] .= dev
        #println(idxdev)
        # for the nect block switch deviant and frequent
      end
      temp =freq
      freq = dev
      dev = temp
    end

    return stimulus, blockonset
end

function initweights((Ne, Ni, p, Jee0, Jei0, Jie, Jii))
	""" initialise the connectivity weight matrix based on initial
	E-to-E Jee0 initial weight e to e plastic
	I-to-E Jei0 initial weight i to e plastic
	E-to-I Jie0 constant weight e to i not plastic
	I-to-I Jii constant weight i to i not plastic"""
	Ncells = Ne+Ni

	# initialise weight matrix
	# w[i,j] is the weight from pre- i to postsynaptic j neuron
	weights = zeros(Float64,Ncells,Ncells)
	weights[1:Ne,1:Ne] .= Jee0
	weights[1:Ne,(1+Ne):Ncells] .= Jie
	weights[(1+Ne):Ncells,1:Ne] .= Jei0
	weights[(1+Ne):Ncells,(1+Ne):Ncells] .= Jii
	# set diagonal elements to 0
	# for cc = 1:Ncells
	# 	weights[cc,cc] = 0
	# end
	weights[diagind(weights)] .= 0.0
	# ensure that the connection probability is only p
	weights = weights.*(rand(Ncells,Ncells) .< p)
	return weights
end

function weightpars(;Ne = 4000, Ni = 1000, p = 0.2 )
	"""Ne, Ni number of excitatory, inhibitory neurons
	p initial connection probability"""
	Jee0 = 2.86 #initial weight e to e plastic
	Jei0 = 48.7 #initial weight i to e plastic
	Jie = 1.27 #constant weight e to i not plastic
	Jii = 16.2 #constant weight i to i not plastic
	return Ne,Ni,p,Jee0,Jei0,Jie,Jii
end

function initassemblymembers(;Nassemblies = 20, pmember = .05, Nmembersmax = 300, Ne = 4000)
	"""Nassemblies number of assemblies
	pmember probability of belonging to any assembly
	Nmembersmax maximum number of neurons in a population (to set size of matrix)

	set up excitatory assemblies"""
	#seed = 1
	#Random.seed!(seed) # ensure same Assemblies when repeating
	assemblymembers = ones(Int,Nassemblies,Nmembersmax)*(-1)
	for pop = 1:Nassemblies
		members = findall(rand(Ne) .< pmember)
		assemblymembers[pop,1:length(members)] = members
	end
	#println(assemblymembers)
	return assemblymembers
end

function initinhibassemblymembers(;Nassemblies = 20, pmember = .15, Nmembersmax = 200, Ne = 4000, Ni = 1000)
	"""Nassemblies number of assemblies
	pmember probability of belonging to any assembly
	Nmembersmax maximum number of neurons in a population (to set size of matrix)

	set up inhibitory assemblies
	higher connection probability"""
	#seed = 1
	#Random.seed!(seed) # ensure same Assemblies when repeating
	inhibmembers = ones(Int,Nassemblies,Nmembersmax)*(-1)
	for pop = 1:Nassemblies
		members = findall(rand(Ni) .< pmember) .+ Ne
		inhibmembers[pop,1:length(members)] = members
	end
	#println(assemblymembers)
	return inhibmembers
end

#include("../simulation/runsimulation_inhibtuning.jl")
#include("../simulation/sequencefunctions.jl")
#include("../simulation/helperfunctions.jl")
#include("../evaluation/evaluationfunctions.jl")

# --------------------- initialise simulation --------------------------------------------

# Define number of excitatory and inhibitory neurons
const Ne = 4000
const Ni = 1000

Ncells = Ne + Ni

weights_au = initweights(weightpars(Ne = Ne, Ni = Ni))
weights_compare = permutedims(weights_au)
heatmap(weights_au)

pweights=0.2
Jee = 2.86 #initial weight e to e plastic
weights_ee = S.make_sparse_weights(Ne,Ne,pweights,Jee)
Jei = 48.7 #initial weight i to e plastic
weights_ei = S.make_sparse_weights(Ne,Ni,pweights,Jei)
Jie = 1.27 #constant weight e to i not plastic
weights_ie = S.make_sparse_weights(Ni,Ne,pweights,Jie)
Jii = 16.2 #constant weight i to i not plastic
weights_ii = S.make_sparse_weights(Ni,Ni,pweights,Jii)

weights_all = Matrix(vcat(hcat(weights_ee,weights_ei), hcat(weights_ie,weights_ii)))

heatmap(weights_all)
heatmap(weights_compare)

##
# Nimg 3, Nreps 20, Nseq 10, Nblocks 1, lenstim 300, strength 12, Ntrain 5, adjustfactor 1.0, adaptive neurons false, inhibfactor 0.1
ARGS = ["3", "20", "10", "1","300", "12", "5", "10", "0", "10"] # unique sequences
# ARGS = ["3", "20", "5", "10","300", "12", "5", "10", "0", "10"] # repeated sequences

# stimulus parameters
Nimg = parse(Int64, ARGS[1]) # number of stimuli per sequence
Nreps = parse(Int64, ARGS[2]) # number of repetitions per sequence cycle
Nseq = parse(Int64, ARGS[3]) # number of sequences
Nblocks = parse(Int64, ARGS[4]) # number of sequence block repetitions
stimstart = 1_000  # 4000 # start time of stimulation in ms
lenstimpre = parse(Int64, ARGS[5]) # duration of assembly stimulation per image in ms
lenpausepre = 0 # duration of no stimulation between two images in ms
strengthpre = parse(Int64, ARGS[6]) # strength of the stimulation in kHz added to baseline input
lenstim = lenstimpre # duration of assembly stimulation per image in ms
lenpause = lenpausepre # duration of no stimulation between two images in ms
strength = strengthpre # strength of the stimulation in kHz added to baseline input
Ntrain = parse(Int64, ARGS[7]) # number of pretraining iterations
Nass = (Nimg) * Nseq +  Nseq * Nblocks # total number of assemblies Nimg * Nseq + all novelty assemblies
##

p_ase = 0.05
p_asi = 0.15
Nassemblies = 40

assemblymembers = initassemblymembers(Nassemblies = Nass,Ne = Ne)
# inhibitory tuning
inhibassemblies = initinhibassemblymembers(Nassemblies = Nass, Ne = Ne, Ni = Ni)

assemblies_e = map(1:Nassemblies) do _
  sort(findall(rand(Ne) .< p_ase))
end
assemblies_i = map(1:Nassemblies) do _
  sort(findall(rand(Ni) .< p_asi))
end

map(r->count(r.>0),eachrow(assemblymembers))
length.(assemblies_e)

map(r->count(r.>0),eachrow(inhibassemblies))
length.(assemblies_i)

##

stimulus = zeros(1, 4)	# initialisation of the stimulus

# ---------------------- generate the stimulus --------------------------------

# specify aspects of the stimulation paradigm
withnovelty 	= true 	# novelty response
pretrainig 		= true	# include pretraining phase

# ----------- standard setting when the following three false -----------------
shuffledstim 	= false	# shuffle the stimulus randomly
lastshuffled 	= false	# shuffle the last two stimuli
reducednovelty 	= false	# reduce the stimulation strength of the novel stimulus

if withnovelty
	if pretrainig
		# --------- main stim generation function ---------------
		stimulus = genstimparadigmpretraining(stimulus, Nass = Nass, Ntrain = Ntrain, stimstart = stimstart, lenstim = lenstimpre, lenpause = lenpausepre, strength = strengthpre)
		lenpretrain = size(stimulus,1)
	end

	if shuffledstim
		if pretrainig
			stimulus, blockonset = genstimparadigmnovelshuffledpretrain(stimulus, lenpretrain, Nimg = Nimg, Nreps = Nreps, Nseq = Nseq, Nblocks = Nblocks, stimstart = stimstart, lenstim = lenstim, lenpause = lenpause, strength = strength)
			tag = "repeatedsequencesshuffled.h5"

		else
			stimulus, blockonset = genstimparadigmnovelshuffled(stimulus, Nimg = Nimg, Nreps = Nreps, Nseq = Nseq, Nblocks = Nblocks, stimstart = stimstart, lenstim = lenstim, lenpause = lenpause, strength = strength )
			tag = "repeatedsequencesshuffled.h5"
		end
	else # if shuffled stim

		# --------- main stim generation function ---------------
		stimulus, blockonset = genstimparadigmnovelcont(stimulus, Nimg = Nimg, Nreps = Nreps, Nseq = Nseq, Nblocks = Nblocks, stimstart = stimstart, lenstim = lenstim, lenpause = lenpause, strength = strength )
		tag = "repeatedsequences.h5"

	end
else # if with novelty
	# even if no novelty assemblies are shown still train with all novelty assemblies
	if pretrainig
		#println(stimulus)
		stimulus = genstimparadigmpretraining(stimulus, Nass = Nass, Ntrain = Ntrain, stimstart = stimstart, lenstim = lenstimpre, lenpause = lenpausepre, strength = strengthpre)
		#println(stimulus)
		lenpretrain = size(stimulus,1)
		tag = "repeatedsequences_imprinting.h5"

	end
	tag = "repeatedsequencesnoNovelty.h5"
	stimulus, blockonset = genstimparadigmNonovel(stimulus, Nimg = Nimg, Nreps = Nreps, Nseq = Nseq, Nblocks = Nblocks, stimstart = stimstart, lenstim = lenstim, lenpause = lenpause, strength = strength )
end

function inputfun(t::R,i::Integer,low::R,
    stimulus_mat::Matrix{Float64},
    assembly_v::Vector{Vector{Int64}};scale_factor::Float64=1.0) where R<:Real
  # catch the first or last
  (!(stimulus_mat[1,2] <= t <= stimulus_mat[end,3])) && (return low)
  stimstarts = view(stimulus_mat,:,2)
  # find row number
  therow = searchsortedfirst(stimstarts,t) - 1
  # account for stim duration
  (t >= stimulus_mat[therow,3]) && (return low)
  # ok, the assembly is active!
  as_v = assembly_v[Int64(stimulus_mat[therow,1])]
  if _is_in_assembly(i,as_v)
    return low + scale_factor*stimulus_mat[therow,4] # add to stimulus
  else
    return low
  end
  return low # in case the compiler is dumb
end

function inputfun_upper(t::R,i::Integer,low::R,stimmat::Matrix{R},asv::Vector;
   scale_factor::Float64=1.0) where R<:Real
  return max(inputfun(t,i,low,stimmat,asv;scale_factor=scale_factor),
              inputfun(t+2dt,i,low,stimmat,asv;scale_factor=scale_factor))
end

function _is_in_assembly(neu_idx::Integer,as_vec::Vector{<:Integer})
  if as_vec[1] <= neu_idx <= as_vec[end]
    idx = searchsortedfirst(as_vec,neu_idx)
    return as_vec[idx] == neu_idx
  else
    return false
  end
end
const assembly_vectors =
    map(eachrow(hcat(assemblymembers,inhibassemblies))) do r
    sort(filter(>(0),r))
end
as_ordering = let Ntot=Ncells
  S.order_by_pattern_idxs(assembly_vectors,Ntot)
end

const stimulus_sec = copy(stimulus)
stimulus_sec[:,[2,3]] ./= 1E3
stimulus_sec[:,4] .*= 1E3

_ = let low_e = 2E3, low_i = 3E3,
  _scal = 0.1
  ratefun = function (t,i)
    if i<=Ne
      inputfun(t,i,low_e,stimulus_sec,assembly_vectors)
    else
      inputfun(t,i,low_i,stimulus_sec,assembly_vectors,scale_factor=_scal)
    end
  end
  times = range(0,100;length=200)
  nshow = Ncells
  neushow = as_ordering[1:nshow]
  vals = hcat([ ratefun.(t,neushow) for t in times]...)
  # vals_toplot = vals[as_ordering,:][1:nshow,:]
  heatmap(times,1:nshow,vals)
end

# now mine!

##
Ttot = 60.0
lowrate_e =2E3
highrate_e =12E3
lowrate_i =3E3
highrate_i = 12E3
scale_fact_i = 0.1
Δt_as = 0.3
t_as_delay = 1.0
Δt_blank = 0.0

inputfun_e = S.PatternPresentation(Δt_as,Ttot,lowrate_e,lowrate_e+highrate_e,Ne,
    assemblies_e;t_pattern_delay=t_as_delay,Δt_pattern_blank=Δt_blank)
inputfun_i = S.PatternPresentation(Δt_as,Ttot,lowrate_i,lowrate_i+scale_fact_i*highrate_i,
  Ni, assemblies_i;t_pattern_delay=t_as_delay,Δt_pattern_blank=Δt_blank)
# force same sequence
inputfun_i.sequence .= inputfun_e.sequence

##

as_ordering_mine = let
  assemblies_joint=map(ei->sort(vcat(ei[1],Ne.+ei[2])),zip(assemblies_e,assemblies_i))
  S.order_by_pattern_idxs(assemblies_joint,Ne+Ni)
end

##

_ = let
  ratefun = function (t,i)
    if i<=Ne
      inputfun_e(t,i)
    else
      inputfun_i(t,i-Ne)
    end
  end
  times = range(0,100;length=200)
  nshow = Ncells
  neushow = as_ordering_mine[1:nshow]
  vals = hcat([ ratefun.(t,neushow) for t in times]...)
  # vals_toplot = vals[as_ordering,:][1:nshow,:]
  heatmap(times,1:nshow,vals)
end


text:
using Test
using DandelionWebSockets: OPCODE_BINARY

mutable struct FakeTLSStream <: IO
    buf::IOBuffer
    write_buf::IOBuffer
    isopen::Bool
    iseof::Bool

    FakeTLSStream() = new(IOBuffer(), IOBuffer(), true, false)
end

test_write(s::FakeTLSStream, frame::Frame) = write(s.buf, frame)
function Base.write(s::FakeTLSStream, frame::Frame)
    mark(s.write_buf)
    write(s.write_buf, frame)
    reset(s.write_buf)
end

Base.write(s::FakeTLSStream, t::UInt8) = write(s.write_buf, t)
Base.write(s::FakeTLSStream, t::UInt16) = write(s.write_buf, t)
Base.write(s::FakeTLSStream, t::UInt64) = write(s.write_buf, t)

Base.read(::FakeTLSStream, ::T) where T = throw(ErrorException())
Base.readavailable(s::FakeTLSStream) = take!(s.buf)
Base.eof(s::FakeTLSStream) = s.iseof

Base.close(s::FakeTLSStream) = s.isopen = false

@testset "TLSBUfferedIO          " begin
    # MbedTLS.SSLContext does not allow you to read bytes from the IO stream. It throws an exception
    # if you try. `BufferedReader` adapts the TLS stream to support byte I/O by reading all
    # available data into a buffer and return data from that.
    @testset "Read two bytes via TLSBufferedIO" begin
        fake_tls = FakeTLSStream()
        write(fake_tls.buf, Vector{UInt8}([1,2]))

        s = DandelionWebSockets.TLSBufferedIO(fake_tls)

        @test read(s, UInt8) == 1
        @test read(s, UInt8) == 2
    end

    @testset "Read frames via TLSBufferedIO" begin
        frame1 = Frame(true, false, false, false, OPCODE_BINARY, false, 126, 256, b"", ones(UInt8, 256))
        frame2 = Frame(true, false, false, false, OPCODE_BINARY, false, 126, 256, b"", zeros(UInt8, 256))
        fake_tls = FakeTLSStream()
        test_write(fake_tls, frame1)
        test_write(fake_tls, frame2)

        s = DandelionWebSockets.TLSBufferedIO(fake_tls)
        @test read(s, Frame) == frame1
        @test read(s, Frame) == frame2
    end

    @testset "Write frames via TLSBufferedIO" begin
        frame1 = Frame(true, false, false, false, OPCODE_BINARY, false, 126, 256, b"", ones(UInt8, 256))
        frame2 = Frame(true, false, false, false, OPCODE_BINARY, false, 126, 256, b"", zeros(UInt8, 256))
        fake_tls = FakeTLSStream()
        s = DandelionWebSockets.TLSBufferedIO(fake_tls)

        mark(fake_tls.write_buf)
        write(s, frame1)
        write(s, frame2)
        reset(fake_tls.write_buf)
        @test read(fake_tls.write_buf, Frame) == frame1
        @test read(fake_tls.write_buf, Frame) == frame2
    end

    @testset "Close a TLSBufferedIO; underlying stream is closed" begin
        fake_tls = FakeTLSStream()
        stream = DandelionWebSockets.TLSBufferedIO(fake_tls)

        close(stream)

        @test fake_tls.isopen == false
    end

    @testset "EOF indicated by eof() method; read throws EOFError" begin
        fake_tls = FakeTLSStream()
        stream = DandelionWebSockets.TLSBufferedIO(fake_tls)
        fake_tls.iseof = true

        @test_throws EOFError read(stream, Frame)
    end
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Types and structures.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    struct ColumnTable

This structure helps to access elements that comply with the column access
specification of Tables.jl.

"""
struct ColumnTable
    data         # .............................................. Original table
    table        # ....................... Table converted using `Tables.column`
    column_names # ................................................ Column names
    size         # ........................................... Size of the table
end

"""
    struct RowTable

This structure helps to access elements that comply with the row access
specification of Tables.jl.

"""
struct RowTable
    data         # .............................................. Original table
    table        # ......................... Table converted using `Tables.rows`
    column_names # ................................................ Column names
    size         # ........................................... Size of the table
end

"""
    PrintInfo{Td,Th,Trn}

This structure stores the information required so that the backends can print
the tables.

"""
@with_kw struct PrintInfo
    data
    header
    id_cols::Vector{Int}
    id_rows::Vector{Int}
    num_rows::Int
    num_cols::Int
    num_printed_cols::Int
    num_printed_rows::Int
    header_num_rows::Int
    header_num_cols::Int
    show_row_names::Bool
    row_names
    row_name_alignment::Symbol
    row_name_column_title::String
    alignment::Vector{Symbol}
    cell_alignment
    formatters
    compact_printing::Bool
    title::String
    title_alignment::Symbol
    header_alignment::Vector{Symbol}
    header_cell_alignment
    cell_first_line_only::Bool
end


text:
module Citation

export deps, cite

using HTTP

import Pkg
PackageSpec = Pkg.Types.PackageSpec

JULIA_CITATION = """
@article{Julia-2017,
    title={Julia: A fresh approach to numerical computing},
    author={Bezanson, Jeff and Edelman, Alan and Karpinski, Stefan and Shah, Viral B},
    journal={SIAM {R}eview},
    volume={59},
    number={1},
    pages={65--98},
    year={2017},
    publisher={SIAM},
    doi={10.1137/141000671}
}
"""

function doi2bib(doi::AbstractString)
    return String(HTTP.get("http://data.crossref.org/$doi", ["Accept" => "application/x-bibtex"]).body)
end

function deps()
    ctx = Pkg.Types.Context()
    pkgs = PackageSpec[]
    Pkg.Operations.load_direct_deps!(ctx, pkgs)
    return pkgs
end

function collect_citations(paths)
    entries = String[]
    for path in paths
        isdir(path) || continue
        project = Pkg.Types.read_project(Pkg.Types.projectfile_path(path))
        if haskey(project.other, "metadata") && haskey(project.other["metadata"], "cite")
            cite = project.other["metadata"]["cite"]
            push!(entries, doi2bib(cite))
        elseif isfile(joinpath(path, "CITATION.bib"))
            cite = read(String, joinpath(path, "CITATION.bib"))
            push!(entries, cite)
        end
    end
    return entries
end

function cite(;cite_julia=false)
    paths = map(Pkg.Operations.source_path, deps())
    entries = collect_citations(paths)
    cite_julia && push!(entries, JULIA_CITATION)
    print(join(entries, "\n"))
end # module

end


text:
using ConservationLawsParticles
using RecursiveArrayTools, DifferentialEquations, Plots

# model
V(t, x) = 3 + (1+sin(3t)) * cos(4x)
Wprime(t, x) = 2 * sign(x) / (abs(x) + 1)
mob(ρ) = max(1 - ρ, 0)
model = SampledModel((V,), ((Wprime,),), (mob,))

# initial condition
x0 = ArrayPartition(gaussian_particles(2, 801))

# time span
tspan = (0., 5.)

# ODE system for the particles
prob = ODEProblem(velocities!, x0, tspan, model)

sol = solve(prob, BS5(); reltol=1e-8, abstol=1e-8)

# plot the particle trajectories
plot(sol, vars=1:20:801; legend=false, color=:blue,
    title="Trajectories", xlabel="time", ylabel="position")
savefig("traffic1.png")

# plot an animation of the density
anim = @animate for t in range(tspan...; step=1/48)
    plot_density(sol(t); color=:blue, legend=false, xlims=(-2,12), ylims=(0,1),
        title="Density", xlabel="position", ylabel="density")
end
gif(anim, "traffic1.gif")


text:
using Test, MPI
using Random
using StaticArrays
using Test

using ClimateMachine
ClimateMachine.init()
using ClimateMachine.Atmos
using ClimateMachine.Orientations
using ClimateMachine.ConfigTypes
using ClimateMachine.Diagnostics
using ClimateMachine.BalanceLaws
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.Mesh.Filters
using ClimateMachine.TemperatureProfiles
using ClimateMachine.Thermodynamics
using ClimateMachine.TurbulenceClosures
using ClimateMachine.VariableTemplates

using CLIMAParameters
using CLIMAParameters.Planet: R_d, cp_d, cv_d, MSLP, grav
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

import ClimateMachine.Mesh.Grids: _x1, _x2, _x3
import ClimateMachine.BalanceLaws: vars_state
import ClimateMachine.VariableTemplates.varsindex

# ------------------------ Description ------------------------- #
# 1) Dry Rising Bubble (circular potential temperature perturbation)
# 2) Boundaries - `All Walls` : Impenetrable(FreeSlip())
#                               Laterally periodic
# 3) Domain - 2500m[horizontal] x 2500m[horizontal] x 2500m[vertical]
# 4) Timeend - 1000s
# 5) Mesh Aspect Ratio (Effective resolution) 1:1
# 7) Overrides defaults for
#               `init_on_cpu`
#               `solver_type`
#               `sources`
#               `C_smag`
# 8) Default settings can be found in `src/Driver/Configurations.jl`
# ------------------------ Description ------------------------- #
function init_risingbubble!(problem, bl, state, aux, localgeo, t)
    (x, y, z) = localgeo.coord

    FT = eltype(state)
    param_set = parameter_set(bl)
    R_gas::FT = R_d(param_set)
    c_p::FT = cp_d(param_set)
    c_v::FT = cv_d(param_set)
    γ::FT = c_p / c_v
    p0::FT = MSLP(param_set)
    _grav::FT = grav(param_set)

    xc::FT = 1250
    yc::FT = 1250
    zc::FT = 1000
    r = sqrt((x - xc)^2 + (y - yc)^2 + (z - zc)^2)
    rc::FT = 500
    θ_ref::FT = 300
    Δθ::FT = 0

    if r <= rc
        Δθ = FT(5) * cospi(r / rc / 2)
    end

    #Perturbed state:
    θ = θ_ref + Δθ # potential temperature
    π_exner = FT(1) - _grav / (c_p * θ) * z # exner pressure
    ρ = p0 / (R_gas * θ) * (π_exner)^(c_v / R_gas) # density
    q_tot = FT(0)
    ts = PhaseEquil_ρθq(param_set, ρ, θ, q_tot)
    q_pt = PhasePartition(ts)

    ρu = SVector(FT(0), FT(0), FT(0))

    #State (prognostic) variable assignment
    e_kin = FT(0)
    e_pot = gravitational_potential(bl.orientation, aux)
    ρe_tot = ρ * total_energy(e_kin, e_pot, ts)
    state.ρ = ρ
    state.ρu = ρu
    state.energy.ρe = ρe_tot
    state.moisture.ρq_tot = ρ * q_pt.tot
end

function config_risingbubble(FT, N, resolution, xmax, ymax, zmax)

    # Choose explicit multirate solver
    ode_solver = ClimateMachine.MultirateSolverType(
        fast_model = AtmosAcousticGravityLinearModel,
        slow_method = LSRK144NiegemannDiehlBusch,
        fast_method = LSRK144NiegemannDiehlBusch,
        timestep_ratio = 10,
    )

    # Set up the model
    T_profile = DryAdiabaticProfile{FT}(param_set)
    C_smag = FT(0.23)
    ref_state = HydrostaticState(T_profile)
    physics = AtmosPhysics{FT}(
        param_set;
        ref_state = ref_state,
        turbulence = SmagorinskyLilly{FT}(C_smag),
    )
    model = AtmosModel{FT}(
        AtmosLESConfigType,
        physics;
        init_state_prognostic = init_risingbubble!,
        source = (Gravity(),),
    )

    # Problem configuration
    config = ClimateMachine.AtmosLESConfiguration(
        "DryRisingBubble",
        N,
        resolution,
        xmax,
        ymax,
        zmax,
        param_set,
        init_risingbubble!,
        solver_type = ode_solver,
        model = model,
    )
    return config
end

function config_diagnostics(driver_config)
    interval = "10000steps"
    dgngrp = setup_atmos_default_diagnostics(
        AtmosLESConfigType(),
        interval,
        driver_config.name,
    )
    return ClimateMachine.DiagnosticsConfiguration([dgngrp])
end
#-------------------------------------------------------------------------
function run_brick_diagostics_fields_test()
    DA = ClimateMachine.array_type()
    mpicomm = MPI.COMM_WORLD
    root = 0
    pid = MPI.Comm_rank(mpicomm)
    npr = MPI.Comm_size(mpicomm)
    toler = Dict(Float64 => 1e-8, Float32 => 1e-4)
    # Working precision
    for FT in (Float32, Float64)
        # DG polynomial order
        N = (4, 5)
        # Domain resolution and size
        Δh = FT(50)
        Δv = FT(50)
        resolution = (Δh, Δh, Δv)
        # Domain extents
        xmax = FT(2500)
        ymax = FT(2500)
        zmax = FT(2500)
        # Simulation time
        t0 = FT(0)
        timeend = FT(1000)

        # Courant number
        CFL = FT(20)

        driver_config = config_risingbubble(FT, N, resolution, xmax, ymax, zmax)
        solver_config = ClimateMachine.SolverConfiguration(
            t0,
            timeend,
            driver_config,
            init_on_cpu = true,
            Courant_number = CFL,
        )
        #------------------------------------------------------------------------

        model = driver_config.bl
        Q = solver_config.Q
        dg = solver_config.dg
        grid = dg.grid
        vgeo = grid.vgeo
        Nel = size(Q.realdata, 3)
        Npl = size(Q.realdata, 1)

        ind = [
            varsindex(vars_state(model, Prognostic(), FT), :ρ)
            varsindex(vars_state(model, Prognostic(), FT), :ρu)
        ]
        _ρ, _ρu, _ρv, _ρw = ind[1], ind[2], ind[3], ind[4]

        x1 = view(vgeo, :, _x1, 1:Nel)
        x2 = view(vgeo, :, _x2, 1:Nel)
        x3 = view(vgeo, :, _x3, 1:Nel)

        fcn0(x, y, z, xmax, ymax, zmax) =
            sin.(pi * x ./ xmax) .* cos.(pi * y ./ ymax) .* cos.(pi * z ./ zmax)    # sample function

        fcnx(x, y, z, xmax, ymax, zmax) =
            cos.(pi * x ./ xmax) .* cos.(pi * y ./ ymax) .*
            cos.(pi * z ./ zmax) .* pi ./ xmax # ∂/∂x
        fcny(x, y, z, xmax, ymax, zmax) =
            -sin.(pi * x ./ xmax) .* sin.(pi * y ./ ymax) .*
            cos.(pi * z ./ zmax) .* pi ./ ymax # ∂/∂y
        fcnz(x, y, z, xmax, ymax, zmax) =
            -sin.(pi * x ./ xmax) .* cos.(pi * y ./ ymax) .*
            sin.(pi * z ./ zmax) .* pi ./ zmax # ∂/∂z

        Q.data[:, _ρ, 1:Nel] .= 1.0 .+ fcn0(x1, x2, x3, xmax, ymax, zmax) * 5.0
        Q.data[:, _ρu, 1:Nel] .=
            Q.data[:, _ρ, 1:Nel] .* fcn0(x1, x2, x3, xmax, ymax, zmax)
        Q.data[:, _ρv, 1:Nel] .=
            Q.data[:, _ρ, 1:Nel] .* fcn0(x1, x2, x3, xmax, ymax, zmax)
        Q.data[:, _ρw, 1:Nel] .=
            Q.data[:, _ρ, 1:Nel] .* fcn0(x1, x2, x3, xmax, ymax, zmax)
        #-----------------------------------------------------------------------
        vgrad = Diagnostics.VectorGradients(dg, Q)
        vort = Diagnostics.Vorticity(dg, vgrad)
        #----------------------------------------------------------------------------
        Ω₁_exact =
            fcny(x1, x2, x3, xmax, ymax, zmax) -
            fcnz(x1, x2, x3, xmax, ymax, zmax)
        Ω₂_exact =
            fcnz(x1, x2, x3, xmax, ymax, zmax) -
            fcnx(x1, x2, x3, xmax, ymax, zmax)
        Ω₃_exact =
            fcnx(x1, x2, x3, xmax, ymax, zmax) -
            fcny(x1, x2, x3, xmax, ymax, zmax)

        err = zeros(FT, 12)

        err[1] = maximum(abs.(fcnx(x1, x2, x3, xmax, ymax, zmax) - vgrad.∂₁u₁))
        err[2] = maximum(abs.(fcny(x1, x2, x3, xmax, ymax, zmax) - vgrad.∂₂u₁))
        err[3] = maximum(abs.(fcnz(x1, x2, x3, xmax, ymax, zmax) - vgrad.∂₃u₁))

        err[4] = maximum(abs.(fcnx(x1, x2, x3, xmax, ymax, zmax) - vgrad.∂₁u₂))
        err[5] = maximum(abs.(fcny(x1, x2, x3, xmax, ymax, zmax) - vgrad.∂₂u₂))
        err[6] = maximum(abs.(fcnz(x1, x2, x3, xmax, ymax, zmax) - vgrad.∂₃u₂))

        err[7] = maximum(abs.(fcnx(x1, x2, x3, xmax, ymax, zmax) - vgrad.∂₁u₃))
        err[8] = maximum(abs.(fcny(x1, x2, x3, xmax, ymax, zmax) - vgrad.∂₂u₃))
        err[9] = maximum(abs.(fcnz(x1, x2, x3, xmax, ymax, zmax) - vgrad.∂₃u₃))

        err[10] = maximum(abs.(vort.Ω₁ - Ω₁_exact))
        err[11] = maximum(abs.(vort.Ω₂ - Ω₂_exact))
        err[12] = maximum(abs.(vort.Ω₃ - Ω₃_exact))

        errg = MPI.Allreduce(err, max, mpicomm)
        @test maximum(errg) < toler[FT]
    end
end
#----------------------------------------------------------------------------
@testset "Diagnostics Fields tests" begin
    run_brick_diagostics_fields_test()
end
#------------------------------------------------


text:
"""
    FDSNRequest

An abstract type representing requests conforming to the FDSN Web Services specification.

Current subtypes of `FDSNRequest`:
- [`FDSNDataSelect`](@ref)
- [`FDSNEvent`](@ref)
- [`FDSNStation`](@ref)
"""
abstract type FDSNRequest <: SeisRequest end

protocol_string(::FDSNRequest) = "fdsnws"

Base.broadcastable(request::FDSNRequest) = Ref(request)

"""
    FDSNEvent(kwargs...)

Create a FDSN web services event query which can be sent to a datacentre which implements
the FDSN Web Services specification.

## Examples
A request asking for events above magnitude 3 within 5° distance of the
Eiffel Tower in the last five years:
```
julia> using Dates

julia> FDSNEvent(starttime=now()-Year(5), longitude=2.295, latitude=48.858, maxradius=5, minmagnitude=3)
```

## Available options

|Parameter|Description|
|:--------|:----------|
|`starttime`|Limit to events on or after the specified start time.|
|`endtime`|Limit to events on or before the specified end time.|
|`minlatitude`|Limit to events with a latitude larger than or equal to the specified minimum.|
|`maxlatitude`|Limit to events with a latitude smaller than or equal to the specified maximum.|
|`minlongitude`|Limit to events with a longitude larger than or equal to the specified minimum.|
|`maxlongitude`|Limit to events with a longitude smaller than or equal to the specified maximum.|
|`latitude`|Specify the latitude to be used for a radius search.|
|`longitude`|Specify the longitude to the used for a radius search.|
|`minradius`|Limit to events within the specified minimum number of degrees from the geographi point defined by the latitude and longitude parameters.|
|`maxradius`|Limit to events within the specified maximum number of degrees from the geographic point defined by the latitude and longitude parameters.|
|`mindepth`|Limit to events with depth more than the specified minimum.|
|`maxdepth`|Limit to events with depth less than the specified maximum.|
|`minmagnitude`|Limit to events with a magnitude larger than the specified minimum.|
|`maxmagnitude`|Limit to events with a magnitude smaller than the specified maximum.|
|`magnitudetype`|Specify a magnitude type to use for testing the minimum and maximum limits.|
|`includeallorigins`|Specify if all origins for the event should be included, default is data center dependent but is suggested to be the preferred origin only.|
|`includeallmagnitudes`|Specify if all magnitudes for the event should be included, default is data center dependent but is suggested to be the preferred magnitude only.|
|`includearrivals`|Specify if phase arrivals should be included.|
|`eventid`|Select a specific event by ID; event identifiers are data center specific.|
|`limit`|Limit the results to the specified number of events.|
|`offset`|Return results starting at the event count specified, starting at 1.|
|`orderby`|Order the result by time or magnitude with the following possibilities: `"time"`: order by origin descending time; `"time-asc"`: order by origin ascending time; `"magnitude"`: order by descending magnitude; `"magnitude-asc"`: order by ascending magnitude|
|`catalog`|Limit to events from a specified catalog|
|`contributor`|Limit to events contributed by a specified contributor.|
|`updatedafter`|Limit to events updated after the specified time.|
|`format`|Specify format of result, either `"xml"` (default) or `"text"` (`"isf"` from the ISC).  If this parameter is not specified the service must return QuakeML.|
|`nodata`|Select status code for “no data”, either `204` (default) or `404`.|

## Specifying dates
Both `Dates.DateTime` objects and `String`s can be passed as values to the
`starttime`, `endtime` and `updatedafter` keywords so long as they are valid.
Examples include `"2000-01-01"` and `"2008-02-03T00:00:30"`.
"""
@with_kw struct FDSNEvent <: FDSNRequest
    starttime::MDateTime = missing
    endtime::MDateTime = missing
    minlatitude::MFloat = missing
    maxlatitude::MFloat = missing
    minlongitude::MFloat = missing
    maxlongitude::MFloat = missing
    latitude::MFloat = missing
    longitude::MFloat = missing
    minradius::MFloat = missing
    maxradius::MFloat = missing
    mindepth::MFloat = missing
    maxdepth::MFloat = missing
    minmagnitude::MFloat = missing
    maxmagnitude::MFloat = missing
    magnitudetype::MString = missing
    includeallorigins::MBool = missing
    includeallmagnitudes::MBool = missing
    includearrivals::MBool = missing
    eventid::MString = missing
    limit::MInt = missing
    offset::MInt = missing
    orderby::MString = missing
    catalog::MString = missing
    contributor::MString = missing
    updatedafter::MDateTime = missing
    format::MString = missing
    nodata::Int = 204
    function FDSNEvent(st, et, minlat, maxlat, minlon, maxlon, lat, lon, minr, maxr, mind, maxd,
                       minmag, maxmag, magtype, allorigins, allmags, arrs, id, limit, offset,
                       orderby, catalog, contributor, updatedafter, format, nodata)
        !ismissing(st) && (st = DateTime(st))
        !ismissing(et) && (et = DateTime(et))
        !ismissing(updatedafter) && (updatedafter = DateTime(updatedafter))
        !ismissing(st) && !ismissing(et) && st > et &&
           throw(ArgumentError("`starttime` must be before endtime"))
        (!ismissing(minr) || !ismissing(maxr)) && (ismissing(lat) || ismissing(lon)) &&
            throw(ArgumentError("Both `latitude` and `longitude` must be present if " *
                                "specifying `minradius` or `maxradius`"))
        coalesce(limit, 1) > 0 || throw(ArgumentError("`limit` cannot be less than 1"))
        coalesce(offset, 1) > 0 || throw(ArgumentError("`offset` cannot be less than 1"))
        !ismissing(orderby) && orderby ∉ ("time", "time-asc", "magnitude", "magnitude-asc")
        -90 ≤ coalesce(minlat, -90) ≤ 90 || throw(ArgumentError("`minlatitude` must be in range -90 to 90"))
        -90 ≤ coalesce(maxlat, 90) ≤ 90 || throw(ArgumentError("`maxlatitude` must be in range -90 to 90"))
        -180 ≤ coalesce(minlon, -180) ≤ 180 || throw(ArgumentError("`minlongitude` must be in range -180 to 180"))
        -180 ≤ coalesce(maxlon, 180) ≤ 180 || throw(ArgumentError("`maxlongitude` must be in range -180 to 180"))
        -90 ≤ coalesce(lat, 0) ≤ 90 || throw(ArgumentError("`latitude` must be in range -90 to 90"))
        -180 ≤ coalesce(lon, 0) ≤ 180 || throw(ArgumentError("`longitude` must be in range -180 to 180"))
        0 ≤ coalesce(minr, 0) ≤ 180 || throw(ArgumentError("`minradius` must be in range 0 to 180"))
        0 ≤ coalesce(maxr, 180) ≤ 180 || throw(ArgumentError("`maxradius` must be in range 0 to 180"))
        coalesce(format, "xml") ∈ ("xml", "isf", "miniseed", "text") ||
            throw(ArgumentError("`format` must be one of \"xml\", \"isf\" \"miniseed\" or \"text\""))
        nodata ∈ (204, 404) || throw(ArgumentError("`nodata` must be 204 or 404"))
        new(st, et, minlat, maxlat, minlon, maxlon, lat, lon, minr, maxr, mind, maxd,
            minmag, maxmag, magtype, allorigins, allmags, arrs, id, limit, offset,
            orderby, catalog, contributor, updatedafter, format, nodata)
    end
end

"""
    FDSNDataSelect(kwargs...)

Create a FDSN web services data selection query which can be sent to a datacentre which
implements the FDSN Web Services specification.

Set options with keyword arguments.

## Example
Request for data in miniseed format (the default) for station A07E in network XM
for all channels on the first day of 2013.
```
julia> using Dates

julia> FDSNDataSelect(network="XM", station="A07E", starttime=DateTime(2013), endtime=DateTime(2013)+Day(1))

## Available options

|Parameter|Description|
|:--------|:----------|
|`starttime`|Limit results to time series samples on or after the specified start time.|
|`endtime`|Limit results to time series samples on or before the specified end time.|
|`network`|Select one or more network codes.  Can be SEED network codes or data center defined codes.  Multiple codes are comma-separated.|
|`station`|Select one or more SEED station codes.  Multiple codes are comma-separated.|
|`location`|Select one or more SEED location identifiers.  Multiple identifiers are comma-separated.  As a special case “--“ (two dashes) will be translated to a string of two space characters to match blank location IDs.|
|`channel`|Select one or more SEED channel codes.  Multiple codes are comma-separated.|
|`code`|SEED-style code, used instead of `network`, `station`, `location` and `channel`.  (Addition to standard in this package.)|
|`quality`|Select a specific SEED quality indicator, handling is data center dependent.|
|`minimumlength`|Limit results to continuous data segments of a minimum length specified in seconds.|
|`longestonly`|Limit results to the longest continuous segment per channel.|
|`format`|Specify format of result, the default value is `"miniseed"`.|
|`nodata`|Select status code for “no data”, either ‘204’ (default) or ‘404’.|

## Specifying dates
Both `Dates.DateTime` objects and `String`s can be passed as values to the
`starttime` and `endtime` keywords so long as they are valid.
Examples include `"2000-01-01"` and `"2008-02-03T00:00:30"`.

## Specifying station channel codes
As well as the options specified by the FDSN standard, the additional option
`code` is provided for ease of use.  This accepts channel codes in the SEED
convention `"⟨network⟩.⟨station⟩.⟨channel⟩.⟨location⟩"`, where each `.`-separated
field can contain wildcards and `,`-separated lists of codes.
"""
struct FDSNDataSelect <: FDSNRequest
    starttime::MDateTime
    endtime::MDateTime
    network::MString
    station::MString
    location::MString
    channel::MString
    quality::MString
    minimumlength::MFloat
    longestonly::MBool
    format::MString
    nodata::Int
    function FDSNDataSelect(st, et, net, sta, loc, cha, q, minlength, longestonly, format, nodata)
        !ismissing(st) && (st = DateTime(st))
        !ismissing(et) && (et = DateTime(et))
        !ismissing(st) && !ismissing(et) && st > et &&
            throw(ArgumentError("`starttime` must be before endtime"))
        !ismissing(net) && !isascii(net) && throw(ArgumentError("`network` must be ASCII"))
        !ismissing(sta) && !isascii(sta) && throw(ArgumentError("`station` must be ASCII"))
        if !ismissing(loc)
            !isascii(loc) && throw(ArgumentError("`location` must be ASCII"))
            loc == "  " && (loc = "--")
        end
        !ismissing(cha) && !isascii(cha) && throw(ArgumentError("`channel` must be ASCII"))
        coalesce(q, "B") ∈ ("D", "R", "Q", "M", "B") ||
            throw(ArgumentError("`quality` must be one of \"D\", \"R\", \"Q\", \"M\" or \"B\""))
        coalesce(minlength, 0) >= 0 || throw(ArgumentError("`minimumlength` must be 0 or more"))
        nodata ∈ (204, 404) || throw(ArgumentError("`nodata` must be 204 or 404"))
        new(st, et, net, sta, loc, cha, q, minlength, longestonly, format, nodata)
    end
end

function FDSNDataSelect(; code=missing, starttime=missing, endtime=missing,
        network=missing, station=missing, location=missing, channel=missing,
        quality=missing, minimumlength=missing, longestonly=missing, format=missing,
        nodata=204)
    if !ismissing(code)
        any(!ismissing, (network, station, location, channel)) &&
            throw(ArgumentError("`code` cannot be provided with any of " *
                                "`network`, `station`, `location`, or `channel`"))
        network, station, location, channel = split_channel_code(code)
    end
    FDSNDataSelect(starttime, endtime, network, station, location, channel, quality,
        minimumlength, longestonly, format, nodata)
end

"""
    FDSNStation(kwargs...)

Create a FDSN web services station query which can be sent to a datacentre which implements
the FDSN Web Services specification.

## Example

Information about all stations called `"ANMO"`, in all networks, at the
default `"station"` level.
```
julia> FDSNStation(station="ANMO")
```

Information about the channels, including sensor response information,
for the broadband, high-sensitivity channels of station JSA in Jersey,
Channel Islands
```
julia> FDSNStation(network="GB", station="JSA", channel="BH?", level="response")
```

Every station known by the datacentre.  Here we use the `format="text"` option
to get only the basic information on where stations are.  (Note that as of
2020, sending this to the default datacentre, IRIS, returns 50,000 stations.
Amazingly, the entire thing takes only a few seconds.)
```
julia> FDSNStation(format="text")
```

Use the `code` option to find all very long period channels in the
Global Seismograph Network active since 2018.
```
julia> using Dates

julia> FDSNStation(code="II,IU.*.*.V??", starttime=DateTime(2018))
```

## Default information level

`FDSNStation` requests by default are returned with detail on the level of
individual stations.  (See table below.)  However, if the `channel` keyword
is specified, then it is assumed that information on the channel level is
desired and the `level` option is set to `"channel"` upon construction.
Set `level` manually to override this default.

## Available options

|Parameter|Description|
|:--------|:----------|
|`starttime`|Limit to metadata epochs starting on or after the specified start time.|
|`endtime`|Limit to metadata epochs ending on or before the specified end time.|
|`startbefore`|Limit to metadata epochs starting before specified time.|
|`startafter`|Limit to metadata epochs starting after specified time.|
|`endbefore`|Limit to metadata epochs ending before specified time.|
|`endafter`|Limit to metadata epochs ending after specified time.|
|`network`|Select one or more network codes.  Can be SEED network codes or data center defined codes.  Multiple codes are comma-separated.|
|`station`|Select one or more SEED station codes.  Multiple codes are comma-separated.|
|`location`|Select one or more SEED location identifiers.  Multiple identifiers are comma-separated.  As a special case “--“ (two dashes) will be translated to a string of two space characters to match blank location IDs.|
|`channel`|Select one or more SEED channel codes.  Multiple codes are comma-separated.|
|`code`|SEED-style code, used instead of `network`, `station`, `location` and `channel`.  (Addition to standard in this package.)|
|`minlatitude`|Limit to stations with a latitude larger than or equal to the specified minimum.|
|`maxlatitude`|Limit to stations with a latitude smaller than or equal to the specified maximum.|
|`minlongitude`|Limit to stations with a longitude larger than or equal to the specified minimum.|
|`maxlongitude`|Limit to stations with a longitude smaller than or equal to the specified maximum.|
|`latitude`|Specify the latitude to be used for a radius search.|
|`longitude`|Specify the longitude to the used for a radius search.|
|`minradius`|Limit results to stations within the specified minimum number of degrees from the geographic point defined by the latitude and longitude parameters.|
|`maxradius`|Limit results to stations within the specified maximum number of degrees from the geographic point defined by the latitude and longitude parameters.|
|`level`|Specify the level of detail for the results, one of `"network`, `"station"` (default), `"channel"` or `"response"`.  If `channel` option is supplied, then `"channel"` is the default.|
|`includerestricted`|Specify if results should include information for restricted stations.|
|`includeavailability`|Specify if results should include information about time series data availability.|
|`updatedafter`|Limit to metadata updated after specified time; updates are data center specific.|
|`matchtimeseries`|Limit to metadata where selection criteria matches time series data availability.|
|`format`|Specify format of result, either `"xml"` (default) or `"text"`.|
|`nodata`|Select status code for “no data”, either `204` (default) or `404`.|

## Specifying dates
Both `Dates.DateTime` objects and `String`s can be passed as values to the
`starttime`, `endtime`, `startbefore`, `startafter`, `endbefore`, `endafter`
and `updatedafter` keywords so long as they are valid.
Examples include `"2000-01-01"` and `"2008-02-03T00:00:30"`.

## Specifying station channel codes
As well as the options specified by the FDSN standard, the additional option
`code` is provided for ease of use.  This accepts channel codes in the SEED
convention `"⟨network⟩.⟨station⟩.⟨channel⟩.⟨location⟩"`, where each `.`-separated
field can contain wildcards and `,`-separated lists of codes.
"""
struct FDSNStation <: FDSNRequest
    starttime::MDateTime
    endtime::MDateTime
    startbefore::MDateTime
    startafter::MDateTime
    endbefore::MDateTime
    endafter::MDateTime
    network::MString
    station::MString
    location::MString
    channel::MString
    minlatitude::MFloat
    maxlatitude::MFloat
    minlongitude::MFloat
    maxlongitude::MFloat
    latitude::MFloat
    longitude::MFloat
    minradius::MFloat
    maxradius::MFloat
    level::MString
    includerestricted::MBool
    includeavailability::MBool
    updatedafter::MDateTime
    matchtimeseries::MBool
    format::MString
    nodata::Int
    function FDSNStation(st, et, sb, sa, eb, ea, net, sta, loc, cha, minlat, maxlat,
                         minlon, maxlon, lat, lon, minr, maxr, level, restricted,
                         availability, updatedafter, matchtimeseries, format, nodata)
        !ismissing(st) && (st = DateTime(st))
        !ismissing(et) && (et = DateTime(et))
        !ismissing(sb) && (sb = DateTime(sb))
        !ismissing(sa) && (sa = DateTime(sa))
        !ismissing(eb) && (eb = DateTime(eb))
        !ismissing(ea) && (ea = DateTime(ea))
        !ismissing(updatedafter) && (updatedafter = DateTime(updatedafter))
        !ismissing(st) && !ismissing(et) && st > et &&
            throw(ArgumentError("`starttime` must be before endtime"))
        !ismissing(sb) && !ismissing(sa) && sb > sa &&
            throw(ArgumentError("`startbefore` must be before startafter"))
        !ismissing(eb) && !ismissing(ea) && eb > ea &&
            throw(ArgumentError("`endbefore` must be before endafter"))
        !ismissing(net) && !isascii(net) && throw(ArgumentError("`network` must be ASCII"))
        !ismissing(sta) && !isascii(sta) && throw(ArgumentError("`station` must be ASCII"))
        if !ismissing(loc)
            !isascii(loc) && throw(ArgumentError("`channel` must be ASCII"))
            loc == "  " && (loc = "--")
        end
        !ismissing(cha) && !isascii(cha) && throw(ArgumentError("`channel` must be ASCII"))
        (!ismissing(minr) || !ismissing(maxr)) && (ismissing(lat) || ismissing(lon)) &&
            throw(ArgumentError("Both `latitude` and `longitude` must be present if " *
                                "specifying `minradius` or `maxradius`"))
        -90 ≤ coalesce(minlat, -90) ≤ 90 || throw(ArgumentError("`minlatitude` must be in range -90 to 90"))
        -90 ≤ coalesce(maxlat, 90) ≤ 90 || throw(ArgumentError("`maxlatitude` must be in range -90 to 90"))
        -180 ≤ coalesce(minlon, -180) ≤ 180 || throw(ArgumentError("`minlongitude` must be in range -180 to 180"))
        -180 ≤ coalesce(maxlon, 180) ≤ 180 || throw(ArgumentError("`maxlongitude` must be in range -180 to 180"))
        -90 ≤ coalesce(lat, 0) ≤ 90 || throw(ArgumentError("`latitude` must be in range -90 to 90"))
        -180 ≤ coalesce(lon, 0) ≤ 180 || throw(ArgumentError("`longitude` must be in range -180 to 180"))
        0 ≤ coalesce(minr, 0) ≤ 180 || throw(ArgumentError("`minradius` must be in range 0 to 180"))
        0 ≤ coalesce(maxr, 180) ≤ 180 || throw(ArgumentError("`maxradius` must be in range 0 to 180"))
        if !ismissing(level)
            !isascii(level) && throw(ArgumentError("`level` must be ASCII"))
            level ∉ ("network", "station", "channel", "response") &&
                throw(ArgumentError("`level` must be one of \"network\", \"station\", \"channel\" or \"response\""))
        else
            # Default to channel level if a channel pattern is given
            ismissing(level) && !ismissing(cha) && (level = "channel")
        end
        nodata ∈ (204, 404) || throw(ArgumentError("`nodata` must be 204 or 404"))
        new(st, et, sb, sa, eb, ea, net, sta, loc, cha, minlat, maxlat,
            minlon, maxlon, lat, lon, minr, maxr, level, restricted,
            availability, updatedafter, matchtimeseries, format, nodata)
    end
end

function FDSNStation(; code=missing,
        starttime=missing, endtime=missing, startbefore=missing, startafter=missing,
        endbefore=missing, endafter=missing, network=missing, station=missing,
        location=missing, channel=missing, minlatitude=missing, maxlatitude=missing,
        minlongitude=missing, maxlongitude=missing, latitude=missing, longitude=missing,
        minradius=missing, maxradius=missing, level=missing, includerestricted=missing,
        includeavailability=missing, updatedafter=missing, matchtimeseries=missing,
        format=missing, nodata=204)
    if !ismissing(code)
        any(!ismissing, (network, station, location, channel)) &&
            throw(ArgumentError("`code` cannot be provided with any of " *
                                "`network`, `station`, `location`, or `channel`"))
        network, station, location, channel = split_channel_code(code)
    end
    FDSNStation(starttime, endtime, startbefore, startafter, endbefore, endafter,
        network, station, location, channel, minlatitude, maxlatitude, minlongitude,
        maxlongitude, latitude, longitude, minradius, maxradius, level,
        includerestricted, includeavailability, updatedafter, matchtimeseries,
        format, nodata)
end

service_string(::FDSNEvent) = "event"
service_string(::FDSNDataSelect) = "dataselect"
service_string(::FDSNStation) = "station"


text:
# This file is a part of RadiationDetectorDSP.jl, licensed under the MIT License (MIT).

const RealOrSIMD{T<:Real} = Union{T,<:Vec{N,<:T} where N}


text:
module BuildBasis

export buildbasisfunc

using TypesParticles, TypesBasis

@doc raw"""
buildbasisfunction(atomicsystem::SystemOfAtoms,basis::Basis)

Returns atom-centered basis functions and number of electrons.

Description: This function will take the system information and create the basis functions
that are atomically centered. Each basis function will contain the parameters for the Gaussian
functions that when linearly combined make up the basis function. Note that the actual basis
function is not yet constructed here, i.e., we are not taking the summation of Gaussians.

Depends on datatypes in TypesParticles.jl and TypesBasis.jl

""" function buildbasisfunc(atomicsystem::SystemOfAtoms,basis::Basis)
    natoms = atomicsystem.natoms;
    nelectrons = 0;

    #Closed-shell system only
    for i=1:natoms
        nelectrons += atomicsystem.atoms[i].Z
    end

    nbasis = basis.nbasisfunc;
    basisfunc = Array{GaussOrbitals,2}(undef,natoms,nbasis);

    #Assign each atom-centered basis functions
    #ie = 1; #total electron counter
    for a=1:natoms
        #For each gaussian store 5 values: c,alpha,xo,yo,zo
        #For Z=1 n-basis per atom is 1 so basisfunc length equal num atoms.
        for j=1:nbasis
            #6-31 type basis set
            if basis.basisfunc[j].info == "6-31"
                basisfunc[a,j] =  GaussOrbitals(atomicsystem.atoms[a],basis.basisfunc[j])
            end
        end # nbasis
    end # atomicsystem.atoms

    return basisfunc,nelectrons
end #buildbasisfunction

end #BuildBasis


text:
using DemoCards: infer_pagedir, is_demosection, is_democard

@testset "DemoPage" begin
    root = joinpath("assets", "page")

    # default DemoPage behavior
    page = DemoPage(joinpath(root, "default"))
    @test page.title == "Default"
    sec1, sec2 = page.sections
    @test sec1.root == joinpath(page.root, "subsection_1")
    @test sec2.root == joinpath(page.root, "subsection_2")

    # title and order
    page = DemoPage(joinpath(root, "title_and_order"))
    @test page.title == "Custom Title"
    sec1, sec2 = page.sections
    @test sec1.root == joinpath(page.root, "subsection_2")
    @test sec2.root == joinpath(page.root, "subsection_1")

    # template
    page = DemoPage(joinpath(root, "template"))
    @test page.title == "Custom Title"

    # template
    page = DemoPage(joinpath(root, "template_2"))
    @test page.title == "Custom Title"

    # template has higher priority
    page = DemoPage(joinpath(root, "suppressed_title"))
    @test page.title == "Custom Title"

    @testset "empty demo section" begin
        page_dir = joinpath("assets", "section", "empty")
        mkpath(joinpath(page_dir, "sec1")) # empty folder cannot be commited into git history

        page = @suppress_err DemoPage(page_dir)
        @test length(page.sections) == 1
        warn_msg = @capture_err DemoPage(page_dir)
        @test occursin("Empty section detected", warn_msg)

        @test_throws ErrorException("Empty demo page, you have to add something.") @suppress_err begin
            preview_demos(joinpath("assets", "section", "empty", "sec1"))
        end
        @test isempty(readdir(joinpath(page_dir, "sec1")))
        rm(joinpath(page_dir, "sec1"))
    end
end

@testset "infer_pagedir" begin
    for pagename in ["default", "suppressed_title", "template", "template_2", "title_and_order"]
        page_dir = joinpath("assets", "page", pagename)

        for (root, dirs, files) in walkdir(page_dir)
            for dir in dirs
                secdir = joinpath(root, dir)
                is_demosection(secdir) && @test infer_pagedir(secdir) == page_dir
            end
            for file in files
                cardpath = joinpath(root, file)
                is_democard(cardpath) && @test infer_pagedir(cardpath) == page_dir
            end
        end
    end

    # when no page structure is detected, return nothing
    @test infer_pagedir("assets") |> isnothing
    @test infer_pagedir(joinpath("assets", "regexes")) |> isnothing
    @test infer_pagedir(joinpath("assets", "regexes", "example_1.jl")) |> isnothing
end


text:
@unroll function max_step(cones::Tuple{Vararg{Cone}}, x)
	maxim = typemin(Float64)
	@unroll for cone in cones
		val = max_step(cone, x)
		if val > maxim
			maxim = val
		end
	end
	return maxim
end

function max_step(cone::POC{dim}, x) where dim
	minim = typemax(Float64)
	for i=cti(cone, 1):cti(cone,dim)
		if x[i] < minim
			minim = x[i]
		end
	end
	return -minim
end

function max_step(cone::SOC{dim}, x) where dim
	sqnrm = 0.0
	for i=cti(cone, 2):cti(cone,dim)
		sqnrm += x[i]^2
	end
	return sqrt(sqnrm) - x[cti(cone, 1)]
end

function compute_step(cones, l, ds, dz)
	mxs = scmax(cones, l, ds) # max_step(cones, scale2(cones, l, ds))
	mxz = scmax(cones, l, dz) # max_step(cones, scale2(cones, l, dz))
	t = max(mxs,mxz,0.0)
	if t == 0.0
		step = 1.0
	else
		step = min(1.0, 1.0/t)
	end
	return step
end

@unroll function scmax(cones::Tuple{Vararg{Cone}}, l, x)
	mxv = typemin(Float64)
	@unroll for cone in cones
		val = scmax(cone, l, x)
		if val > mxv
			mxv = val
		end
	end
	return mxv
end

function scmax(c::POC{dim}, l, x) where dim
	mxv = typemin(Float64)
	for i=cti(c, 1):cti(c, dim)
		val = -x[i]/l[i]
		if val > mxv
			mxv = val
		end
	end
	return mxv
end

function scmax(c::SOC{dim}, li, xi) where dim
	# a = 1/sqrt(ln[1]^2 - sum(ln[2:end].^2))
	i1 = cti(c,1)
	ai = li[i1]^2
	for ii = cti(c,2):cti(c,dim)
		ai -= li[ii]^2
	end
	a = 1/sqrt(ai)

	# r1 = a*ln[1]*x[1] - dot(a*ln[2:end], x[2:end])
	r1 = a*li[i1]*xi[i1]
	for ii = cti(c,2):cti(c,dim)
		r1 -= a*li[ii]*xi[ii]
	end

	# r2 = a .* (x[2:end] - cst * a*ln[2:end])
	cst = (r1 + xi[i1])/(a*li[i1] + 1)
	r2s = 0.0
	for ii = cti(c,2):cti(c,dim)
		r2s += (a * (xi[ii] - cst * a * li[ii]))^2
	end
	return sqrt(r2s) - a * r1 # norm(r2) - a*r1
end


text:
using MLStyle

struct TokenView
    current :: Int
    length  :: Int
    source  :: Vector{Token}
end

TokenView(src :: Vector{Token}) = TokenView(1, length(src), src)

struct CtxTokens{State}
    tokens :: TokenView
    maxfetched :: Int
    state :: State
end

CtxTokens{A}(tokens :: TokenView) where A = CtxTokens(tokens, 1, crate(A))

function inherit_effect(ctx0 :: CtxTokens{A}, ctx1 :: CtxTokens{B}) where {A, B}
    CtxTokens(ctx0.tokens, max(ctx1.maxfetched, ctx0.maxfetched), ctx0.state)
end

orparser(pa, pb) = function (ctx_tokens)
    @match pa(ctx_tokens) begin
        (nothing, ctx_tokens) => pb(ctx_tokens)
        _ && a =>  a
    end
end

tupleparser(pa, pb) = function (ctx_tokens)
    @match pa(ctx_tokens) begin
        (nothing, _) && a => a
            (a, remained) =>
            @match pb(remained) begin
                (nothing, _) && a => a
                (b, remained) => ((a, b), remained)
            end
    end
end

manyparser(p) = function (ctx_tokens)
    res = []
    remained = ctx_tokens
    while true
        (elt, remained) = p(remained)
        if elt === nothing
            break
        end
        push!(res, elt)
    end

    (res, remained)
end

hlistparser(ps) = function (ctx_tokens)
    hlist = Vector{Union{T, Nothing} where T}(nothing, length(ps))
    done = false
    remained = ctx_tokens
    for (i, p) in enumerate(ps)
        @match p(remained) begin
            (nothing, a) =>
                begin
                    hlist = nothing
                    remained = a
                    done = true
                end
            (elt, a) =>
                begin
                    hlist[i] = elt
                    remained = a
                end
        end
        if done
            break
        end
    end
    (hlist, remained)
end

# what is a htuple? ... In fact a tuple with multielements is called a hlist,
# but in dynamic lang things get confusing..
htupleparser(ps) = function (ctx_tokens)
    hlist = Vector{Union{T, Nothing} where T}(nothing, length(ps))
    done = false
    remained = ctx_tokens
    for (i, p) in enumerate(ps)
        @match p(remained) begin
            (nothing, a) =>
                begin
                    hlist = nothing
                    remained = a
                    done = true
                end
            (elt, a) =>
                begin
                    hlist[i] = elt
                    remained = a
                end
        end
        if done
            break
        end
    end
    (hlist !== nothing ? Tuple(hlist) : nothing, remained)
end

trans(f, p) = function (ctx_tokens)
    @match p(ctx_tokens) begin
        (nothing, _) && a => a
        (arg, remained) => (f(arg), remained)
    end
end

optparser(p) = function (ctx_tokens)
    @match p(ctx_tokens) begin
        (nothing, _) => (Some(nothing), ctx_tokens)
        (a, remained) => (Some(a), remained)
    end
end

"""
State.put
"""
updateparser(p, f) = function (ctx_tokens)
    @match p(ctx_tokens) begin
        (nothing, _) && a =>  a
        (a, remained) =>
            let state = f(remained.state, a)
                (a, update_state(remained, state))
            end
    end
end

"""
State.get
"""
getparser(p, f) = function (ctx_tokens)
    @match p(ctx_tokens) begin
        (nothing, _) && a =>  a
        (_, remained) => (f(remained.state), remained)
    end
end

tokenparser(f) = function (ctx_tokens)
    let tokens = ctx_tokens.tokens
        if tokens.current <= tokens.length
            token = tokens.source[tokens.current]
            if f(token)
                new_tokens = TokenView(tokens.current + 1, tokens.length, tokens.source)
                max_fetched = max(new_tokens.current, ctx_tokens.maxfetched)
                new_ctx = CtxTokens(new_tokens, max_fetched, ctx_tokens.state)
                return (token, new_ctx)
            end
        end
        (nothing, ctx_tokens)
    end
end

function direct_lr(init, trailer, reducer)
    function (ctx_tokens)
        res = nothing
        remained = ctx_tokens
        res, remained = init(remained)
        if res === nothing
            return (nothing, remained)
        end

        while true
            miscellaneous, remained = trailer(remained)
            if miscellaneous === nothing
                break
            end
            res = reducer(res, miscellaneous)
        end
        (res, remained)
    end
end

function update_state(ctx:: CtxTokens{A}, state::B) where {A, B}
    CtxTokens(ctx.tokens, ctx.maxfetched, state)
end

function crate
end

function crate(::Type{Vector{T}}) where T
    T[]
end

function crate(::Type{Vector{T} where T})
    []
end

function crate(::Type{Nothing})
    nothing
end

function crate(::Type{Union{T, Nothing}}) where T
    nothing
end

function crate(::Type{Bool}) where T
    false
end

function crate(::Type{T}) where T <: Number
    zero(T)
end

function crate(::Type{String})
    ""
end

function crate(::Type{Any}) where T
    nothing
end

function crate(::Type{LinkedList})
    nil(Any)
end

function crate(::Type{LinkedList{T}}) where T
    nil(T)
end

# number = mklexer(LexerSpec(r"\G\d+"))
# dio = mklexer(LexerSpec("dio"))
# a = mklexer(LexerSpec('a'))

# tables = [:number => number, :dio => dio, :a => a]
# lexer = @genlex tables

# preda(::Token{:a}) = true
# preda(_) = false

# prednum(::Token{:number}) = true
# prednum(_) = false

# preddio(::Token{:dio}) = true
# preddio(_) = false

# p = hlistparser([tokenparser(preddio), tokenparser(prednum), tokenparser(preda)])
# source = lexer("dio111a")
# println(source)
    # tokens = TokenView(source)
    # ctx = CtxTokens{Nothing}(tokens)
    # println(p(ctx))