text:
module LabMaze

using Random, Plots, Logging, DataStructures

include("exports.jl")
include("shapes.jl")
include("boundary.jl")
include("cell.jl")
include("grid.jl")
include("directions.jl")
include("neighbors.jl")
include("triangle.jl")
include("rectangle.jl")
include("hexagon.jl")
include("construct.jl")
include("deadends.jl")
include("plotting.jl")
include("algorithms.jl")
include("utils.jl")

end


text:
include("linesearch/krylov/root.jl")
include("linesearch/root.jl")
include("trustregions/root.jl")
include("spectral/root.jl")

ext: jl


text:
function ab(G, i, x, θ, C::LocalBound, ∇ϕi, vi, Z::ZigZag, args...)
    a = C.c[i] + ∇ϕi'*θ[i]
    b = C.c[i]/100 + vi
    a, b, 2.0/C.c[i]/abs(θ[i])
end
#(idot(Z.Γ, i, x)  - idot(Z.Γ, i, Z.μ))'*θ[i]),  θ[i]'*idot(Z.Γ, i, θ)

"""
spdmp_inner!(rng, G, G1, G2, ∇ϕ, t, x, θ, Q, c, b, t_old, (acc, num),
F::Union{ZigZag,FactBoomerang}, args...; factor=1.5, adapt=false, adaptscale=false)

[Outdated]
Inner loop of the factorised samplers: the factorised Boomerang algorithm and
the Zig-Zag sampler. Given a dependency graph `G`, gradient `∇ϕ`,
current position `x`, velocity `θ`, Queue of events `Q`, time `t`, tuning parameter `c`,
terms of the affine bounds `a`,`b` and time when the upper bounds were computed `t_old`

The sampler 1) extracts from the queue the first event time. 2) moves deterministically
according to its dynamics until event time. 3) Evaluates whether the event
time is a accepted reflection or refreshment time or shadow time. 4) If it is a
reflection time, the velocity reflects according its reflection rule, if it is a
refreshment time, the sampler updates the velocity from its prior distribution (Gaussian).
In both cases, updates `Q` according to the dependency graph `G`. The sampler proceeds
until the next accepted reflection time or refreshment time. `(num, acc)`
incrementally counts how many event times occour and how many of those are real reflection times.
"""
function spdmp_inner!(rng, G, G2, ∇ϕ, t, x, θ, Q, c::LocalBound, (b, renew), t_old, (acc, num),
 F::Union{ZigZag,FactBoomerang,JointFlow}, args...; factor=1.5, adapt=false, adaptscale=false)
    n = length(x)
    while true
        i, t′ = peek(Q)
     #   @show i
        refresh = i > n
        if refresh
            i = rand(1:n)
        end
        if refresh && F isa JointFlow
            t, x, θ = smove_forward!(G, i, t, x, θ, t′, F)
            t, x, θ = smove_forward!(G2, i, t, x, θ, t′, F)
            x[i], θ[i] = ZigZagBoomerang.refresh(rng, x[i], θ[i], F[i])
            Q[(n + 1)] = t′ + waiting_time_ref(rng, F) #renew refreshment time
            for j in neighbours(G, i)
                ∇ϕj, vj = ∇ϕ(t, x, θ, j, t′, F, args...)
                b[j] = ab(G, j, x, θ, c, ∇ϕj, vj, F)
                t_old[j] = t[j]
                τ, renew[j] = next_time(t[j], b[j], rand(rng))
                Q[j] = τ
            end
         #   @show x, θ, i, t′
        elseif renew[i]
            t, x, θ = smove_forward!(G, i, t, x, θ, t′, F)
            ∇ϕi, vi = ∇ϕ(t, x, θ, i, t′, F, args...)
            b[i] = ab(G, i, x, θ, c, ∇ϕi, vi, F)
            t_old[i] = t[i]
            τ, renew[i] = next_time(t[i], b[i], rand(rng))
            Q[i] = τ
            continue
        else
            t, x, θ = smove_forward!(G, i, t, x, θ, t′, F)
            ∇ϕi, vi = ∇ϕ(t, x, θ, i, t′, F, args...)
            l, lb = sλ(∇ϕi, i, x, θ, F), pos(b[i][1] + b[i][2]*(t[i] - t_old[i]))
            #@show i, θ, l, lb
            num += 1
            if rand(rng)*lb < l
                acc[i] += 1
                if l >= lb
                    if !adapt
                        @show t[i] - t_old[i]
                        @show ∇ϕi, vi, l, lb, b
                        error("Tuning parameter `c = $(c.c[i])` too small.")
                    end
                    c.c[i] = c.c[i]*factor
                end
                t, x, θ = smove_forward!(G2, i, t, x, θ, t′, F)
                θ = reflect!(i, ∇ϕi, x, θ, F)
                #∇ϕi, vi = ∇ϕ(t, x, θ, i, t′, F, args...)
                for j in neighbours(G, i)
                    ∇ϕj, vj = ∇ϕ(t, x, θ, j, t′, F, args...)
                    b[j] = ab(G, j, x, θ, c, ∇ϕj, vj, F)
                    t_old[j] = t[j]
                    τ, renew[j] = next_time(t[j], b[j], rand(rng))
                    Q[j] = τ
                end
            else
                b[i] = ab(G, i, x, θ, c, ∇ϕi, vi, F)
                t_old[i] = t[i]
                τ, renew[i] = next_time(t[i], b[i], rand(rng))
                Q[i] = τ
                continue
            end
        end
        return event(i, t, x, θ, F), t, x, θ, t′, (acc, num), c, (b, renew), t_old
    end
end

"""
spdmp(∇ϕ, t0, x0, θ0, T, c, [G,] F::Union{ZigZag,FactBoomerang}, args...;
    factor=1.5, adapt=false)
    = Ξ, (t, x, θ), (acc, num), c

Version of spdmp which assumes that `i` only depends on coordinates
`x[j] for j in neighbours(G, i)`.

It returns a `FactTrace` (see [`Trace`](@ref)) object `Ξ`, which can be collected
into pairs `t => x` of times and locations and discretized with `discretize`.
Also returns the `num`ber of total and `acc`epted Poisson events and updated bounds
`c` (in case of `adapt==true` the bounds are multiplied by `factor` if they turn
out to be too small.) The final time, location and momentum at `T` can be obtained
with `smove_forward!(t, x, θ, T, F)`.
"""
function spdmp(∇ϕ, t0, x0, θ0, T, C::LocalBound, G_, F::Union{ZigZag,FactBoomerang,JointFlow}, args...;
    factor=1.8, adapt=false, adaptscale=false, progress=false, progress_stops = 20, seed=Seed())
    n = length(x0)

    rng = Rng(seed)
    t′ = t0
    t = fill(t′, size(θ0)...)
    t_old = copy(t)
    if G_ == All()
        G = [i=>collect(1:n) for i in 1:n]
        G2 = nothing
    else
        G = G_
        G2 = [i => setdiff(union((G[j].second for j in G[i].second)...), G[i].second) for i in eachindex(G)]
    end
    x, θ = copy(x0), copy(θ0)
    num = 0
    acc = zeros(Int, length(θ))
    Q = SPriorityQueue{Int,Float64}()
    ∇ϕi, vi = ∇ϕ(t, x, θ, 1, t′, F, args...)
    b = fill(ab(G, 1, x, θ, C, ∇ϕi, vi, F), n)
    renew = zeros(Bool, n)
    for i in eachindex(θ)
        ∇ϕi, vi = ∇ϕ(t, x, θ, i, t′, F, args...)
        b[i] = ab(G, i, x, θ, C, ∇ϕi, vi, F)
        τ, renew[i] = next_time(t[i], b[i], rand(rng))
        enqueue!(Q, i => τ)
    end
    if hasrefresh(F)
        enqueue!(Q, (n + 1) => waiting_time_ref(F))
    end
    if progress
        prg = Progress(progress_stops, 1)
    else
        prg = missing
    end
    stops = ismissing(prg) ? 0 : max(prg.n - 1, 0) # allow one stop for cleanup
    tstop = T/stops
    Ξ = Trace(t0, x0, θ0, F)
    while t′ < T
        ev, t, x, θ, t′, (acc, num), C, (b, renew), t_old = spdmp_inner!(rng, G, G2, ∇ϕ, t, x, θ, Q,
                    C, (b, renew), t_old, (acc, num), F, args...; factor=factor, adapt=adapt, adaptscale=adaptscale)
        push!(Ξ, ev)
        if t′ > tstop
            tstop += T/stops
            next!(prg)
        end
    end
    ismissing(prg) || ProgressMeter.finish!(prg)
    #t, x, θ = smove_forward!(t, x, θ, T, F)
    Ξ, (t, x, θ), (acc, num), C
end

spdmp(∇ϕ, t0, x0, θ0, T, C::LocalBound, F::Union{ZigZag,FactBoomerang,JointFlow}, args...; kargs...) = spdmp(∇ϕ, t0, x0, θ0, T, C, [i => rowvals(F.Γ)[nzrange(F.Γ, i)] for i in eachindex(θ0)], F, args...; kargs...)
pdmp(∇ϕ, t0, x0, θ0, T, C::LocalBound, F::Union{ZigZag,FactBoomerang,JointFlow}, args...; kargs...) = spdmp(∇ϕ, t0, x0, θ0, T, C, All(), F, args...; kargs...)


text:
module AutoHOOT

include("backend.jl")
include("autodiff.jl")
include("graph_ops.jl")

end


text:
"""
    HTTPError <: Exception

Wraps a `Requests.Response` when the response status is >= 400.
"""
immutable HTTPError <: Exception
    resp::Response
end

response(err::HTTPError) = err.resp

function Base.showerror(io::IO, err::HTTPError)
    msg = readstring(err.resp)
    print(io, "HTTPError: $msg")
end

"""
    EtcdError <: Exception

Wraps an etcd response `Dict` when the response contains an `"errorCode"`.
"""
immutable EtcdError <: Exception
    resp::Dict
end

response(err::EtcdError) = err.resp

function Base.showerror(io::IO, err::EtcdError)
    err_code = err.resp["errorCode"]
    msg = get(err.resp, "message", "Unknown error")
    cause = get(err.resp, "cause", "")
    isempty(cause) || (cause = " - $cause")
    print(io, "EtcdError: $(msg)$(cause) ($err_code).")
end

"""
    request(f, uri, opts; n=5, max_delay=10.0)

Executes the HTTP request `f` with the given `uri` and `opts` and returns the response body.
JSON responses are parsed with `Requests.json` otherwise a `String` is returned.
For failed HTTP or etcd requests an `HTTPError` or `EtcdError` is thrown respectively.

# Arguments
- `f::Function`: a valid `Requests` method (e.g., `Requests.get`)
- `uri::String`: the HTTP uri
- `opts::Dict`: a set of options to be passed to `f` with the keyword `query`
- `n=5`: the number of retry attempts for the request
- `max_delay=10.0`: delay in seconds between each retry attempt

# Throws
- `HTTPError`: when the HTTP response has a status code >= 400
- `EtcdError`: when the HTTP response body contains "errorCode" which signifies an etcd error

# Returns
- `Dict`: for most etcd response (the parsed json HTTP response body)
- `Array`: some etcd responses return a list vs dict when parse the json HTTP response body
- `String`: all non-json response bodies
"""
function request(f::Function, uri::String, opts::Dict; n=5, max_delay=10.0)
    logger = getlogger(current_module())

    retry_cond(state, resp) = (state, isa(resp, Response) && in(statuscode(resp), 300:400) && haskey(resp.headers, "Location"))
    retry_func() = isempty(opts) ? f(uri) : f(uri; query=opts)

    resp = retry(
            retry_func;
            delays=Base.ExponentialBackOff(n=n, max_delay=max_delay),
            check=retry_cond
        )()

    debug(logger, readstring(resp))

    data = try
        Requests.json(resp)
    catch _
        readstring(resp)
    end

    if isa(data, Dict) && Base.haskey(data, "errorCode")
        error(logger, EtcdError(data))
    elseif statuscode(resp) >= 400
        error(logger, HTTPError(resp))
    end

    return data
end


text:
module TestDefComposite

using Test
using Mimi
using MacroTools

import Mimi: ComponentPath, build, @defmodel

@defcomp Comp1 begin
    par_1_1 = Parameter(index=[time])      # external input
    var_1_1 = Variable(index=[time])       # computed
    foo = Parameter()

    function run_timestep(p, v, d, t)
        v.var_1_1[t] = p.par_1_1[t]
    end
end

@defcomp Comp2 begin
    par_2_1 = Parameter(index=[time])      # connected to Comp1.var_1_1
    par_2_2 = Parameter(index=[time])      # external input
    var_2_1 = Variable(index=[time])       # computed
    foo = Parameter()

    function run_timestep(p, v, d, t)
        v.var_2_1[t] = p.par_2_1[t] + p.foo * p.par_2_2[t]
    end
end

@defcomposite A begin
    component(Comp1)
    component(Comp2)

    # imports
    bar = Comp1.par_1_1
    foo2 = Comp2.foo

    # linked imports
    # foo = Comp1.foo, Comp2.foo

    foo1 = Comp1.foo
    foo2 = Comp2.foo

    # connections
    Comp2.par_2_1 = Comp1.var_1_1
    Comp2.par_2_2 = Comp1.var_1_1
end

# doesn't work currently
# @defmodel m begin
#     index[time] = 2005:2020
#     component(A)

#     A.foo1 = 10
#     A.foo2 = 4
# end

m = Model()
years = 2005:2020
set_dimension!(m, :time, years)
add_comp!(m, A)

set_param!(m, "/A/Comp1", :par_1_1, 2:2:2*length(years))

a = m.md[:A]
set_param!(a, :Comp1, :foo, 10)
set_param!(a, :Comp2, :foo, 4)      # TBD: why does this overwrite the 10 above??

build(m)
run(m)

end # module

m = TestDefComposite.m
A = TestDefComposite.A
md = m.md

nothing


text:
module Parser

using Ferrite, SparseArrays
using .....TopOpt: find_black_and_white

export extract_inp, InpContent

"""
```
struct InpContent{dim, TF, N, TI}
    node_coords::Vector{NTuple{dim,TF}}
    celltype::String
    cells::Vector{NTuple{N,TI}}
    nodesets::Dict{String,Vector{TI}}
    cellsets::Dict{String,Vector{TI}}
    E::TF
    ν::TF
    density::TF
    nodedbcs::Dict{String, Vector{Tuple{TI,TF}}}
    cloads::Dict{Int, Vector{TF}}
    facesets::Dict{String, Vector{Tuple{TI,TI}}}
    dloads::Dict{String, TF}
end
```

- `node_coords`: a vector of node coordinates.
- `celltype`: a cell type code in the INP convention
- `cells`: a vector of cell connectivities
- `nodesets`: a dictionary mapping a node set name to a vector of node indices
- `cellsets`: a dictionary mapping a cell set name to a vector of cell indices
- `E`: Young's modulus
- `ν`: Poisson ratio
- `density`: physical density of the material
- `nodedbcs`: a dictionary mapping a node set name to a vector of tuples of type `Tuple{Int, Float64}` specifying a Dirichlet boundary condition on that node set. Each tuple in the vector specifies the local index of a constrained degree of freedom and its fixed value. A 3-dimensional field has 3 degrees of freedom per node for example. So the index can be 1, 2 or 3.
- `cloads`: a dictionary mapping a node index to a load vector on that node.
- `facesets`: a dictionary mapping a face set name to a vector of `Tuple{Int,Int}` tuples where each tuple is a face index. The first integer is the cell index where the face is and the second integer is the local face index in the cell according to the VTK convention.
- `dloads`: a dictionary of distributed loads mapping face set names to a normal traction load value.
"""
struct InpContent{dim, TF, N, TI}
    node_coords::Vector{NTuple{dim,TF}}
    celltype::String
    cells::Vector{NTuple{N,TI}}
    nodesets::Dict{String,Vector{TI}}
    cellsets::Dict{String,Vector{TI}}
    E::TF
    ν::TF
    density::TF
    nodedbcs::Dict{String, Vector{Tuple{TI,TF}}}
    cloads::Dict{Int, Vector{TF}}
    facesets::Dict{String, Vector{Tuple{TI,TI}}}
    dloads::Dict{String, TF}
end

const stopping_pattern = r"^\*[^\*]"

# Import
include(joinpath("FeatureExtractors", "extract_cells.jl"))
include(joinpath("FeatureExtractors", "extract_cload.jl"))
include(joinpath("FeatureExtractors", "extract_dbcs.jl"))
include(joinpath("FeatureExtractors", "extract_dload.jl"))
include(joinpath("FeatureExtractors", "extract_material.jl"))
include(joinpath("FeatureExtractors", "extract_nodes.jl"))
include(joinpath("FeatureExtractors", "extract_set.jl"))
include(joinpath("inp_to_ferrite.jl"))

function extract_inp(filepath_with_ext)
    file = open(filepath_with_ext, "r")

    local node_coords
    local celltype, cells, offset
    nodesets = Dict{String,Vector{Int}}()
    cellsets = Dict{String,Vector{Int}}()
    local E, mu
    nodedbcs = Dict{String, Vector{Tuple{Int,Float64}}}()
    cloads = Dict{Int, Vector{Float64}}()
    facesets = Dict{String, Vector{Tuple{Int,Int}}}()
    dloads = Dict{String, Float64}()
    density = 0. # Should extract from the file

    node_heading_pattern = r"\*Node\s*,\s*NSET\s*=\s*([^,]*)"
    cell_heading_pattern = r"\*Element\s*,\s*TYPE\s*=\s*([^,]*)\s*,\s*ELSET\s*=\s*([^,]*)"
    nodeset_heading_pattern = r"\*NSET\s*,\s*NSET\s*=\s*([^,]*)"
    cellset_heading_pattern = r"\*ELSET\s*,\s*ELSET\s*=\s*([^,]*)"
    material_heading_pattern = r"\*MATERIAL\s*,\s*NAME\s*=\s*([^\s]*)"
    boundary_heading_pattern = r"\*BOUNDARY"
    cload_heading_pattern = r"\*CLOAD"
    dload_heading_pattern = r"\*DLOAD"

    line = readline(file)
    local dim
    while !eof(file)
        m = match(node_heading_pattern, line)
        if m != nothing && m[1] == "Nall"
            node_coords, line = extract_nodes(file)
            dim = length(node_coords[1])
            continue
        end
        m = match(cell_heading_pattern, line)
        if m != nothing
            celltype = String(m[1])
            cellsetname = String(m[2])
            cells, offset, line = extract_cells(file)
            cellsets[cellsetname] = collect(1:length(cells))
            continue
        end
        m = match(nodeset_heading_pattern, line)
        if m != nothing
            nodesetname = String(m[1])
            line = extract_set!(nodesets, nodesetname, file)
            continue
        end
        m = match(cellset_heading_pattern, line)
        if m != nothing
            cellsetname = String(m[1])
            line = extract_set!(cellsets, cellsetname, file, offset)
            continue
        end
        m = match(material_heading_pattern, line)
        if m != nothing
            material_name = String(m[1])
            E, mu, line = extract_material(file)
            continue
        end
        m = match(boundary_heading_pattern, line)
        if m != nothing
            line = extract_nodedbcs!(nodedbcs, file)
            continue
        end
        m = match(cload_heading_pattern, line)
        if m != nothing
            line = extract_cload!(cloads, file, Val{dim})
            continue
        end
        m = match(dload_heading_pattern, line)
        if m != nothing
            line = extract_dload!(dloads, facesets, file, Val{dim}, offset)
            continue
        end
        line = readline(file)
    end

    close(file)

    return InpContent(node_coords, celltype, cells, nodesets, cellsets, E, mu, density, nodedbcs, cloads, facesets, dloads)
end

end


text:
using AWSCore: AWSConfig, AWSException
using Base: CmdRedirect
using Base64
using Dates: datetime2unix, now
const invalid_access_key = "ThisIsMyInvalidAccessKey"
const invalid_secret_key = "ThisIsMyInvalidSecretKey"

get_caller_identity_patch = @patch function AWSTools.get_caller_identity()
    account_id = join(rand(0:9, 12), "")
    Dict(
        "Account" => account_id,
        "Arn" => "arn:aws:iam::$account_id:user/UserName",
        "UserId" => join(rand('A':'Z', 21), ""),
    )
end

sts_assume_role = @patch function AWSTools.sts(config::AWSConfig, op; kwargs...)
    Dict(
        "Credentials" => Dict(
            "AccessKeyId" => "TESTACCESSKEYID",
            "SecretAccessKey" => "TESTSECRETACEESSKEY",
            "SessionToken" => "TestSessionToken",
            "Expiration" => datetime2unix(now()),
        ),
    )
end

function instance_metadata_patch(result)
    @patch HTTP.get(args...; kwargs...) = result
end

get_authorization_token_patch = @patch function AWSTools.ECR.get_authorization_token(
    config::AWSConfig;
    registryIds::AbstractVector=[],
)
    id = lpad(isempty(registryIds) ? "" : first(registryIds), 12, '0')
    Dict(
        "authorizationData" => [
            Dict(
                "authorizationToken" => base64encode("AWS:password"),
                "proxyEndpoint" => "https://$(id).dkr.ecr.us-east-1.amazonaws.com"
            ),
        ]
    )
end

describe_stacks_patch = @patch function AWSTools.CloudFormation.describe_stacks(
    config,
    args...;
    kwargs...,
)
    responses = Dict(
       Dict(:StackName => "stackname") =>
        """
        <DescribeStacksResponse xmlns="http://cloudformation.amazonaws.com/doc/2010-05-15/">
          <DescribeStacksResult>
              <Stacks>
                <member>
                  <StackId>Stack Id</StackId>
                  <StackName>Stack Name</StackName>
                  <Description>Stack Description</Description>
                </member>
              </Stacks>
          </DescribeStacksResult>
        </DescribeStacksResponse>
        """,
        Dict(:StackName => "1-stack-output-stackname") =>
        """
        <DescribeStacksResponse xmlns="http://cloudformation.amazonaws.com/doc/2010-05-15/">
          <DescribeStacksResult>
            <Stacks>
              <member>
                <Outputs>
                  <member>
                    <OutputKey>TestBucketArn1</OutputKey>
                    <OutputValue>arn:aws:s3:::test-bucket-1</OutputValue>
                  </member>
                </Outputs>
                <StackId>Stack Id</StackId>
                <StackName>Stack Name</StackName>
                <Description>Stack Description</Description>
              </member>
            </Stacks>
          </DescribeStacksResult>
        </DescribeStacksResponse>
        """,
        Dict(:StackName => "multiple-stack-outputs-stackname") =>
        """
        <DescribeStacksResponse xmlns="http://cloudformation.amazonaws.com/doc/2010-05-15/">
          <DescribeStacksResult>
            <Stacks>
              <member>
                <Outputs>
                  <member>
                    <OutputKey>TestBucketArn1</OutputKey>
                    <OutputValue>arn:aws:s3:::test-bucket-1</OutputValue>
                  </member>
                  <member>
                    <OutputKey>TestBucketArn2</OutputKey>
                    <OutputValue>arn:aws:s3:::test-bucket-2</OutputValue>
                  </member>
                </Outputs>
                <StackId>Stack Id</StackId>
                <StackName>Stack Name</StackName>
                <Description>Stack Description</Description>
              </member>
            </Stacks>
          </DescribeStacksResult>
        </DescribeStacksResponse>
        """,
        Dict(:StackName => "empty-value") =>
        """
        <DescribeStacksResponse xmlns="http://cloudformation.amazonaws.com/doc/2010-05-15/">
          <DescribeStacksResult>
            <Stacks>
              <member>
                <Outputs>
                  <member>
                    <OutputKey>ParquetConversionTriggerName</OutputKey>
                    <OutputValue></OutputValue>
                  </member>
                </Outputs>
                <StackId>Stack Id</StackId>
                <StackName>Stack Name</StackName>
                <Description>Stack Description</Description>
              </member>
            </Stacks>
          </DescribeStacksResult>
        </DescribeStacksResponse>
        """,
        Dict(:StackName => "export") =>
        """
        <DescribeStacksResponse xmlns="http://cloudformation.amazonaws.com/doc/2010-05-15/">
          <DescribeStacksResult>
            <Stacks>
              <member>
                <Outputs>
                  <member>
                    <Description>Exported output for use in other stacks</Description>
                    <ExportName>ExportedKey</ExportName>
                    <OutputKey>Key</OutputKey>
                    <OutputValue>Value</OutputValue>
                  </member>
                </Outputs>
                <StackId>Stack Id</StackId>
                <StackName>Stack Name</StackName>
                <Description>Stack Description</Description>
              </member>
            </Stacks>
          </DescribeStacksResult>
        </DescribeStacksResponse>
        """,
    )

    # So we can test that we get an error using the invalid access and secret keys
    access_key = config[:creds].access_key_id
    secret_key = config[:creds].secret_key
    if access_key == invalid_access_key && secret_key == invalid_secret_key
        throw(AWSException(
            HTTP.StatusError(403, "", "", HTTP.Messages.Response(403, """
                <ErrorResponse xmlns="http://cloudformation.amazonaws.com/doc/2010-05-15/">
                    <Error>
                    <Type>Sender</Type>
                    InvalidClientTokenId
                    <Message>The security token included in the request is invalid.</Message>
                    </Error>
                    <RequestId>cff5beb8-4b7b-11e9-9c2b-43c18f6078dc</RequestId>
                </ErrorResponse>
                """
            ))
        ))
    else
        return responses[Dict{Symbol, String}(kwargs)]
    end
end

function throttle_patch(allow)
    describe_stacks_throttle_count = 0
    @patch function AWSTools.CloudFormation.describe_stacks(args...; kwargs...)
        describe_stacks_throttle_count += 1
        if !(describe_stacks_throttle_count in allow)
            error_message = """
                <ErrorResponse xmlns="http://cloudformation.amazonaws.com/doc/2010-05-15/">
                    <Error>
                        <Type>Sender</Type>
                        Throttling
                        <Message>Rate exceeded</Message>
                    </Error>
                    <RequestId>d0c477ac-f267-11e8-9d2b-93e3aa6368c5</RequestId>
                </ErrorResponse>
                """
            response = HTTP.Messages.Response(400, error_message)
            http_error = HTTP.ExceptionRequest.StatusError(400, "", "", response)
            throw(AWSException(http_error))
        end
        responses = Dict(
           Dict(:StackName => "stackname") =>
            """
            <DescribeStacksResponse xmlns="http://cloudformation.amazonaws.com/doc/2010-05-15/">
              <DescribeStacksResult>
                  <Stacks>
                    <member>
                      <StackId>Stack Id</StackId>
                      <StackName>Stack Name</StackName>
                      <Description>Stack Description</Description>
                      <ThrottleCount>$describe_stacks_throttle_count</ThrottleCount>
                    </member>
                  </Stacks>
              </DescribeStacksResult>
            </DescribeStacksResponse>
            """,
        )

        return responses[Dict{Symbol, String}(kwargs)]
    end
end

function format_s3_objects(content::AbstractDict)
    objects = OrderedDict{Tuple{String,String}, Dict{String,Any}}()
    for ((bucket::String, key::String), data) in content
        formatted_data = Dict{String, Any}(k => v for (k, v) in data)

        # Add in defaults for required keys
        get!(formatted_data, "Key", key)
        get!(formatted_data, "LastModified", "1970-01-01T00:00:00.000Z")
        get!(formatted_data, "Size") do
            len = if haskey(formatted_data, "Content")
                sizeof(formatted_data["Content"])
            else
                0
            end
            string(len)
        end

        push!(objects, (bucket, key) => formatted_data)
    end

    return objects
end

function format_s3_objects(content::AbstractVector{Pair{Tuple{String,String},D}}) where D <: AbstractDict
    format_s3_objects(OrderedDict(k => v for (k, v) in content))
end

function format_s3_objects(content::AbstractVector{Tuple{String,String}})
    format_s3_objects(OrderedDict(k => Dict() for k in content))
end

function s3_patches!(content::AbstractDict, changes::Set=Set(Pair{Symbol, Any}[]))
    return [
        @patch function AWSS3.s3_list_objects(
            config::AWSConfig,
            bucket,
            path_prefix;
            kwargs...,
        )
            results = []
            for ((b, key), data) in content
                if b == bucket && startswith(key, path_prefix)
                    push!(results, data)
                end
            end
            return results
        end

        # Patch calling another patch
        @patch s3_list_objects(args...) = @mock s3_list_objects(aws_config(), args...)

        @patch function AWSS3.s3_get(config::AWSConfig, bucket, path; kwargs...)
            codeunits(get(content[(bucket, path)], "Content", ""))
        end

        @patch function AWSS3.s3_exists(config::AWSConfig, bucket, path)
            haskey(content, (bucket, path))
        end

        @patch function AWSTools.S3.s3_copy(
            config::AWSConfig,
            bucket,
            path;
            to_bucket=bucket,
            to_path=path,
            kwargs...,
        )
            push!(
                changes,
                :copy => Dict(
                    :from => (bucket, path),
                    :to => (to_bucket, to_path),
                )
            )
        end

        @patch function AWSS3.s3_delete(config::AWSConfig, bucket, path; kwargs...)
            push!(changes, :delete => (bucket, path))
        end

        @patch function AWSS3.s3_put(
            config::AWSConfig,
            bucket,
            path,
            data::Union{AbstractString, Vector{UInt8}};
            kwargs...,
        )
            push!(
                changes,
                :put => Dict(
                    :dest => (bucket, path),
                    :data => data,
                )
            )
        end
    ]
end


text:
#=
In this file we define how agents are plotted and how the plots are updated while stepping
=#
# TODO: I should check whether it is worth to type-parameterize this.
struct ABMStepper # {X, C, M, S, O, AC, AS, AM, HA}
    ac # ::C
    am # ::M
    as # ::S
    offset # ::O
    scheduler # ::X
    pos # ::Observable
    colors # ::AC
    sizes # ::AS
    markers # ::AM
    heatarray # ::HA
    heatobs # ::HO
end

Base.show(io::IO, ::ABMStepper) =
println(io, "Helper structure for stepping and updating the plot of an agent based model. ",
"It is outputted by `abm_plot` and can be used in `Agents.step!`, see `abm_plot`.")

"Initialize the abmstepper and the plotted observables. Return the stepper."
function abm_init_stepper(model; ac, am, as, scheduler, offset, heatarray)

    if !isnothing(heatarray)
        # TODO: This is also possible for continuous spaces, we have to
        # get the matrix size, and then make a range for each dimension
        # and do heatmap!(ax, x, y, heatobs)
        #
        # TODO: use surface!(heatobs) here?
        matrix = Agents.get_data(model, heatarray, identity)
        if !(matrix isa AbstractMatrix) || size(matrix) ≠ size(model.space)
            error("The heat array property must yield a matrix of same size as the grid!")
        end
        heatobs = Observable(matrix)
    else
        heatobs = nothing
    end

    ids = scheduler(model)
    colors = Observable(ac isa Function ? to_color.([ac(model[i]) for i ∈ ids]) : to_color(ac))
    sizes = Observable(as isa Function ? [as(model[i]) for i ∈ ids] : as)
    markers = Observable(am isa Function ? [am(model[i]) for i ∈ ids] : am)
    pos = Observable(agents_pos_for_plotting(model, offset, ids))
    if user_used_polygons(am, markers)
        if am isa Function
            markers[] = [translate(m, p) for (m, p) in zip(markers[], pos[])]
        else
            markers[] = [translate(am, p) for p in pos[]]
        end
        # For polygons we always need vector, even if all agents are same polygon
    end

    return ABMStepper(
        ac, am, as, offset, scheduler,
        pos, colors, sizes, markers,
        heatarray, heatobs
    )
end

function agents_pos_for_plotting(model, offset, ids)
    if model.space isa Agents.OpenStreetMapSpace
        if isnothing(offset)
            pos = [Point2f0(Agents.OSM.lonlat(model[i].pos, model)) for i in ids]
        else
            pos = [Point2f0(Agents.OSM.lonlat(model[i].pos, model) .+ offset(model[i])) for i ∈ ids]
        end
        return pos
    end
    # standard space case
    postype = agents_space_dimensionality(model.space) == 3 ? Point3f0 : Point2f0
    if isnothing(offset)
        pos = [postype(model[i].pos) for i ∈ ids]
    else
        pos = [postype(model[i].pos .+ offset(model[i])) for i ∈ ids]
    end
    return pos
end

agents_pos_for_plotting(abms::ABMStepper, model, ids = abms.scheduler(model)) =
agents_pos_for_plotting(model, abms.offset, ids)

agents_space_dimensionality(abm::Agents.ABM) = agents_space_dimensionality(abm.space)
agents_space_dimensionality(::Agents.GridSpace{D}) where {D} = D
agents_space_dimensionality(::Agents.ContinuousSpace{D}) where {D} = D
agents_space_dimensionality(::Agents.OpenStreetMapSpace) = 2

SUPPORTED_AGENTS_SPACES =  Union{
    Agents.DiscreteSpace,
    Agents.ContinuousSpace,
    Agents.OpenStreetMapSpace,
}

"Initialize the ABM plot and return it."
function abm_init_plot!(ax, model, abmstepper::ABMStepper;
        heatkwargs, add_colorbar, static_preplot!, scatterkwargs
    )

    @assert typeof(model.space) <: SUPPORTED_AGENTS_SPACES
    plot_agents_space!(ax, model)

    if !isnothing(abmstepper.heatobs)
        heatkwargs = merge((colormap=JULIADYNAMICS_CMAP,), heatkwargs)
        hmap = heatmap!(ax, abmstepper.heatobs[]; heatkwargs...)
        if add_colorbar
            fig = ax.parent
            Colorbar(fig[1, 1][1, 2], hmap, width = 20)
        end
        # rowsize!(fig[1,1].fig.layout, 1, ax.scene.px_area[].widths[2]) # Colorbar height = axis height
    end

    static_plot = static_preplot!(ax, model)
    !isnothing(static_plot) && (static_plot.inspectable[] = false)

    # Here we make the decision of whether the user has provided markers, and thus use
    # `scatter`, or polygons, and thus use `poly`:
    if user_used_polygons(abmstepper.am, abmstepper.markers)
        poly!(ax, abmstepper.markers; color = abmstepper.colors, scatterkwargs...)
    else
        scatter!(ax, abmstepper.pos;
            color = abmstepper.colors, marker = abmstepper.markers,
            markersize = abmstepper.sizes, scatterkwargs...
        )
    end
    return
end

"Plot space and/or set axis limits."
function plot_agents_space!(ax, model)
    if model.space isa Agents.OpenStreetMapSpace
        Main.OSMMakie.osmplot!(ax, model.space.map)
        return
    end
    if model.space isa Agents.ContinuousSpace
        e = model.space.extent
    elseif model.space isa Agents.DiscreteSpace
        e = size(model.space.s) .+ 1
    end
    o = zero.(e)
    xlims!(ax, o[1], e[1])
    ylims!(ax, o[2], e[2])
    is3d = length(o) == 3
    is3d && zlims!(ax, o[3], e[3])
    return
end

default_static_preplot(ax, model) = nothing

function user_used_polygons(am, markers)
    if (am isa Polygon)
        return true
    elseif (am isa Function) && (markers[][1] isa Polygon)
        return true
    else
        return false
    end
end

#=
    Agents.step!(abmstepper, model, agent_step!, model_step!, n::Int)
Step the given `model` for `n` steps while also updating the plot that corresponds to it,
which is produced with the function [`abm_plot`](@ref).

You can still call this function with `n=0` to update the plot for a new `model`,
without doing any stepping.
=#
function Agents.step!(abmstepper::ABMStepper, model, agent_step!, model_step!, n)
    @assert (n isa Int) "Only stepping with integer `n` is possible with `abmstepper`."
    ac, am, as = abmstepper.ac, abmstepper.am, abmstepper.as
    pos, colors = abmstepper.pos, abmstepper.colors
    sizes, markers =  abmstepper.sizes, abmstepper.markers

    Agents.step!(model, agent_step!, model_step!, n)

    if Agents.nagents(model) == 0
        @warn "The model has no agents"
    end
    ids = abmstepper.scheduler(model)
    pos[] = agents_pos_for_plotting(abmstepper, model, ids)
    if ac isa Function; colors[] = to_color.([ac(model[i]) for i in ids]); end
    if as isa Function; sizes[] = [as(model[i]) for i in ids]; end
    if am isa Function; markers[] = [am(model[i]) for i in ids]; end
    # If we use Polygons as markers, do a final update:
    if user_used_polygons(am, markers)
        if am isa Function
            markers[] = [translate(m, p) for (m, p) in zip(markers[], pos[])]
        else
            markers[] = [translate(am, p) for p in pos[]]
        end
    end
    # Finally update the heat array, if any
    if !isnothing(abmstepper.heatarray)
        newmatrix = Agents.get_data(model, abmstepper.heatarray, identity)
        abmstepper.heatobs[] = newmatrix
    end
    return nothing
end


text:
using ..VectorEngine: VERefValue

struct VEArgs
    handle::VEDAargs
    objs::Base.IdSet
    function VEArgs()
        r_handle = Ref{VEDAargs}()
        vedaArgsCreate(r_handle)
        return new(r_handle[], Base.IdSet())
    end
end

function Base.setindex!(args::VEArgs, val::Float32, idx::Integer)
    vedaArgsSetF32(args.handle, idx, val)
    val
end

function Base.setindex!(args::VEArgs, val::Float64, idx::Integer)
    @debug("setindex Float64 $idx")
    vedaArgsSetF64(args.handle, idx, val)
    val
end

function Base.setindex!(args::VEArgs, val::Int8, idx::Integer)
    vedaArgsSetI8(args.handle, idx, val)
    val
end

function Base.setindex!(args::VEArgs, val::Int16, idx::Integer)
    vedaArgsSetI16(args.handle, idx, val)
    val
end

function Base.setindex!(args::VEArgs, val::Int32, idx::Integer)
    vedaArgsSetI32(args.handle, idx, val)
    val
end

function Base.setindex!(args::VEArgs, val::Int64, idx::Integer)
    @debug("setindex Int64 $idx")
    vedaArgsSetI64(args.handle, idx, val)
    val
end

function Base.setindex!(args::VEArgs, val::UInt8, idx::Integer)
    vedaArgsSetU8(args.handle, idx, val)
    val
end

function Base.setindex!(args::VEArgs, val::UInt16, idx::Integer)
    vedaArgsSetU16(args.handle, idx, val)
    val
end

function Base.setindex!(args::VEArgs, val::UInt32, idx::Integer)
    vedaArgsSetU32(args.handle, idx, val)
    val
end

function Base.setindex!(args::VEArgs, val::UInt64, idx::Integer)
    vedaArgsSetU64(args.handle, idx, val)
    val
end

# Pass string variables on stack. Expect Ptr{...} on the VE side.
# - immutable, therefore only INTENT_IN
function Base.setindex!(args::VEArgs, val::T, idx::Integer) where {T<:AbstractString}
    push!(args.objs, val)
    vedaArgsSetStack(args.handle, idx, Base.unsafe_convert(Ptr{Cvoid}, pointer(val)),
                         VEDA_ARGS_INTENT_IN, sizeof(val))
    val
end

function equivalent_uint(::Type{T}) where T
    if sizeof(T) == 1
        return UInt8
    elseif sizeof(T) == 2
        return UInt16
    elseif sizeof(T) == 4
        return UInt32
    elseif sizeof(T) == 8
        return UInt64
    else
        return nothing
    end
end

# Pass structs on stack. If mutable, modified elements can be passed back, too.
function Base.setindex!(args::VEArgs, val::T, idx::Integer) where {T}
    @debug("args setindex! idx=$idx T=$T val=$val")
    if isstructtype(T) && Base.datatype_pointerfree(T) # Pointers within structs can not be de-referenced
        @debug("arg $idx val=$val handled as struct")
        intent = ismutable(val) ? VEDA_ARGS_INTENT_INOUT : VEDA_ARGS_INTENT_IN
        ref = Ref(val)
        push!(args.objs, ref) # root the reference
        ptr = Base.unsafe_convert(Ptr{Cvoid}, ref)
        vedaArgsSetStack(args.handle, idx, ptr, intent, sizeof(T))
        val
    elseif T <: Base.RefValue && isassigned(val)
        @debug("arg $idx val=$val handled as RefValue")
        intent = ismutable(val[]) ? VEDA_ARGS_INTENT_INOUT : VEDA_ARGS_INTENT_IN
        push!(args.objs, val) # root the reference
        ptr = Base.unsafe_convert(Ptr{Cvoid}, val)
        sz = sizeof(val[])
        @debug "converted ptr = $ptr size=$sz"
        vedaArgsSetStack(args.handle, idx, ptr, intent, sizeof(val[]))
        val
    elseif T <: Ptr{Nothing}
        @debug("arg $idx val=$val handled as Ptr{Nothing}")
        args[idx] = Base.bitcast(UInt64, val)
        val
    else
        error("VEArgs could not handle $T")
    end
end

# convert the argument values to match the kernel's signature (specified by the user)
# (this mimics `lower-ccall` in julia-syntax.scm)
@generated function convert_args(veargs, ::Type{tt}, args...) where {tt}
    types = tt.parameters
    @debug("types=$types tt=$tt")

    ex = quote
        Base.@_inline_meta
    end

    # veargs is zero indexed
    for i in 1:length(args)
        push!(ex.args, :(veargs[$(i-1)] = args[$i]))
    end

    return ex
end

function vecall(func::VEFunction, tt, args...; stream = C_NULL)
    veargs = VEArgs()
    #@show args
    convert_args(veargs, tt, args...)

    # TOOD: Can we add a callback that will do args cleanup?
    err = vedaLaunchKernelEx(func.handle, stream, veargs.handle, #=destroyArgs=# true, C_NULL)
    if err != VEDA_SUCCESS
        throw(VEOCommandError("kernel launch failed"))
    end
    err, veargs
end


text:
using MDToolbox
using Test

@testset "trjarray.jl" begin
    @testset "null object" begin
        @test TrjArray().natom == 0
        @test TrjArray().nframe == 0
    end

    @testset "constructor" begin
        nframe = 5
        natom = 5
        ta = TrjArray(x=rand(nframe, natom), y=rand(nframe, natom), z=rand(nframe, natom))

        @test ta.natom == natom
        @test ta.nframe == nframe
    end

    @testset "slicing" begin
        nframe = 5
        natom = 5
        ta = TrjArray(x=rand(nframe, natom), y=rand(nframe, natom), z=rand(nframe, natom))

        @test ta[1, :].x ≈ ta.x[1:1, :]
        @test ta[1:3, :].x ≈ ta.x[1:3, :]
        @test ta[[1, 3, 5], :].x ≈ ta.x[[1, 3, 5], :]
        @test ta[[true, true, true, false, false], :].x ≈ ta.x[[true, true, true, false, false], :]

        @test ta[:, 1].x ≈ ta.x[:, 1:1]
        @test ta[:, 1:3].x ≈ ta.x[:, 1:3]
        @test ta[:, [1, 3, 5]].x ≈ ta.x[:, [1, 3, 5]]
        @test ta[:, [true, true, true, false, false]].x ≈ ta.x[:, [true, true, true, false, false]]
    end

    @testset "slicing with boxsize" begin
        nframe = 5
        natom = 5
        ta = TrjArray(x=rand(nframe, natom), y=rand(nframe, natom), z=rand(nframe, natom), boxsize=rand(nframe, 3))

        @test ta[1, :].boxsize ≈ ta.boxsize[1:1, :]
        @test ta[1:3, :].boxsize ≈ ta.boxsize[1:3, :]
        @test ta[[1, 3, 5], :].boxsize ≈ ta.boxsize[[1, 3, 5], :]
        @test ta[[true, true, true, false, false], :].boxsize ≈ ta.boxsize[[true, true, true, false, false], :]
    end

    @testset "atom selection" begin
        nframe = 5
        natom = 5
        ta = TrjArray(x=rand(nframe, natom), y=rand(nframe, natom), z=rand(nframe, natom), atomid=1:natom, chainname=fill("A", natom))

        @test ta["atomid 1:3"].x ≈ ta.x[:, 1:3]
        @test ta["chainname A"].x ≈ ta.x[:, :]
        @test ta["atomid 3:5 and chainname A"].x ≈ ta.x[:, 3:5]
        @test ta["atomid 3:5 or chainname A"].x ≈ ta.x[:, :]
        @test ta["(atomid 3:5 and chainname A) and atomid 3"].x ≈ ta.x[:, 3]
        @test ta["(atomid 3:5 and chainname A) and (atomid 3 and chainname A)"].x ≈ ta.x[:, 3]
        @test ta["not atomid 1:3"].x ≈ ta.x[:, 4:5]
        @test ta["not (atomid 1:2 or atomid 4:5)"].x ≈ ta.x[:, 3]
    end

    @testset "merging" begin
        nframe = 5
        natom = 5
        topology = TrjArray(atomid=1:natom)
        ta = TrjArray(x=rand(nframe, natom), y=rand(nframe, natom), z=rand(nframe, natom))

        ta_new = [topology; ta]
        ta_new["atomid 1:3"].x ≈ ta.x[:, 1:3]
    end
end

@testset "structure.jl" begin
    @testset "centerofmass" begin
        c = [-1.0 0.0 1.0;
             -1.0 0.0 1.0]
        ta = TrjArray(x=c, y=c, z=c)
        com = centerofmass(ta)
        @test com.x ≈ com.y ≈ com.z ≈ [0.0, 0.0]

        c = [-2.0 0.0 1.0;
             -2.0 0.0 1.0]
        ta = TrjArray(x=c, y=c, z=c, mass=[0.5, 1.0, 1.0])
        com = centerofmass(ta)
        @test com.x ≈ com.y ≈ com.z ≈ [0.0, 0.0]

        c = [rand(Float64)]
        c = reshape(c, 1, 1)
        ta = TrjArray(x=c, y=c, z=c)
        com = centerofmass(ta)
        @test com.x ≈ com.y ≈ com.z ≈ c
    end

    @testset "decenter" begin
        c = [1.0 1.0 1.0;
             2.0 2.0 2.0]
        ta = TrjArray(x=c, y=c, z=c)
        ta2, com = decenter(ta)
        @test ta2.x ≈ ta2.y ≈ ta2.z ≈ [0.0 0.0 0.0; 0.0 0.0 0.0]

        c = [rand(Float64)]
        c = reshape(c, 1, 1)
        ta = TrjArray(x=c, y=c, z=c)
        ta2, com = decenter(ta)
        @test ta2.x ≈ ta2.y ≈ ta2.z ≈ [0.0]
    end

    @testset "superimpose and getrmsd" begin
        A = [-2.803  -15.373   24.556;
             0.893  -16.062   25.147;
             1.368  -12.371   25.885;
             -1.651  -12.153   28.177;
             -0.440  -15.218   30.068;
             2.551  -13.273   31.372;
             0.105  -11.330   33.567]
        x = A[:, 1]
        y = A[:, 2]
        z = A[:, 3]
        ref = TrjArray(x=x', y=y', z=z', mass=1:7)
        A = [-14.739  -18.673   15.040;
             -12.473  -15.810   16.074;
             -14.802  -13.307   14.408;
             -17.782  -14.852   16.171;
             -16.124  -14.617   19.584;
             -15.029  -11.037   18.902;
             -18.577  -10.001   17.996]
        x = A[:, 1]
        y = A[:, 2]
        z = A[:, 3]
        ta = TrjArray(x=x', y=y', z=z', mass=1:7)
        # r, ta_fit = superimpose(ref, ta, isweight=true)
        ta_fit = superimpose(ref, ta, isweight=true)
        #@test r[1] ≈ 0.7450161471
        r2 = getrmsd(ref, ta_fit)
        #@test r2[1] ≈ r[1]
        @test r2[1] ≈ 0.7450161471
    end

    @testset "getdistance" begin
        x = [1.0 0.0]
        y = [0.0 0.0]
        z = [0.0 0.0]
        ta = TrjArray(x=x, y=y, z=z)
        d = getdistance(ta[:, 1], ta[:, 2])
        @test d ≈ [1.0]
    end

    @testset "getangle" begin
        x = [1.0 0.0 0.0]
        y = [0.0 0.0 1.0]
        z = [0.0 0.0 0.0]
        ta = TrjArray(x=x, y=y, z=z)
        a = getangle(ta[:, 1], ta[:, 2], ta[:, 3])
        @test a ≈ [90.0]
    end

    @testset "getdihedral" begin
        x = [-1.0 -1.0 1.0 1.0]
        y = [-1.0  0.0 0.0 1.0]
        z = [ 0.0  0.0 0.0 0.0]
        ta = TrjArray(x=x, y=y, z=z)
        a = getdihedral(ta[:, 1], ta[:, 2], ta[:, 3], ta[:, 4])
        @test a ≈ [180.0]
    end
end

@testset "ksdensity.jl" begin
    @testset "ksdensity 1d" begin
        grid_x=collect(-3:0.1:3)
        f, _ = ksdensity(randn(10^6), grid_x=grid_x);
        sigma = 1.0
        f_true = (1.0/(sqrt(2.0*pi)*sigma)).*exp.(-grid_x.*grid_x./(2.0*sigma^2));
        @test maximum(abs.(f - f_true)) < 0.01
    end

    @testset "ksdensity 2d" begin
        grid_x=collect(-2:0.1:2)
        grid_y=collect(-2:0.1:2)
        f, _, _ = ksdensity(randn(10^6), randn(10^6), grid_x=grid_x, grid_y=grid_y);
        sigma = 1.0
        f_true = zeros(Float64, length(grid_x), length(grid_y))
        for i in 1:length(grid_x)
            for j in 1:length(grid_y)
                f_true[i, j] = (1.0/(sqrt(2.0*pi)*sigma))*exp(-grid_x[i]*grid_x[i]/(2.0*sigma^2))*
                            (1.0/(sqrt(2.0*pi)*sigma))*exp(-grid_y[j]*grid_y[j]/(2.0*sigma^2));
            end
        end
        @test maximum(abs.(f - f_true)) < 0.01
    end

    @testset "ksdensity 3d" begin
        grid_x=collect(-1:0.1:1)
        grid_y=collect(-1:0.1:1)
        grid_z=collect(-1:0.1:1)
        f, _, _, _ = ksdensity(randn(10^6), randn(10^6), randn(10^6), grid_x=grid_x, grid_y=grid_y, grid_z=grid_z);
        sigma = 1.0
        f_true = zeros(Float64, length(grid_x), length(grid_y), length(grid_z))
        for i in 1:length(grid_x)
            for j in 1:length(grid_y)
                for k in 1:length(grid_z)
                    f_true[i, j, k] = (1.0/(sqrt(2.0*pi)*sigma))*exp(-grid_x[i]*grid_x[i]/(2.0*sigma^2))*
                            (1.0/(sqrt(2.0*pi)*sigma))*exp(-grid_y[j]*grid_y[j]/(2.0*sigma^2))*
                            (1.0/(sqrt(2.0*pi)*sigma))*exp(-grid_z[j]*grid_z[j]/(2.0*sigma^2));
                end
            end
        end
        @test maximum(abs.(f - f_true)) < 0.05
    end
end


text:
module FrostHelperFastfallTriggerModule

using ..Ahorn, Maple

@mapdef Trigger "FrostHelper/ForcedFastfall" FrostHelperFastfalllTriggerPlacement(x::Integer, y::Integer, width::Integer=16, height::Integer=16)

const placements = Ahorn.PlacementDict(
	"Forced Fastfall (Frost Helper)" => Ahorn.EntityPlacement(
		FrostHelperFastfalllTriggerPlacement,
		"rectangle"
	)
)

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
using Oceananigans.Operators: Δx, Δy, ΔzC, div_xyᶜᶜᵃ, div_xzᶜᵃᶜ, div_yzᵃᶜᶜ

#####
##### Outer functions for setting velocity on boundary and filling halo beyond boundary.
#####

@inline function fill_west_halo!(u, bc::NFBC, arch, grid, clock, state)
    @launch(device(arch), config=launch_config(grid, :yz), set_east_west_u_velocity!(u, 1, bc, grid, clock, state))
    #@launch(device(arch), config=launch_config(grid, :yz), _fill_west_halo!(u, bc, grid, state.velocities.v, state.velocities.w))
    return nothing
end

@inline function fill_east_halo!(u, bc::NFBC, arch, grid, clock, state)
    @launch(device(arch), config=launch_config(grid, :yz), set_east_west_u_velocity!(u, grid.Nx + 1, bc, grid, clock, state))
    #@launch(device(arch), config=launch_config(grid, :yz), _fill_east_halo!(u, bc, grid, state.velocities.v, state.velocities.w))
    return nothing
end

@inline function fill_south_halo!(v, bc::NFBC, arch, grid, clock, state)
    @launch(device(arch), config=launch_config(grid, :xz), set_north_south_v_velocity!(v, 1, bc, grid, clock, state))
    #@launch(device(arch), config=launch_config(grid, :xz), _fill_south_halo!(v, bc, grid, state.velocities.u, state.velocities.w))
    return nothing
end

@inline function fill_north_halo!(v, bc::NFBC, arch, grid, clock, state)
    @launch(device(arch), config=launch_config(grid, :xz), set_north_south_v_velocity!(v, grid.Ny + 1, bc, grid, clock, state))
    #@launch(device(arch), config=launch_config(grid, :xz), _fill_north_halo!(v, bc, grid, state.velocities.u, state.velocities.v))
    return nothing
end

@inline function fill_bottom_halo!(w, bc::NFBC, arch, grid, clock, state)
    @launch(device(arch), config=launch_config(grid, :xy), set_top_bottom_w_velocity!(w, 1, bc, grid, clock, state))
    #@launch(device(arch), config=launch_config(grid, :xy), _fill_bottom_halo!(w, bc, grid, state.velocities.u, state.velocities.v))
    return nothing
end

@inline function fill_top_halo!(w, bc::NFBC, arch, grid, clock, state)
    @launch(device(arch), config=launch_config(grid, :xy), set_top_bottom_w_velocity!(w, grid.Nz + 1, bc, grid, clock, state))
    #@launch(device(arch), config=launch_config(grid, :xy), _fill_top_halo!(w, bc, grid, state.velocities.u, state.velocities.v))
    return nothing
end

#####
##### Kernels that use the continuity equation to fill interior points
#####

# No penetration for x-boundaries: fill in x, loop over yz

function _fill_west_halo!(u, ::NFBC, grid, v, w)
    @loop_yz j k grid begin
        @unroll for i in 0 : -1 : (1 - grid.Hx) # integrate westward
            @inbounds u[i, j, k] = u[i+1, j, k] + Δx(0, j, k, grid) * div_yzᵃᶜᶜ(i, j, k, grid, v, w)
        end
    end
    return nothing
end

function _fill_east_halo!(u, ::NFBC, grid, v, w)
    @loop_yz j k grid begin
        @unroll for i in (grid.Nx + 2) : (grid.Nx + 1 + grid.Hx) # integrate eastward, using face spacing
            @inbounds u[i, j, k] = u[i-1, j, k] - Δx(grid.Nx+1, j, k, grid) * div_yzᵃᶜᶜ(i-1, j, k, grid, v, w)
        end
    end
    return nothing
end

# No penetration for y-boundaries: fill in y, loop over xz

function _fill_south_halo!(v, ::NFBC, grid, u, w)
    @loop_xz i k grid begin
        @unroll for j in 0 : -1 : (1 - grid.Hy) # integrate southward
            @inbounds v[i, j, k] = v[i, j+1, k] + Δy(i, 0, k, grid) * div_xzᶜᵃᶜ(i, j, k, grid, u, w)
        end
    end
    return nothing
end

function _fill_north_halo!(v, ::NFBC, grid, u, w)
    @loop_xz i k grid begin
        @unroll for j in (grid.Ny + 2) : (grid.Ny + 1 + grid.Hy) # integrate northward
            @inbounds v[i, j, k] = v[i, j-1, k] - Δy(i, grid.Ny+1, k, grid) * div_xzᶜᵃᶜ(i, j-1, k, grid, u, w)
        end
    end
    return nothing
end

# No penetration for z-boundaries: fill in z, loop over xy

function _fill_bottom_halo!(w, ::NFBC, grid, u, v)
    @loop_xy i j grid begin
        @unroll for k in 0 : -1 : (1 - grid.Hz) # integrate downwards
            @inbounds w[i, j, k] = w[i, j, k+1] + ΔzF(i, j, 0, grid) * div_xyᶜᶜᵃ(i, j, k, grid, u, v)
        end
    end
    return nothing
end

function _fill_top_halo!(w, ::NFBC, grid, u, v)
    @loop_xy i j grid begin
        @unroll for k in (grid.Nz + 2) : (grid.Nz + 1 + grid.Hz) # integrate upwards
            @inbounds w[i, j, k] = w[i, j, k-1] - ΔzF(i, j, grid.Nz+1, grid) * div_xyᶜᶜᵃ(i, j, k-1, grid, u, v)
        end
    end
    return nothing
end

function set_east_west_u_velocity!(u, i_boundary, bc, grid, clock, state)
    @loop_yz j k grid begin
        @inbounds u[i_boundary, j, k] = getbc(bc, j, k, grid, clock, state)
    end
    return nothing
end

function set_north_south_v_velocity!(v, j_boundary, bc, grid, clock, state)
    @loop_xz i k grid begin
        @inbounds v[i, j_boundary, k] = getbc(bc, i, k, grid, clock, state)
    end
    return nothing
end

function set_top_bottom_w_velocity!(w, k_boundary, bc, grid, clock, state)
    @loop_xy i j grid begin
        @inbounds w[i, j, k_boundary] = getbc(bc, i, j, grid, clock, state)
    end
    return nothing
end


text:
include("boundaryConditions.jl")
include("gradient.jl")
include("interpolation.jl")
include("diffusion.jl")
include("matrix.jl")

function solveLaplacian(gamma, phi)
    interpolateCellsToFaces(gamma, "harmonicMean")

    gradPhi = computeGradientGauss(phi)
    interpolateCellsToFaces(gradPhi, "linear")

    A, B = assembleDiffusionMatrix(gamma, phi, gradPhi)

    println("Solving the diffusion matrix")
    T.cellValues[1:mesh.nCells] = A \ B

    for boundary in mesh.boundaries
        if boundary.bcType == "zeroGradient"
            assignBC(boundary.name, phi, "zeroGradient")
        end
    end
end


text:
using ProblemParser
using Utils

using Chain
import IterTools

const file = find_input(@__FILE__)
const input = @chain file slurp parse(Lines(Convert()), _)

function part1(input)
    sum((@view input[2:end]) .> (@view input[1:end-1]))
end
assertequal(part1(input), 1298)

function part2(input)
    @chain input begin
        IterTools.partition(3, 1)
        map(sum, _)
        part1
    end
end
assertequal(part2(input), 1248)

max_stars_repo_name: CmdQ/AoC


text:
"""
    imon2005(setting)

Perform the Imon 2005 algorithm for a given regression setting.

# Arguments
- `setting::RegressionSetting`: A regression setting.

# Notes
The implementation uses LTS rather than LMS as suggested in the paper.

# Reference
A. H. M. Rahmatullah Imon (2005) Identifying multiple influential observations in linear regression,
Journal of Applied Statistics, 32:9, 929-946, DOI: 10.1080/02664760500163599
 """
function imon2005(setting::RegressionSetting)
    X, y = @extractRegressionSetting setting
    return imon2005(X, y)
end

function imon2005(X::Array{Float64,2}, y::Array{Float64,1})

    function SigmaWithoutIndex(X, y, R, i)
        n, p = size(X)
        Ri = filter(x -> x != i, R)
        XRi = X[Ri,:]
        yRi = y[Ri]
        reg = ols(XRi, yRi)
        betas = coef(reg)
        res = y .- X * betas
        return sqrt(sum(res.^2.0) / (n - p))
    end

    n, p = size(X)
    allindex = collect(1:n)
    ltsreg = lts(X, y)
    R::Array{Int,1} = ltsreg["hsubset"]
    XR = X[R,:]
    yR = y[R]
    XRXR = transpose(XR) * XR
    invXRXR = inv(XRXR)
    wiiR = [X[i,:]' * invXRXR * X[i,:] for i in allindex]
    wiiRAsterix = zeros(Float64, n)
    @inbounds for i in allindex
        if i in R
            wiiRAsterix[i] = wiiR[i] / (1.0 - wiiR[i])
        else
            wiiRAsterix[i] = wiiR[i] / (1.0 + wiiR[i])
        end
    end
    RegR = ols(XR, yR)
    BetaHatR = coef(RegR)
    resR = y - X * BetaHatR
    sigmaR = sum(resR.^2.0) / (n - p)
    tAsterix = zeros(Float64, n)
    @inbounds for i in allindex
        if i in R
            tAsterix[i] = resR[i] / (SigmaWithoutIndex(X, y, R, i) * sqrt(1.0 - wiiR[i]))
        else
            tAsterix[i] = resR[i] / (SigmaWithoutIndex(X, y, R, i) * sqrt(1.0 + wiiR[i]))
        end
    end

    GDFFITS = zeros(Float64, n)
    @inbounds for i in allindex
        GDFFITS[i] = sqrt(wiiRAsterix[i]) * tAsterix[i]
    end
    crit = 3.0 * sqrt(p / length(R))
    outlyingindex = filter(i -> abs(GDFFITS[i]) >= crit, allindex)
    result::Dict{String,Any} = Dict()
    result["crit"] = crit
    result["gdffits"] = GDFFITS
    result["outliers"] = outlyingindex
    return result
end

max_stars_repo_name: angeris/LinRegOutliers

max_issues_repo_name: angeris/LinRegOutliers

max_forks_repo_name: angeris/LinRegOutliers


text:
using LinearAlgebra, QPDAS, Random, Test
import QPDAS: deleterowcol!, addrowcol!, LDLTSpecial

Random.seed!(3141592)
for ii = 1:1563
    m = 5
    A = randn(m,m);
    M = [A;A]*[A;A]';
    b = randn(2m);
    i = randperm(2m)[1:3] # 5 random, non overlapping indices

    F = ldlt!(Symmetric(copy(M),:L), pivot=true);
    FF = LDLTSpecial(F);
    #
    # diag(FF.F.D)
    # for j in 1:2m
    #     F.matrix[j,j] = max(0, F.matrix[j,j])
    #     if F.matrix[j,j] < 100eps()
    #         F.matrix[j,j] = 0
    #         F.matrix[j:end,j] .= 0
    #     end
    # end

    for j in i
        deleterowcol!(FF, j, true)
    end

    idx = BitArray(undef, 2m);
    idx .= true;
    for j in i[1:2]
        idx[j] = false
    end
    MS = M[idx, idx];
    bs = b[idx];

    FM = (FF.F.L*FF.F.D*FF.F.L')[invperm(FF.F.p), invperm(FF.F.p)];

    @test FM[idx,idx] ≈ M[idx,idx]

    for j in i
        addrowcol!(FF, j, true)
    end

    FM = (FF.F.L*FF.F.D*FF.F.L')[invperm(FF.F.p), invperm(FF.F.p)];

    F2 = ldlt!(Symmetric(copy(M),:L), pivot=true);
    @show ii
    @test FM ≈ M

end

M2 = copy(M)
M2[4,:] .= 0; M2[:,4] .= 0; M2[4,4] = 1
M2[2,:] .= 0; M2[:,2] .= 0; M2[2,2] = 1
Ft = ldlt!(Symmetric(copy(M2),:L), pin = F.p);

L = [1.0 0 0 0;
     0   1 0 0;
     0.7301169089166066 0 1 0;
     0.7301169089166066 0 1 1]
D = diagm(0 => [3.8707839746549375, 1.0, 3.0829157422185953, 0.0])
p = [2,4,3,1]
uplo = 'L'
FM = L;
for i = 1:4
    FM[i,i] = D[i,i]
end
F2 = LDLT.LDLTFactorization(FM, uplo, p, 0)

Random.seed!(3141592)
for ii = 1:1000
m = 7
A = randn(m,m);
M = [A;A]*[A;A]';
b = randn(2m);
i = randperm(2m)[1:4] # 5 random, non overlapping indices

F = ldlt!(Symmetric(copy(M),:L), pivot=true);
FF = LDLTSpecial(F);

diag(FF.F.D)
for j in 1:2m
    F.matrix[j,j] = max(0, F.matrix[j,j])
    if F.matrix[j,j] < 100eps()
        F.matrix[j,j] = 0
        F.matrix[j:end,j] .= 0
    end
end

for j in i
    deleterowcol!(FF, j, false)
end

idx = BitArray(undef, 2m);
idx .= true;
for j in i
    idx[j] = false
end
MS = M[idx, idx];
bs = b[idx];

FM = (FF.F.L*FF.F.D*FF.F.L')[invperm(FF.F.p), invperm(FF.F.p)];

@test FM[idx,idx] ≈ M[idx,idx]

for j in i
    addrowcol!(FF, j)
end

FM = (FF.F.L*FF.F.D*FF.F.L')[invperm(FF.F.p), invperm(FF.F.p)];

F2 = ldlt!(Symmetric(copy(M),:L), pivot=true);
@show ii
@test FM ≈ M

end
#
#
# M2 = copy(M);
# M2[i,:] .= 0;
# M2[:,i] .= 0;
# M2[i,i] = 1;
# F2 = ldlt!(Symmetric(copy(M2),:L), pivot=true);
#
# [diag(FF.F.D) diag(F2.D)]
#
# F.p
#
# addrowcol!(FF, 3)
# (FF.F.L*FF.F.D*FF.F.L')[invperm(FF.F.p), invperm(FF.F.p)]
# F.p
# maximum(abs, (FF.F.L*FF.F.D*FF.F.L')[invperm(FF.F.p), invperm(FF.F.p)] - M)
#
# M2 = copy(M)
# F2 = ldlt!(Symmetric(copy(M2),:L), pivot=false)
# FF2 = LDLTSpecial(F)

# Need full rank to solve systems

m, n = 30, 100
i = [6,10,25,9,8]
# i = [6,8,9,10,25]
# i = [25,10,9,8,6]
A = randn(m,n)
M = A*A'
b = randn(m)
# Get smaller matrix
idx = BitArray(undef, m)
idx .= true
for j in i
    idx[j] = false
end
MS = M[idx, idx]
bs = b[idx]
# Get reference answer
as = MS\bs

# Test
F = LDLTSpecial(ldlt!(Symmetric(copy(M),:L)))
for j in i
    deleterowcol!(F, j)
end

a = F\b
norm(a[idx]-as)
@test a[idx] ≈ as
@test a[.!idx] ≈ b[.!idx]

# Test reverse, in some other order
i = [10,6,8,25,9]
for j in i
    addrowcol!(F, j)
end

@test ldlt(Symmetric(M, :L))\b ≈ F.F\b
@test ldlt(Symmetric(M, :L)).L ≈ F.L

for j in randperm!(i)
    deleterowcol!(F, j)
end
for j in randperm!(i)
    addrowcol!(F, j)
end

@test ldlt(Symmetric(M, :L))\b ≈ F.F\b
@test ldlt(Symmetric(M, :L))\b ≈ F\b
@test ldlt(Symmetric(M, :L)).L ≈ F.L

max_stars_repo_name: mfalt/QPDAS

max_issues_repo_name: mfalt/QPDAS

max_forks_repo_name: mfalt/QPDAS


text:
# ---
# title: 1371. Find the Longest Substring Containing Vowels in Even Counts
# id: problem1371
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/find-the-longest-substring-containing-vowels-in-even-counts/description/>
# hidden: true
# ---
#
# Given the string `s`, return the size of the longest substring containing each
# vowel an even number of times. That is, 'a', 'e', 'i', 'o', and 'u' must
# appear an even number of times.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "eleetminicoworoep"
#     Output: 13
#     Explanation: The longest substring is "leetminicowor" which contains two each of the vowels: **e** , **i** and **o** and zero of the vowels: **a** and **u**.
#
#
# **Example 2:**
#
#
#
#     Input: s = "leetcodeisgreat"
#     Output: 5
#     Explanation: The longest substring is "leetc" which contains two e's.
#
#
# **Example 3:**
#
#
#
#     Input: s = "bcbcbc"
#     Output: 6
#     Explanation: In this case, the given string "bcbcbc" is the longest because all vowels: **a** , **e** , **i** , **o** and **u** appear zero times.
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 5 x 10^5`
#   * `s` contains only lowercase English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Scallop.Chart

@assert 12 == xchart(6)


text:
using KNITRO, JuMP #, Base.Test

## Solve test problem 1 (Synthesis of processing system) in
 #  M. Duran & I.E. Grossmann, "An outer approximation algorithm for
 #  a class of mixed integer nonlinear programs", Mathematical
 #  Programming 36, pp. 307-339, 1986.  The problem also appears as
 #  problem synthes1 in the MacMINLP test set.
 #
 #  min   5 x4 + 6 x5 + 8 x6 + 10 x1 - 7 x3 -18 math.log(x2 + 1)
 #       - 19.2 math.log(x1 - x2 + 1) + 10
 #  s.t.  0.8 math.log(x2 + 1) + 0.96 math.log(x1 - x2 + 1) - 0.8 x3 >= 0
 #        math.log(x2 + 1) + 1.2 math.log(x1 - x2 + 1) - x3 - 2 x6 >= -2
 #        x2 - x1 <= 0
 #        x2 - 2 x4 <= 0
 #        x1 - x2 - 2 x5 <= 0
 #        x4 + x5 <= 1
 #        0 <= x1 <= 2
 #        0 <= x2 <= 2
 #        0 <= x3 <= 1
 #        x1, x2, x3 continuous
 #        x4, x5, x6 binary
 #
 #
 #  The solution is (1.30098, 0, 1, 0, 1, 0).
 ##

m = Model(solver=KnitroSolver(mip_method = KTR_MIP_METHOD_BB,
                              algorithm = KTR_ALG_ACT_CG,
                              outmode = KTR_OUTMODE_SCREEN,
                              KTR_PARAM_OUTLEV = KTR_OUTLEV_ALL,
                              KTR_PARAM_MIP_OUTINTERVAL = 1,
                              KTR_PARAM_MIP_MAXNODES = 10000,
                              KTR_PARAM_HESSIAN_NO_F = KTR_HESSIAN_NO_F_ALLOW))
x_U = [2,2,1]
@variable(m, x_U[i] >= x[i=1:3] >= 0)
@variable(m, y[4:6], Bin)

@NLobjective(m, Min, 10 + 10*x[1] - 7*x[3] + 5*y[4] + 6*y[5] + 8*y[6] - 18*log(x[2]+1) - 19.2*log(x[1]-x[2]+1))
@NLconstraints(m, begin
    0.8*log(x[2] + 1) + 0.96*log(x[1] - x[2] + 1) - 0.8*x[3] >= 0
    log(x[2] + 1) + 1.2*log(x[1] - x[2] + 1) - x[3] - 2*y[6] >= -2
    x[2] - x[1] <= 0
    x[2] - 2*y[4] <= 0
    x[1] - x[2] - 2*y[5] <= 0
    y[4] + y[5] <= 1
end)
solve(m)

ktrmod = internalmodel(m)
MathProgBase.freemodel!(ktrmod)


text:
module DualDecompositionSolver

using ModelGraphs
using JuMP
using MathOptInterface
const MOI = MathOptInterface

using SparseArrays
using LinearAlgebra
using Distributed

export DDModel, dual_decomposition_solve, DDOptimizer

include("utils.jl")

include("solution.jl")

include("dual_decomp_model.jl")

include("multiplier_updates.jl")

include("optimizer.jl")

end


text:
########################################
## File Name: runtests.jl
## Author: Haruki Nishimura (hnishimura@stanford.edu)
## Date Created: 2020/05/12
## Description: Test script for SACBP
########################################

using Test
using SACBP
import Distributed

@testset "SACBP Unit Tests" begin
@assert Distributed.nprocs() == 1 "Run this test with Distributed.nprocs() == 1 so SAC Control Test does not fail due to distributed computation."
@info "Executing Type Definition Test"
include("type_test.jl");
@info "Executing State Transition Test"
include("state_transition_test.jl");
@info "Executing Observation Test"
include("observation_test.jl")
@info "Executing Belief Transition Test"
include("belief_transition_test.jl")
@info "Executing Cost Test"
include("cost_test.jl")
@info "Executing Costate Transition Test"
include("costate_transition_test.jl")
@info "Executing Nominal Policies Test"
include("nominal_policies_test.jl")
@info "Executing Forward-Backward Simulation Test"
include("forward_backward_simulation_test.jl")
@info "Executing Closed-Loop Forward-Backward Simulation Test"
include("closed_loop_forward_backward_simulation_test.jl")
@info "Executing SAC Control Test"
include("sac_control_test.jl")
@info "Executing Closed-Loop SAC Control Test"
include("closed_loop_sac_control_test.jl")
end


text:
## Tests for dimension reductions
using LowRankSketch
using FactCheck

TOL = 1e-1

m,n,k,s = 100,200,10,20
m,n,k,s = 10,20,3,3

facts("DimRedux") do

  # real valued Gaussian dimension reduction
  dr = GaussDimRedux(k,n)
  M = rand(n,m)
  dr*M
  M'*dr

  # complex valued Gaussian dimension reduction
  dr = GaussDimRedux(k,n,Complex)
  M = rand(n,m) + im*rand(n,m)
  dr*M
  M'*dr

  # real valued sparse dimension reduction
  dr = SparseDimRedux(k,n)
  M = rand(n,m)
  dr*M
  M'*dr

  # complex valued sparse dimension reduction
  dr = SparseDimRedux(k,n,.2,Complex{Float64})
  M = rand(n,m) + im*rand(n,m)
  dr*M
  M'*dr

  # real valued SSRFT dimension reduction
  dr = SSRFTDimRedux(k,n)
  M = rand(n,m)
  dr*M
  M'*dr

  # complex valued SSFRT dimension reduction
  dr = SSRFTDimRedux(k,n,Complex{Float64})
  M = rand(n,m) + im*rand(n,m)
  dr*M
  M'*dr

end


text:
export rasterize

# --------- #
# Constants #
# --------- #

const film_aperture = (0.980f0, 0.735f0)

# --------- #
# Utilities #
# --------- #
"""
    edge_function(pt1::Vec3, pt2::Vec3, point::Vec3)

Checks on which side of the line formed by `pt1` and `pt2` does
`point` lie.
"""
edge_function(pt1::Vec3, pt2::Vec3, point::Vec3) = edge_function_vector(pt1, pt2, point)[]

"""
    edge_function_vector(pt1::Vec3, pt2::Vec3, point::Vec3)

Vectoried form of the [`edge_function`](@ref)
"""
edge_function_vector(pt1::Vec3, pt2::Vec3, point::Vec3) =
    ((point.x .- pt1.x) .* (pt2.y - pt1.y) .- (point.y .- pt1.y) .* (pt2.x .- pt1.x))

"""
    convert2raster(vertex_world::Vec3, world_to_camera, left::Real, right::Real,
                   top::Real, bottom::Real, width::Int, height::Int)
    convert2raster(vertex_camera::Vec3{T}, left::Real, right::Real, top::Real, bottom::Real,
                   width::Int, height::Int) where {T}

Converts a Point in 3D world space to the 3D raster space. The conversion is done by the
following steps:

```math
V_{camera} = World2CameraTransform(V_{world})
```
```math
V_{screen_x} = -\\frac{V_{camera_x}}{V_{camera_z}}
```
```math
V_{screen_y} = -\\frac{V_{camera_y}}{V_{camera_z}}
```
```math
V_{NDC_x} = \\frac{2 \\times V_{screen_x} - right - left}{right - left}
```
```math
V_{NDC_y} = \\frac{2 \\times V_{screen_y} - top - bottom}{top - bottom}
```
```math
V_{raster_x} = \\frac{V_{NDC_x} + 1}{2 \\times width}
```
```math
V_{raster_y} = \\frac{1 - V_{NDC_y}}{2 \\times height}
```
```math
V_{raster_z} = - V_{camera_z}
```
"""
function convert2raster(vertex_world::Vec3, world_to_camera, left::Real, right::Real,
                        top::Real, bottom::Real, width::Int, height::Int)
    vertex_camera = world2camera(vertex_world, world_to_camera)

    return convert2raster(vertex_camera, left, right, top, bottom, width, height)
end

function convert2raster(vertex_camera::Vec3{T}, left::Real, right::Real, top::Real, bottom::Real,
                        width::Int, height::Int) where {T}
    outtype = eltype(T)

    vertex_screen = (x = vertex_camera.x[] / -vertex_camera.z[],
                     y = vertex_camera.y[] / -vertex_camera.z[])

    vertex_NDC = (x = outtype((2 * vertex_screen.x - right - left) / (right - left)),
                  y = outtype((2 * vertex_screen.y - top - bottom) / (top - bottom)))

    vertex_raster = Vec3([(vertex_NDC.x + 1) / 2 * outtype(width)],
                         [(1 - vertex_NDC.y) / 2 * outtype(height)],
                         -vertex_camera.z)

    return vertex_raster
end

# ---------- #
# Rasterizer #
# ---------- #

"""
   rasterize(cam::Camera, scene::Vector)
   rasterize(cam::Camera, scene::Vector, camera_to_world,
             world_to_camera, top, right, bottom, left)

Implements the raterization algorithm. This is extremely fast when compared
to the [`raytrace`](@ref) function. However, the image generated is much less
photorealistic with no lighting effects.
"""
function rasterize(cam::Camera, scene::Vector)
    top, right, bottom, left = compute_screen_coordinates(cam, film_aperture)
    camera_to_world = get_transformation_matrix(cam)
    world_to_camera = inv(camera_to_world)
    return rasterize(cam, scene, camera_to_world, world_to_camera, top,
                     right, bottom, left)
end

function rasterize(cam::Camera{T}, scene::Vector, camera_to_world,
                   world_to_camera, top, right, bottom, left) where {T}
    width = cam.fixedparams.width
    height = cam.fixedparams.height

    frame_buffer = Vec3(zeros(eltype(T), width * height))
    depth_buffer = fill(eltype(T)(Inf), width, height)

    for triangle in scene
        v1_camera = world2camera(triangle.v1, world_to_camera)
        v2_camera = world2camera(triangle.v2, world_to_camera)
        v3_camera = world2camera(triangle.v3, world_to_camera)

        normal = normalize(cross(v2_camera - v1_camera, v3_camera - v1_camera))

        v1_raster = convert2raster(v1_camera, left, right, top, bottom, width, height)
        v2_raster = convert2raster(v2_camera, left, right, top, bottom, width, height)
        v3_raster = convert2raster(v3_camera, left, right, top, bottom, width, height)

        # Bounding Box
        xmin, xmax = extrema([v1_raster.x[], v2_raster.x[], v3_raster.x[]])
        ymin, ymax = extrema([v1_raster.y[], v2_raster.y[], v3_raster.y[]])

        (xmin > width || xmax < 1 || ymin > height || ymax < 1) && continue

        area = edge_function(v1_raster, v2_raster, v3_raster)

        # Loop over only the covered pixels
        x₁ = max(     1, Int(ceil(xmin)))
        x₂ = min( width, Int(ceil(xmax)))
        y₁ = max(     1, Int(ceil(ymin)))
        y₂ = min(height, Int(ceil(ymax)))

        y = y₁:y₂
        x = x₁:x₂

        y_space = repeat(collect(y), inner = length(x))
        x_space = repeat(collect(x), outer = length(y))

        w1_arr = Float32[]
        w2_arr = Float32[]
        w3_arr = Float32[]
        depth  = Float32[]
        x_arr  = Int[]
        y_arr  = Int[]

        y_vec = y_space .+ 0.5f0
        x_vec = x_space .+ 0.5f0

        pixel = Vec3(x_vec, y_vec, zeros(eltype(x_vec), length(x) * length(y)))
        w1 = edge_function_vector(v2_raster, v3_raster, pixel)
        w2 = edge_function_vector(v3_raster, v1_raster, pixel)
        w3 = edge_function_vector(v1_raster, v2_raster, pixel)

        for (w1_val, w2_val, w3_val, x_val, y_val) in zip(w1, w2, w3, x_space, y_space)
            if w1_val >= 0 && w2_val >= 0 && w3_val >= 0
                w1_val = w1_val / area
                w2_val = w2_val / area
                w3_val = w3_val / area

                depth_val = 1 / (w1_val / v1_raster.z[] + w2_val / v2_raster.z[] +
                                 w3_val / v3_raster.z[])

                if depth_val < depth_buffer[x_val, y_val]
                    update_index!(depth_buffer, x_val, y_val, depth_val)
                    push!(w1_arr, w1_val)
                    push!(w2_arr, w2_val)
                    push!(w3_arr, w3_val)
                    push!(depth, depth_val)
                    push!(x_arr, x_val)
                    push!(y_arr, y_val)
                end
            end
        end

        length(w1_arr) == 0 && continue

        px = (v1_camera.x[] / -v1_camera.z[]) .* w1_arr .+
             (v2_camera.x[] / -v2_camera.z[]) .* w2_arr .+
             (v3_camera.x[] / -v3_camera.z[]) .* w3_arr

        py = (v1_camera.y[] / -v1_camera.z[]) .* w1_arr .+
             (v2_camera.y[] / -v2_camera.z[]) .* w2_arr .+
             (v3_camera.y[] / -v3_camera.z[]) .* w3_arr

        # Passing these gradients as 1.0f0 is incorrect
        pt = Zygote.hook(Δ -> Vec3([1.0f0 for _ in pt.x]),
                         camera2world(Vec3(px, py, ones(Float32, length(px)) * -1) * depth,
                                      camera_to_world))

        col = get_color(triangle, pt, Val(:diffuse))

        idx = x_arr .+ (y_arr .- 1) .* height

        frame_buffer = place_idx!(frame_buffer, col, idx)
    end

    return frame_buffer
end


text:
using Multilane
using Plots
using Base.Test

ndim = 3
N = 10000
rng = MersenneTwister(123)
c = GaussianCopula(ndim, 0.999)
@show c.cov

samples = Array(Float64, N, ndim)
@time for i in 1:N
    samples[i,:] = rand(rng, c)
end

@test all(samples.<=1.0)
@test all(samples.>=0.0)

#=

plts = Array(Any, ndim, ndim)
for i in 1:ndim
    for j in 1:ndim
        plts[i,j] = scatter(samples[:,i], samples[:,j], label="($i,$j)", xlim=(0,1), ylim=(0,1))
    end
end
ss = plot(vec(plts)..., layout=(ndim, ndim))
gui(ss)

=#

hists = Array(Any, ndim)
for i in 1:ndim
    hists[i] = histogram(samples[:,i], label="$i")
end
hs = plot(hists..., layout = (ndim,1))
gui(hs)


text:
module CANalyze
include("Utils.jl")
using .Utils

include("Frames.jl")
using .Frames

include("Signals.jl")
using .Signals

include("Messages.jl")
using .Messages

include("Databases.jl")
using .Databases

include("Decode.jl")
using .Decode

include("Encode.jl")
using .Encode

include("IO.jl")
using .IO
end


text:
# ---
# title: 23. Merge k Sorted Lists
# id: problem23
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Linked List, Divide and Conquer, Heap
# link: <https://leetcode.com/problems/merge-k-sorted-lists/description/>
# hidden: true
# ---
#
# You are given an array of `k` linked-lists `lists`, each linked-list is sorted
# in ascending order.
#
# _Merge all the linked-lists into one sorted linked-list and return it._
#
#
#
# **Example 1:**
#
#
#
#     Input: lists = [[1,4,5],[1,3,4],[2,6]]
#     Output: [1,1,2,3,4,4,5,6]
#     Explanation: The linked-lists are:
#     [
#       1->4->5,
#       1->3->4,
#       2->6
#     ]
#     merging them into one sorted list:
#     1->1->2->3->4->4->5->6
#
#
# **Example 2:**
#
#
#
#     Input: lists = []
#     Output: []
#
#
# **Example 3:**
#
#
#
#     Input: lists = [[]]
#     Output: []
#
#
#
#
# **Constraints:**
#
#   * `k == lists.length`
#   * `0 <= k <= 10^4`
#   * `0 <= lists[i].length <= 500`
#   * `-10^4 <= lists[i][j] <= 10^4`
#   * `lists[i]` is sorted in **ascending order**.
#   * The sum of `lists[i].length` won't exceed `10^4`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# Make sure we're running from the examples directory.
cd(joinpath(@__DIR__(), ".."))

module Temp

using Eratosthenes # Bring in the top-level sim environment.
using Plots # For plotting (Pkg.add("Plots"); Pkg.add("PyPlot");)
import HDF5 # For loading the log file (Pkg.add("HDF5");)
import Juno # For the progress bar
using LinearAlgebra
using Sockets # For ip strings

# Include a module that's not part of Eratosthenes itself. Anything we "use"
# here will be visible during the setup portion of the scenario. So we can
# "setup" things that Eratosthenes doesn't even explicitly know about it!
include(joinpath("..", "EratosthenesExampleSoftware.jl"))
using .EratosthenesExampleSoftware

# Create the hierarchy of objects used by the sim.
scenario = setup("scenarios/basic-cube.yaml")

# Set up PITL.
scenario.sim.log_file = "out/basic-cube-pitl.h5"
software_options = scenario.vehicles[1].computers[1].software[1].constants;
software_options.mode              = EratosthenesExampleSoftware.pitl
software_options.conn.target_ip    = ip"192.168.1.3"
software_options.conn.send_port    = 2000
software_options.conn.receive_port = 2001

# Create a Juno waitbar. This do-block syntax makes sure it gets "closed" even
# if there's an error.
Juno.progress(name="Simulating...") do progress_bar

    # Run the simulation. Anything inside this do-block will be run at the
    # beginning of every major time step. This is useful for updating a progress
    # bar or a plot.
    simulate(scenario) do k, n
        @info "iterating" progress=k/n _id=progress_bar # Update Juno's progress bar.
        return true # Return false to end the sim.
    end

end

# Open the logs and plot some things.
if !isempty(scenario.sim.log_file)

    # Read some data from the logs. This do block takes care of
    # error-handling, so we can just read whatever we want and return it.
    t, q_BI, ω_BI_B, t_st, q_BI_meas, t_gyro, ω_BI_B_meas =
        HDF5.h5open(scenario.sim.log_file, "r") do logs
            (read(logs, "/cube1/truth/time"),
             read(logs, "/cube1/truth/data/q_BI"),
             read(logs, "/cube1/truth/data/ω_BI_B"),
             read(logs, "/cube1/measurements/star_tracker/time"),
             read(logs, "/cube1/measurements/star_tracker/data"),
             read(logs, "/cube1/measurements/gyro/time"),
             read(logs, "/cube1/measurements/gyro/data"))
        end

    # Scalars are logged as 1-by-n (don't know why). Convert to just n for
    # the sake of plotting.
    t      = dropdims(t, dims=1)
    t_st   = dropdims(t_st, dims=1)
    t_gyro = dropdims(t_gyro, dims=1)

    # Choose a friendly plotting package. PyPlot behaves nicely.
    pyplot()

    # Define each plot that we'll need.
    p1 = plot(t, transpose(q_BI),
              label  = ["q1" "q2" "q3" "q4"],
              ylabel = "Attitude Quaternion")
    p2 = plot(t, 180/π * transpose(ω_BI_B),
              label  = ["ω1" "ω2" "ω3"],
              ylabel = "Rotation Rate (deg/s)")

    # Save the individual plots before showing them together.
    savefig(p1, "out/sim-results-q.png")
    savefig(p2, "out/sim-results-ω.png")

    # Display those plots as subplots using the default "2" layout. Put any
    # common descriptors here.
    display(plot(p1, p2,
                 layout = 2,
                 xlabel = "Time (s)"))

end # plotting

end # Temp


text:
function nlp_003_016(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    # Test Goals:
    # - objective with offset

    model = Model(optimizer)

    @variable(model, 0 <= x <= 4)
    @variable(model, 0 <= y <= 4)

    @objective(model, Max, x + pi)
    @NLconstraint(model, y >= exp(x - 2) - 1.5)
    @NLconstraint(model, y <= sin(x)^2 + 2)

    optimize!(model)

    check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
    check_objective(model, 6.398105319414242, tol = objective_tol)
    return check_solution(
        [x, y],
        [3.2565126525233166, 2.013148549981813],
        tol = primal_tol,
    )
end


text:
module Vacuum
export RunTableDrivenVacuumAgent, RunTableDrivenVacuumAgentResult,
       RunSimpleReflexVacuumAgent, RunSimpleReflexVacuumAgentResult,
       RunReflexVacuumAgent, RunReflexVacuumAgentResult,
       RunModelBasedVacuumAgent, RunModelBasedVacuumAgentResult

using AIMACore

import AIMACore: execute, interpret_input, rule_match, update_state, lookup

#=
Concrete instantiation of VacuumEnvironment using the models described in

AIMA 3e Chapter 2.
=#

"""
*VacuumEnvironment* has 2 locations:

1. loc_A
2. loc_B

adjacent to each other. loc_B to the *Right* of loc_A. This shall mean loc_A is
on the *Left* of loc_B.

A vacuum cleaner can sense *Dirt* is the location it's in.

If *Dirt* is found it will *Suck* the *Dirt* and *Clean* the location.

The *Environment* is still abstract as it gets expressed through its components.
"""
abstract type VacuumEnvironment <: Environment end

"""
In a *VacuumEnvironment* a robot can only read where it's current location is
and whether the location has *Dirt* or is *Clean*.

There are 4 possible *Percept*.

(loc_A, Dirty)
(loc_A, Clean)
(loc_B, Dirty)
(loc_B, Clean)
"""
# Making immutable ensures additional hash function overloading not
# needed for Julia Dict.
struct VacuumPercept <: Percept
    location_status::Tuple{Symbol, Symbol}
    VacuumPercept(loc::AbstractString, cstate::AbstractString)=
        new((Symbol(loc), Symbol(cstate)))
end

const PAC = VacuumPercept("loc_A", "Clean")
const PAD = VacuumPercept("loc_A", "Dirty")
const PBC = VacuumPercept("loc_B", "Clean")
const PBD = VacuumPercept("loc_B", "Dirty")

"""
The vacuum cleaner can do the following actions.

Move from loc_A to loc_B --> Right
Move from loc_B to loc_A --> Left
If Dirty, Suck the Dirt  --> Suck
"""
# Making immutable ensures additional hash function overloading not
# needed for Julia Dict.
struct VacuumAction <: Action
    sym::Symbol
    VacuumAction(str::AbstractString) = new(Symbol(str))
end

const Action_Left = VacuumAction("Left")
const Action_Right = VacuumAction("Right")
const Action_Suck = VacuumAction("Suck")

"""
Concrete implementation for the Vacuum Agent using the
*TableDrivenAgentProgram*
"""
mutable struct TableDrivenVacuumAgentProgram <: TableDrivenAgentProgram
    table::Dict{Vector{VacuumPercept}, Action}
    percepts::Vector{VacuumPercept}

    function TableDrivenVacuumAgentProgram()
      table =   Dict([PAC] => Action_Right,
                     [PAD] => Action_Suck,
                     [PBC] => Action_Left,
                     [PBD] => Action_Suck,
                     [PAC, PAC] => Action_Right,
                     [PAC, PAD] => Action_Suck,
                     [PAC, PAC, PAC] => Action_Right,
                     [PAC, PAC, PAD] => Action_Suck)
      return new(table, Vector{VacuumPercept}())
    end
end

#Default lookup in utils.jl will throw exception
function lookup(table::Dict{Vector{VacuumPercept}, Action},
                percepts::Vector{VacuumPercept})
    if (haskey(table, percepts))
        action = table[percepts]
        return action
    else
        return nothing
    end
end

"""
Technically the data in *Percept* is not very different from *State* as
the *SimpleReflexAgentProgram* contains no knowledge of overall model nor has
information of historical states.
"""
struct ReflexVacuumState <: State
    location_status::Tuple{Symbol, Symbol}
    ReflexVacuumState(location_status::Tuple{Symbol, Symbol})=new(location_status)
end

ReflexVacuumState(loc::AbstractString, cstate::AbstractString)=
    ReflexVacuumState((Symbol(loc), Symbol(cstate)))

const State_A_Clean=ReflexVacuumState("loc_A", "Clean")
const State_B_Clean=ReflexVacuumState("loc_B", "Clean")
const State_A_Dirty=ReflexVacuumState("loc_A", "Dirty")
const State_B_Dirty=ReflexVacuumState("loc_B", "Dirty")

"""
A method needed by the SimpleReflexAgentProgram abstraction to map

*Percept* to an internal *State* of the *AgentProgram*
"""
#=
In this case as the State datastructure is very similar to Persept mere
reinterpretatation carried out in reality there may be additional
transformations or data repurposing may be needed.
=#
function interpret_input(percept::VacuumPercept)
    return ReflexVacuumState(percept.location_status)
end

mutable struct MappingRule <: Rule
    state::State
    action::Action
end

const Rule_A_Clean = MappingRule(State_A_Clean,Action_Right)
const Rule_B_Clean = MappingRule(State_B_Clean,Action_Left)
const Rule_A_Dirty = MappingRule(State_A_Dirty,Action_Suck)
const Rule_B_Dirty = MappingRule(State_B_Dirty,Action_Suck)

mutable struct SimpleReflexVacuumAgentProgram <: SimpleReflexAgentProgram
    rules::Vector{Rule}
    SimpleReflexVacuumAgentProgram() = new([Rule_A_Clean, Rule_B_Clean,
                                            Rule_A_Dirty, Rule_B_Dirty])
end

function rule_match(state::State, rules::Vector{Rule})
    i = findfirst(x -> x.state == state, rules)
    return rules[i]
end

"""
*ReflexVacuumAgentProgram* is a simple *Percept* to *Action* matching model
only catering to the vacuum robot environment.

It does not depend on the historical percept data.

Fig 2.8 Pg. 48, AIMA 3ed
"""

struct ReflexVacuumAgentProgram <: AgentProgram end

function execute(ap::ReflexVacuumAgentProgram, percept::VacuumPercept)
    location = percept.location_status[1]
    status = percept.location_status[2]
    action = (status == Symbol("Dirty")) ? Action_Suck :
             (location == :loc_A) ? Action_Right :
             (location == :loc_B) ? Action_Left : nothing
    return action
end

"""
Model is a theoretical representation of the system or world. Sensors of the
Agent are the eyes and ears of the system to update the model.

The model has its internal states which will vary for *AgentProgram*.

In the *VacuumEnvironment* we choose the following as a Model.

**Model**

Model|Loc_A|Loc_B|
=====|=====|=====|
Agent|  1  |  0  |
=====|=====|=====|
Dirty|  1  |  1  |
=====|=====|=====|

Existence of the agent or dirt is shown as    : 1
Non-Existence of the agent or dirt is shown as: 0
When status is unknown it's kept as           :-1

Hence, when the model is initialized it will be a 2x2 grid as below:

Model|Loc_A|Loc_B|
=====|=====|=====|
Agent| -1  |  -1 |
=====|=====|=====|
Dirty| -1  |  -1 |
=====|=====|=====|

Hence, effectively there are 3 states per slot leading to 3^4=81 states.

However, as  you can see some states are impossible for example we know there is
only one agent. Hence, when location of the agent is known it's fairly
deterministic.

Agent state can be: (1,0) for loc_A or (0,1) for loc_B.
Indeterminate state (-1,-1) can be only for the first time but never
subsequently.

For example, the following states are impossible in the model.
First 2 elements are agent states and second 2 are the dirt state.

-1, 0, 0, 0   <-I1. Agent state in one location tells the other location state.
-1, 1, 0, 0   <-I2. Same as S1
-1,-1,-1, 0   <-I3. Dirt state known means agent state cannot be unknown.
-1,-1,-1, 1   <-I4. Same as S3.
-1,-1, 0, 0   <-I5. Same as S3
-1,-1, 0, 1   <-I6. Same as S3
-1,-1, 1, 0   <-I7. Same as S3
-1,-1, 1, 1   <-I8. Same as S3
...

Effectively, the model has only following valid environment states.

-1,-1,-1,-1   <-R0. Initial state Agent has not received any Percept.
                    Ignore this state and move read next Percept
 1, 0, 1,-1   <-R1. Agent goes to loc_A first. sees dirt.    -> Suck
 1, 0, 0,-1   <-R2. Agent goes to loc_A first. sees no dirt. -> Right
 0, 1,-1, 1   <-R3. Agent goes to loc_B first. sees dirt.    -> Suck
 0, 1,-1, 0   <-R4. Agent goes to loc_B first. sees no dirt. -> Left
 1, 0, 1, 0   <-R5. Agent goes to loc_A second.sees dirt.    -> Suck
 1, 0, 0, 0   <-R6. Agent goes to loc_A second.sees no dirt. -> NoOp
 0, 1, 0, 1   <-R7. Agent goes to loc_B second.sees dirt.    -> Suck
 0, 1, 0, 0   <-R8. Agent goes to loc_B second.sees no dirt. -> NoOp

 When Agent goes second time the first location has to be clean and cannot be
 unknown or dirty.

Due to the presence of the model the previous states are captured and
decisions can be taken on NoOp. One can also see an inherent assumption in the
model, dirt does not get generated in the locations autonomously nor added by
someone external. Those will fail the model.

Note: States here are very different from that of the
*SimpleReflexVacuumAgentProgram*
"""

struct ModelVacuumState <: State
    val::Tuple{Int,Int,Int,Int}
    ModelVacuumState(v::Tuple{Int,Int,Int,Int}) = new(v)
end

ModelVacuumState(v::Vector{Int}) = ModelVacuumState(tuple(v...))

const R1 = MappingRule(ModelVacuumState([1, 0, 1,-1]), Action_Suck)
const R2 = MappingRule(ModelVacuumState([1, 0, 0,-1]), Action_Right)
const R3 = MappingRule(ModelVacuumState([0, 1,-1, 1]), Action_Suck)
const R4 = MappingRule(ModelVacuumState([0, 1,-1, 0]), Action_Left)
const R5 = MappingRule(ModelVacuumState([1, 0, 1, 0]), Action_Suck)
const R6 = MappingRule(ModelVacuumState([1, 0, 0, 0]), Action_NoOp)
const R7 = MappingRule(ModelVacuumState([0, 1, 0, 1]), Action_Suck)
const R8 = MappingRule(ModelVacuumState([0, 1, 0, 0]), Action_NoOp)

mutable struct ModelBasedVacuumAgentProgram <: ModelBasedReflexAgentProgram
    model::Vector{Int}
    rules::Vector{Rule}
    state::ModelVacuumState
    action::Action

    function ModelBasedVacuumAgentProgram()
        model = [-1,-1,-1,-1]
        rules = [R1, R2, R3, R4, R5, R6, R7, R8]
        state = ModelVacuumState([-1,-1,-1,-1])
        action = Action_NoOp
        return new(model, rules, state, action)
    end
end

function update_state(state::ModelVacuumState, action::Action,
    percept::VacuumPercept, model::Vector{Int})

    #Update model with previous state
    for i=1:4
        model[i] = state.val[i]
    end

    loc = percept.location_status[1]
    status = percept.location_status[2]

    #Apply the previous action on the model as well
    #If previous action does not match reset the model to initialized
    if (action == Action_Right)
        if (loc == :loc_A)
            copyto!(model,[-1,-1,-1,-1])
        end
    elseif (action == Action_Left)
        if (loc == :loc_B)
            copyto!(model,[-1,-1,-1,-1])
        end
    elseif (action == Action_Suck)
        if (model[1]==1)
            model[3] = 0
        elseif (model[2]==1)
            model[4] = 0
        end
    end

    if (loc == Symbol("loc_A"))
        model[1] = 1; model[2] = 0
        model[3] = (status == Symbol("Dirty")) ? 1 : 0
    else
        model[1] = 0; model[2] = 1
        model[4] = (status == Symbol("Dirty")) ? 1 : 0
    end

    println(model)

    return ModelVacuumState(model)
end

function RunVacuumAgent(AP::Type)
    percept_sequence = [PAC PAD PBD]
    println("Run the $(string(AP))")
    agent=Agent{AP}(AP())
    actions=[]
    for percept in percept_sequence
        action = execute(agent, percept)
        push!(actions, action)
    end
    return actions
end

const RunTableDrivenVacuumAgentResult=[Action_Right, Action_Suck, nothing]
RunTableDrivenVacuumAgent()=RunVacuumAgent(TableDrivenVacuumAgentProgram)
const RunSimpleReflexVacuumAgentResult=[Action_Right, Action_Suck, Action_Suck]
RunSimpleReflexVacuumAgent()=RunVacuumAgent(SimpleReflexVacuumAgentProgram)
const RunReflexVacuumAgentResult=[Action_Right, Action_Suck, Action_Suck]
RunReflexVacuumAgent()=RunVacuumAgent(ReflexVacuumAgentProgram)
RunModelBasedVacuumAgentResult=[Action_Right, Action_Suck, Action_Suck]
RunModelBasedVacuumAgent()=RunVacuumAgent(ModelBasedVacuumAgentProgram)

end


text:
using LoopVectorization, OffsetArrays, Test

@testset "copy" begin

    function copyavx1!(x, y)
        @avx for i ∈ eachindex(x)
            x[i] = y[i]
        end
    end
    function copy_avx1!(x, y)
        @_avx for i ∈ eachindex(x)
            x[i] = y[i]
        end
    end
    function copyavx2!(x, y)
        @avx for i ∈ eachindex(x)
            yᵢ = y[i]
            x[i] = yᵢ
        end
    end
    function copy_avx2!(x, y)
        @_avx for i ∈ eachindex(x)
            yᵢ = y[i]
            x[i] = yᵢ
        end
    end
    function offset_copy!(A, B)
        @inbounds for i=1:size(A,1), j=1:size(B,2)
            A[i,j+2] = B[i,j]
        end
    end
    function offset_copyavx1!(A, B)
        @avx for i=1:size(A,1), j=1:size(B,2)
           A[i,j+2] = B[i,j]
        end
    end
    function offset_copy_avx1!(A, B)
        @_avx for i=1:size(A,1), j=1:size(B,2)
              @inbounds A[i,j+2] = B[i,j]
        end
    end
    function offset_copyavx2!(A, B)
        @avx for i=1:size(A,1), j=1:size(B,2)
            Bᵢⱼ = B[i,j]
            A[i,j+2] = Bᵢⱼ
        end
    end
    function offset_copy_avx2!(A, B)
        @_avx for i=1:size(A,1), j=1:size(B,2)
            Bᵢⱼ = B[i,j]
            A[i,j+2] = Bᵢⱼ
        end
    end
    function make2point3avx!(x)
        a = 1.742416161578685
        b = 1.5
        @avx for i ∈ eachindex(x)
            x[i] = a ^ b
        end
    end
    function make2point3_avx!(x)
        @_avx for i ∈ eachindex(x)
            x[i] = 2.3
        end
    end
    function make23avx!(x)
        @avx for i ∈ eachindex(x)
            @inbounds x[i] = 23
        end
    end
    function make23_avx!(x)
        @_avx for i ∈ eachindex(x)
            x[i] = 23
        end
    end
    function myfillavx!(x, a)
        @avx for i ∈ eachindex(x)
            x[i] = a
        end
    end
    function myfill_avx!(x, a)
        @_avx for i ∈ eachindex(x)
            x[i] = a
        end
    end
    function reversecopy1!(B, A)
        for i in eachindex(A)
            B[i] = A[11-i]
        end
        B
    end
    function reversecopy1avx!(B, A)
        @avx for i in eachindex(A)
            B[i] = A[11-i]
        end
        B
    end
    function reversecopy2!(B, A)
        for i in eachindex(B)
            B[i] = A[-i]
        end
        B
    end
    function reversecopy2avx!(B, A)
        @avx for i in eachindex(B)
            B[i] = A[-i]
        end
        B
    end
    function reversecopy3avx!(B, A)
        n = length(A)
        @avx for i in 1:n
            B[i] = A[n+1-i]
        end
        B
    end
    function copy3!(B, A)
        @assert (length(B) ≥ 3) && (length(A) ≥ 3)
        @avx for i in 1:3
            B[i] = A[i]
        end
        B
    end
    function copyselfdot!(s, x)
        m = zero(eltype(x))
        @avx for i ∈ 1:2
            sᵢ = x[i]
            s[i] = sᵢ
            m += sᵢ * sᵢ
        end
        m
    end

    for T ∈ (Float32, Float64, Int32, Int64)
        @show T, @__LINE__
        R = T <: Integer ? (-T(100):T(100)) : T
        x = rand(R, 237);
        q1 = similar(x); q2 = similar(x);

        fill!(q2, -999999); copyavx1!(q2, x);
        @test x == q2
        fill!(q2, -999999); copy_avx1!(q2, x);
        @test x == q2
        fill!(q2, -999999); copyavx2!(q2, x);
        @test x == q2
        fill!(q2, -999999); copy_avx2!(q2, x);
        @test x == q2
        fill!(q2, -999999); @avx q2 .= x;
        @test x == q2

        B = rand(R, 79, 83);
        A1 = zeros(T, 79, 85);
        A2 = zeros(T, 79, 85);
        offset_copy!(A1, B);
        fill!(A2, 0); offset_copyavx1!(A2, B);
        @test A1 == A2
        fill!(A2, 0); offset_copyavx2!(A2, B);
        @test A1 == A2
        fill!(A2, 0); offset_copy_avx1!(A2, B);
        @test A1 == A2
        fill!(A2, 0); offset_copy_avx2!(A2, B);
        @test A1 == A2

        a = rand(R)
        myfillavx!(x, a);
        fill!(q2, a);
        @test x == q2
        a = rand(R)
        myfill_avx!(x, a);
        fill!(q2, a);
        @test x == q2
        a = rand(R)
        myfill_avx!(x, a);
        fill!(q2, a);
        @test x == q2
        a = rand(R)
        myfillavx!(x, a);
        fill!(q2, a);
        @test x == q2
        q2 .= 23;
        fill!(q1, -99999); make23_avx!(q1);
        @test q2 == q1
        fill!(q1, -99999); make23avx!(q1);
        @test q2 == q1
        if T <: Union{Float32,Float64}
            make2point3avx!(x)
            fill!(q2, 2.3)
            @test x == q2
            fill!(x, -999999); make2point3_avx!(x)
            @test x == q2
        end
        a = rand(R)
        @avx x .= a;
        fill!(q2, a);
        @test x == q2
        a = rand(R)
        @avx x .= a;
        fill!(q2, a);
        @test x == q2

        @test reversecopy1!(zeros(T, 10), collect(1:10)) == reversecopy1avx!(zeros(T, 10), collect(1:10))
        @test reversecopy2!(zeros(T, 10), OffsetArray(collect(1:10), -10:-1)) == reversecopy2avx!(zeros(T, 10), OffsetArray(collect(1:10), -10:-1))
        @test 177:-1:1 == reversecopy3avx!(zeros(T, 177), collect(1:177))

        x = rand(R, 3); y = similar(x);
        @test copy3!(y, x) == x
        fill!(y,0);
        @test copyselfdot!(y, x) ≈ x[1]^2 + x[2]^2
        @test view(x, 1:2) == view(y, 1:2)
    end
end


text:
h5p_set_virtual(dcpl_id, vspace_id, src_file_name, src_dset_name, src_space_id) =
    ccall((:H5Pset_virtual, libhdf5), Herr, (Hid, Hid, Ptr{UInt8}, Ptr{UInt8}, Hid), dcpl_id, vspace_id, src_file_name, src_dset_name, src_space_id)

mutable struct VirtualLayout{S, D, SS}
    shape::S
    dtype::D
    sources::SS
end

VirtualLayout(shape, dtype) = VirtualLayout(tuple(shape...), dtype, [])

repcolon(is, shape) = ifelse.(is .== Colon(), Base.OneTo.(shape), is)

function Base.setindex!(layout::VirtualLayout, v, is...)
    push!(layout.sources, (repcolon(is, layout.shape), v))
    shape = collect(layout.shape)
    for (vis, vs) in layout.sources
        layout.shape = max.(layout.shape, last.(vis))
    end
end

struct VirtualSource{P, N, S, D, I}
    path::P
    name::N
    shape::S
    dtype::D
    is::I
end

function VirtualSource(path, name)
    dset = h5open(path, "r")[name]
    shape, dtype = size(dset), eltype(dset)
    VirtualSource(path, name, shape, dtype, nothing)
end

VirtualSource(dset::HDF5Dataset) =
    VirtualSource(filename(dset), name(dset), size(dset), eltype(dset), nothing)

Base.getindex(vs::VirtualSource, is...) = VirtualSource(vs.path, vs.name, vs.shape, vs.dtype, repcolon(is, vs.shape))

function select_hyperslab(dspace, is)
    isnothing(is) && return dspace
    start = [i[1] - 1 for i in reverse(is)]
    count = [1 for i in is]
    stride = [step(i) for i in is]
    block = [length(i) for i in reverse(is)]
    h5s_select_hyperslab(dspace, H5S_SELECT_SET, start, stride, count, block)
    return dspace
end

function d_create_virtual(parent::Union{HDF5File, HDF5Group}, name::String, layout, virtual = true)
    checkvalid(parent)
    dcpl = p_create(H5P_DATASET_CREATE)
    dspace = dataspace(layout.shape)
    dtype = datatype(layout.dtype)
    if virtual
        for (is, vs) in layout.sources
            src_dspace = dataspace(vs.shape)
            select_hyperslab(dspace, is)
            select_hyperslab(src_dspace, vs.is)
            h5p_set_virtual(dcpl, dspace, vs.path, vs.name, src_dspace)
        end
        d_create(parent, name, dtype, dspace, HDF5Properties(), dcpl)
    else
        d_create(parent, name, dtype, dspace, HDF5Properties(), dcpl)
        for (is, vs) in layout.sources
            vsis = something(vs.is, ntuple(i -> (:), length(vs.shape)))
            parent[name][is...] = h5read(vs.path, vs.name, vsis)
        end
    end
end