text:
"""
    reduce(V::Additive, C::Counts{OrderedCategory}, i, B)

Reduce the instance given by the pair (V, C) to a new instance by giving the
supplied agent, `i`, the supplied bundle, `B`. Returns a reduction, where the
transformation, in addition to converting the allocation to one for the original
instance, allocates `B` to `i`.
"""
function reduce(V::Additive, C::Counts{OrderedCategory}, i, B)
    N, M = agents(V), items(V)
    n′, m′ = na(V) - 1, ni(V) - length(B)

    Vs = zeros(n′, m′)
    λg = Int[]

    for g in M
        g ∈ B && continue

        push!(λg, g)
        g′ = length(λg)

        Vs[1:n′,g′] = [value(V, j, g) for j in N if i != j]
    end

    # Create new ordered categories
    Cs = OrderedCategory[]
    index = 1
    for c in C
        newlength = length(c ∩ λg)
        push!(Cs, OrderedCategory(index, newlength, c.threshold))
        index += newlength
    end

    λi = [j for j in agents(V) if i != j]
    V′, C′ = Additive(Vs), Counts(Cs)

    return Reduction(V′, C′, λi, λg, A -> revert(λg, i, B, A))
end

"""
    reduce(V::Submodular, i, B)

Reduce the instance given by `V` to a new instance by giving the specified
bundle, `B`, to agent `i`. Returns a reduction, where the transformation, in
addition to converting the allocation to one for the original instance,
allocates `B` to `i`.
"""
function reduce(V::Submodular, i, B)

    N, M = agents(V), items(V)

    # Translation table
    λg = [g for g in M if g ∉ B]
    λi = [j for j in N if j != i]

    Vf′ = [B -> value(V, j, Set(λg[g] for g in B)) for j in λi]
    V′ = Submodular(Vf′, length(λg))

    return Reduction(V′, λi, λg, (A) -> revert(λg, i, B, A))

end

"""
    revert(translate::Vector{Int}, i, B, A::Allocation)

Convert an allocation for a reduced instance to one for the original instance,
including giving the removed bundle, `B`, to the removed agent, `i`.
"""
function revert(λg::Vector{Int}, i, B, A::Allocation)
    A′ = Allocation(na(A) + 1, ni(A) + length(B))

    for j in 1:na(A)
        j′ = j + (j >= i)
        give!(A′, j′, [λg[g] for g in bundle(A, j)])
    end

    give!(A′, i, B)

    return A′
end

"""
    reduce(V::Additive, C::Vector{OrderedCategory}, α)

Reduce an ordered instance by normalizing the values and giving any agent that
value an individual item greater than or equal to α the item and any low value
items required to reduce to a valid instance. This reduction is performed
recursively until no more such items exist. The reduction does not decrease the
MMS guarantee of any remaining agents and all agents that are allocated a
bundle in the reduction is guaranteed to value their bundle at least α of their
MMS guarantee.
"""
function reduce(V::Additive, C::Counts{OrderedCategory}, α::Real)
    N, n, M = agents(V), na(V), items(V)

    if n == 1 return reduce(V, C, 1, M) end

    V = normalize(V)

    # Find any agent who values a single item at least α and give that item to
    # the agent along with the least valuable items in the remaining bundles
    # that must be given to the agent to guarantee a valid instance.
    for i in N, g in M
        if value(V, i, g) >= α
            bundle = union(
                Set{Int}(g),
                [c[end - required(c, n) + (g in c) + 1:end] for c in C]...)

            R = reduce(V, C, i, bundle)

            V, C = profile(R), constraint(R)
            # Recursive application, as the removed items may cause items worth
            # less than α to be worth α or more after a new normalization.
            R′ = reduce(V, C, α)

            # Combine the reductions
            return chain(R, R′)
        end
    end

    # Return an empty reduction
    return Reduction(V, C)
end

"""
    reduce(V::Profile, α::Real)

Produce a reduced instance by giving an item to any agent that values it at `α`
or more. This reduction is performed repeatedly, until no such item exists.
"""
function reduce(V::Profile, α::Real)
    N, M = agents(V), items(V)

    for i in N, g in M
        if value(V, i, g) ≥ α
            R = reduce(V, i, Set(g))
            V = profile(R)

            # Recursive application on the new instance
            R′ = reduce(V, α)

            # Combine the reductions
            return chain(R, R′)
        end
    end

    return Reduction(V)
end

"""
    order(V::Additive, C::Counts)

Create an ordered instance for the given weights and categories. The items are
reorded such that each category has a continous range of indices for its items.
Returns a reduction, with a transformation that converts an allocation to one in
the original instance where each agent gets at least the same value as in the
ordered instance.
"""
function order(V::Additive, C::Counts)
    N = agents(V)
    Vo = zeros(na(V), ni(V))
    Co = OrderedCategory[]

    m = 1
    for c in C
        for i in N
            Vo[i,m:m+length(c)-1] = sort([value(V, i, g) for g in c], rev=true)
        end

        push!(Co, OrderedCategory(m, length(c), threshold(c)))
        m += length(c)
    end

    V′, C′ = Additive(Vo), Counts(Co)

    return Reduction(V′, C′, agents(V), items(V), A -> revert(V, C, C′, A))
end

"""
    revert(V::Additive, C::Counts, C′::Counts, A)

Convert an allocation for the ordered instance (`C′`) to one for the original
instance `(V, C)`.
"""
function revert(V::Additive, C::Counts, C′::Counts, A)
    A′ = Allocation(na(A), ni(A))

    for (unordered_category, ordered_category) in zip(C, C′)
        items = [g for g in unordered_category]
        for g in ordered_category
            i = owner(A, g)

            # Select the good in the category with the highest value
            g′ = maximum((g) -> (value(V, i, g), g), items)[2]

            give!(A′, i, g′)
            items = setdiff(items, g′)
        end
    end

    return A′
end


text:
"""
FEniCS tutorial demo program: Poisson equation with Dirichlet conditions.
Test problem is chosen to give an exact solution at all nodes of the mesh.
  -Laplace(u) = f    in the unit square
            u = u_D  on the boundary
  u = 1 + x^2 + 2y^2 = u_D
  f = -6
This is an extended version of the demo program poisson.py which
encapsulates the solver as a Python function.
"""
module ft12

using FenicsPy

function solver(f, u_D, Nx, Ny, degree=1)
    """
    Solve -Laplace(u) = f on [0,1] x [0,1] with 2*Nx*Ny Lagrange
    elements of specified degree and u = u_D (Expresssion) on
    the boundary.
    """

    # Create mesh and define function space
    mesh = UnitSquareMesh(Nx, Ny)
    V = FunctionSpace(mesh, "P", degree)

    bc = DirichletBC(V, u_D, "on_boundary")

    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    a = dot(grad(u), grad(v))*dx
    L = f*v*dx

    # Compute solution
    u = FeFunction(V)
    solve(a == L, u, bc)

    return u
end

function run_solver()
    "Run solver to compute and post-process solution"

    # Set up problem parameters and call solver
    u_D = Expression("1 + x[0]*x[0] + 2*x[1]*x[1]", degree=2)
    f = Constant(-6.0)
    u = solver(f, u_D, 8, 8, 1)

    # Plot solution and mesh
    #plot(u)
    #plot(u.function_space().mesh())

    # Save solution to file in VTK format
    vtkfile = File("poisson_solver/solution.pvd")
    vtkfile << u
end

function test_solver()
    "Test solver by reproducing u = 1 + x^2 + 2y^2"

    # Set up parameters for testing
    tol = 1E-10
    u_D = Expression("1 + x[0]*x[0] + 2*x[1]*x[1]", degree=2)
    f = Constant(-6.0)

    # Iterate over mesh sizes and degrees
    for Nxy in [(3, 3), (3, 5), (5, 3), (20, 20)]
        Nx, Ny = Nxy
        for degree in [1, 2, 3]
            println("Solving on a 2 x (", Nx," x ", Ny, ") mesh with P", degree, " elements.")

            # Compute solution
            u = solver(f, u_D, Nx, Ny, degree)

            # Extract the mesh
            mesh = u.function_space().mesh()

            # Compute maximum error at vertices
            vertex_values_u_D = u_D.compute_vertex_values(mesh)
            vertex_values_u  = u.compute_vertex_values(mesh)
            error_max = max(abs.(vertex_values_u_D - vertex_values_u)...)

            # Check maximum error
            msg = string("error_max = " , error_max)
            @assert error_max<tol msg
        end
    end
end

run_solver()
println("------------------")
test_solver()

end # module ft12

ext: jl


text:
## Handling of the topological archive internally used to reckognize topologies.
import Serialization
import Pkg

const CRYSTAL_NETS_VERSION = VersionNumber(Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"])
const arc_location = normpath(joinpath(@__DIR__, "..", "archives/"))
const CRYSTAL_NETS_ARCHIVE = if isdir(arc_location) && !isempty(readdir(arc_location))
    flag, parsed = parse_arcs(arc_location)
    if !flag
        error("""CrystalNets.jl appears to have a broken installation (the archive version is older than that package's).
        Please rebuild CrystalNets.jl with `import Pkg; Pkg.build("CrystalNets")`.
        """)
    end
    parsed
else
    error("""CrystalNets.jl appears to have a broken installation (missing default archive).
          Please rebuild CrystalNets.jl with `import Pkg; Pkg.build("CrystalNets")`.

          This issue may be due to an error while using `CrystalNets.clean_default_archive!(...)`.
          If you have not encountered any such error and you did not modify or erase the archive located at $arc_location please open an issue at https://github.com/coudertlab/CrystalNets.jl/issues/new
          """)
end
const REVERSE_CRYSTAL_NETS_ARCHIVE = Dict{String,String}(last(x) => first(x) for x in CRYSTAL_NETS_ARCHIVE)

function validate_archive(custom_arc)::Dict{String,String}
    arc = try
        Serialization.deserialize(custom_arc)
    catch
        nothing
    end
    if arc isa Dict{String,String} && !isempty(arc) && isnumeric(first(first(keys(arc))))
        @info "Processing input as a serialized CrystalNets.jl archive"
    else
        try
            flag, parsed = parse_arc(custom_arc)
            if flag
                arc = parsed
            else
                @info "Processing input as a topological .arc file"
                if occursin("Made by CrystalNets.jl", readline(custom_arc))
                    @info "This archive was generated by an older version of CrystalNets."
                end
                @info "Keys will be converted to the topological genome used by CrystalNets. This may take a while."
                arc = Dict{String,String}()
                Threads.@threads for (key, id) in collect(parsed)
                    PeriodicGraph(key) isa PeriodicGraph3D || continue # TODO should eventually be removed
                    genome::String = try
                        string(topological_genome(PeriodicGraph3D(key)))
                    catch e
                        if e isa InterruptException ||
                          (e isa TaskFailedException && e.task.result isa InterruptException)
                            rethrow()
                        end
                        println(stderr, "Failed for net "*id*" with error:")
                        Base.display_error(stderr, e, catch_backtrace())
                        println(stderr)
                        ""
                    end
                    if !isempty(genome)
                        arc[genome] = id
                    end
                end
            end
        catch e
            print(stderr, """
            Impossible to parse input as a topological archive. Please use a format
            similar to that of the RCSR Systre .arc file, with for each entry at
            least the "id" and "key" fields.

            This error may also occur if the given archive was empty. If you
            wish to set an empty archive, use `CrystalNets.empty_default_archive!()`
            """, '\n', "Encountered error: ")
            rethrow()
            # showerror(stderr, e)
            # Base.display_error(stderr, e, catch_backtrace())
        end
    end
    return arc
end

"""
    clean_default_archive!(custom_arc=nothing; validate=true, refresh=true)

Erase the default archive used by CrystalNets.jl to reckognize known topologies.

If no parameter is provided, the new default archive is the RCSR Systre archive.
Otherwise it is replaced by the file located at `custom_arc`. In this case, the
`validate` parameter controls whether the new file is checked and converted to a
format usable by CrystalNets.jl. If unsure, leave it set.

The `refresh` optional parameter controls whether the current archive should be
replaced by the new default one.

!!! warning
    This archive will be kept and used for subsequent runs of CrystalNets.jl, even
    if you restart your Julia session.

    To only change the archive for the current session, use `CrystalNets.change_current_archive!(custom_arc)`.

    See also `refresh_current_archive!` for similar uses.

!!! warning
    Using an invalid archive will make CrystalNets.jl unusable. If this happens,
    simply run `CrystalNets.clean_default_archive!()` to revert to the RCSR Systre
    default archive.

See also `add_to_current_archive!`, `change_current_archive!`, `refresh_current_archive!`,
`set_default_archive!`, `empty_default_archive!`
"""
function clean_default_archive!(custom_arc=nothing; validate=true, refresh=true)
    rm(arc_location; force=true)
    if custom_arc isa Nothing
        cp(joinpath(@__DIR__, "..", "deps", "RCSR.arc"), arc_location)
    else
        if validate
            arc = validate_archive(custom_arc)
            export_arc(arc_location, false, arc)
        else
            cp(custom_arc, arc_location)
        end
    end
    if refresh
        refresh_current_archive!()
    end
    nothing
end

"""
    set_default_archive!()

Set the current archive as the new default archive.

!!! warning
    This archive will be kept and used for subsequent runs of CrystalNets.jl, even
    if you restart your Julia session.

See also `add_to_current_archive!`, `change_current_archive!`, `refresh_current_archive!`,
`clean_default_archive!`, `empty_default_archive!`
"""
function set_default_archive!()
    global CRYSTAL_NETS_ARCHIVE
    export_arc(arc_location)
end

"""
    empty_default_archive!(; refresh=true)

Empty the default archive. This will prevent CrystalNets from reckognizing any
topology before they are explicitly added.

The `refresh` optional parameter controls whether the current archive should also
be emptied.

!!! warning
    This empty archive will be kept and used for subsequent runs of CrystalNets.jl, even
    if you restart your Julia session. If you only want to empty the current archive,
    use `empty!(CrystalNets.CRYSTAL_NETS_ARCHIVE)`.

See also `add_to_current_archive!`, `change_current_archive!`, `refresh_current_archive!`,
`clean_default_archive!`, `set_default_archive!`
"""
function empty_default_archive!(; refresh=true)
    global CRYSTAL_NETS_ARCHIVE
    global REVERSE_CRYSTAL_NETS_ARCHIVE
    export_arc(arc_location, true)
    if refresh
        empty!(CRYSTAL_NETS_ARCHIVE)
        empty!(REVERSE_CRYSTAL_NETS_ARCHIVE)
    end
    nothing
end

function _change_current_archive!(newarc)
    global CRYSTAL_NETS_ARCHIVE
    global REVERSE_CRYSTAL_NETS_ARCHIVE
    empty!(CRYSTAL_NETS_ARCHIVE)
    empty!(REVERSE_CRYSTAL_NETS_ARCHIVE)
    merge!(CRYSTAL_NETS_ARCHIVE, newarc)
    merge!(REVERSE_CRYSTAL_NETS_ARCHIVE,
           Dict{String,String}(last(x) => first(x) for x in CRYSTAL_NETS_ARCHIVE))
    nothing
end

"""
    change_current_archive!(custom_arc; validate=true)

Erase the current archive used by CrystalNets.jl to reckognize known topologies and
replace it with the archive stored in the file located at `custom_arc`.

The `validate` optional parameter controls whether the new file is checked and converted
to a format usable by CrystalNets.jl. If unsure, leave it set.

!!! note
    This modification will only last for the duration of this Julia session.

    If you wish to change the default archive and use it for subsequent runs, use
    `CrystalNets.clean_default_archive!`.

!!! warning
    Using an invalid archive will make CrystalNets.jl unusable. If this happens,
    simply run `CrystalNets.refresh_current_archive!()` to revert to the
    default archive.

See also `add_to_current_archive!`, `refresh_current_archive!`, `clean_default_archive!`,
`set_default_archive!`, `empty_default_archive!`
"""
function change_current_archive!(custom_arc; validate=true)
    arc::Dict{String,String} = if validate
         validate_archive(custom_arc)
    else
        last(parse_arc(custom_arc))
    end
    _change_current_archive!(arc)
end

"""
    refresh_current_archive!()

Revert the current topological archive to the default one.

See also `add_to_current_archive!`, `change_current_archive!`, `clean_default_archive!`,
`set_default_archive!`, `empty_default_archive!`
"""
function refresh_current_archive!()
    _change_current_archive!(last(parse_arc(arc_location)))
end

function _update_archive!(id, genome)
    global CRYSTAL_NETS_ARCHIVE
    global REVERSE_CRYSTAL_NETS_ARCHIVE
    CRYSTAL_NETS_ARCHIVE[genome] = id
    REVERSE_CRYSTAL_NETS_ARCHIVE[id] = genome
    nothing
end

"""
    add_to_current_archive!(id, genome)

Mark `genome` as the topological genome associated with the name `id` in the
current archive.

The input `id` and `genome` are not modified by this operation.

!!! note
    This modification will only last for the duration of this Julia session.

    If you wish to save the archive and use it for subsequent runs, use
    `CrystalNets.set_default_archive!` after calling this function.

See also `change_current_archive!`, `refresh_current_archive!`, `clean_default_archive!`,
`set_default_archive!`, `empty_default_archive!`
"""
function add_to_current_archive!(id::AbstractString, genome::AbstractString)
    if !isnumeric(first(genome))
        throw(ArgumentError("""
            This genome ("$genome") does not look like a genome. Are you sure you did not mix `id` and `genome`?

            If you really want to associate this id with this genome, use `CrystalNets._update_archive!(id, genome)`
            """))
    end
    global CRYSTAL_NETS_ARCHIVE
    for (x,y) in CRYSTAL_NETS_ARCHIVE
        if x == genome
            y == id && return
            throw(ArgumentError("""
                This genome is already registered under the name "$y".

                If you really want to change the name associated with it, use `CrystalNets._update_archive!(id, genome)`
                """))
        end
        if y == id
            throw(ArgumentError("""
                The name $id already corresponds to a different genome: "$x"

                If you really want to store another genome with the same name, use `CrystalNets._update_archive!(id, genome)`
                """))
        end
    end
    _update_archive!(id, genome)
end

"""
    make_archive(path, destination=nothing)

Make an archive from the files located in the directory given by `path` and export
it to `destination`, if specified.

The archive can then be used with `change_current_archive!(destination; validate=false)`
for instance.
"""
function make_archive(path, destination)
    arc = Dict{String,String}()
    Threads.@threads for f in readdir(path; join=false)
        name = splitext(f)[1]
        print("Handling "*name*"... ")
        genome::String = try
            x = topological_genome(CrystalNet(parse_chemfile(path*f)))
            println(name*" done.")
            string(x)
        catch e
            if e isa InterruptException ||
              (e isa TaskFailedException && e.task.result isa InterruptException)
                rethrow()
            end
            println(stderr, "Failed for file "*f*" with error:")
            showerror(stderr, e)
            Base.display_error(stderr, e, catch_backtrace())
            println(stderr)
            "DISREGARD"
        end
        arc[string(genome)] = name
    end
    delete!(arc, "DISREGARD")
    if !(destination isa Nothing)
        export_arc(destination, false, arc)
    end
    return arc
end


text:
"""
    LazyModel{Name, T} <: AbstractMillModel

A model node for processing [`LazyNode`](@ref)s. It applies a (sub)model `m` stored in it to data of
the [`LazyNode`](@ref) after calling [`Mill.unpack2mill`](@ref).

# Examples
```jldoctest unpack2mill; output=false
function Mill.unpack2mill(ds::LazyNode{:Sentence})
    s = split.(ds.data, " ")
    x = NGramMatrix(reduce(vcat, s))
    BagNode(ArrayNode(x), Mill.length2bags(length.(s)))
end
# output
```

```jldoctest unpack2mill
julia> Random.seed!(0);

julia> n = LazyNode{:Sentence}(["foo bar", "baz"])
LazyNode{Sentence} 	# 2 obs, 90 bytes

julia> m = LazyModel{:Sentence}(BagModel(Dense(2053, 3), SegmentedMean(3), identity))
LazyModel{Sentence}
  └── BagModel ↦ SegmentedMean(3) ↦ ArrayModel(identity) 	# 1 arrays, 3 params (all zero), 52 bytes
        └── ArrayModel(Dense(2053, 3)) 	# 2 arrays, 6_162 params, 24.148 KiB
```

```jldoctest unpack2mill; filter=$(DOCTEST_FILTER)
julia> m(n)
3×2 ArrayNode{Matrix{Float32}, Nothing}:
 -0.006...  -0.022...
  0.034...   0.055...
 -0.062...   0.071...
```

See also: [`AbstractMillModel`](@ref), [`LazyNode`](@ref), [`Mill.unpack2mill`](@ref).
"""
struct LazyModel{Name, T <: AbstractMillModel} <: AbstractMillModel
    m::T
end

"""
    LazyModel([Name::Symbol], m::AbstractMillModel)
    LazyModel{Name}(m::AbstractMillModel)

Construct a new [`LazyModel`](@ref) with name `Name`, and model `m`.

# Examples
```jldoctest
julia> LazyModel{:Sentence}(ArrayModel(Dense(2, 2)))
LazyModel{Sentence}
  └── ArrayModel(Dense(2, 2)) 	# 2 arrays, 6 params, 104 bytes
```

See also: [`AbstractMillModel`](@ref), [`LazyNode`](@ref), [`Mill.unpack2mill`](@ref).
"""
LazyModel(Name::Symbol, m::T) where T <: AbstractMillModel = LazyNode{Name, T}(m)
LazyModel{Name}(m::M) where {Name, M} = LazyModel{Name, M}(m)

Flux.@functor LazyModel

(m::LazyModel{Name})(x::LazyNode{Name}) where {Name} = m.m(unpack2mill(x))

function HiddenLayerModel(m::LazyModel{N}, ds::LazyNode{N}, n) where {N}
    hm, o = HiddenLayerModel(m.m, unpack2mill(ds), n)
    return(LazyModel{N}(hm), o )
end

function mapactivations(hm::LazyModel{N}, x::LazyNode{N}, m::LazyModel{N}) where {N}
    ho, o = mapactivations(hm.m, unpack2mill(x), m.m)
end

# Base.hash(m::LazyModel{T}, h::UInt) where {T} = hash((T, m.m), h)
# (m1::LazyModel{T} == m2::LazyModel{T}) where {T} = m1.m == m2.m


text:
using FileIO
using LazIO
using LasIO
using Test
using Tables

workdir = @__DIR__
lasio_testdir = joinpath(dirname(pathof(LasIO)), "..", "test")
# source: http://www.liblas.org/samples/
fname_las = "libLAS_1.2.las" # point format 0
testfile_str = joinpath(lasio_testdir, fname_las)
laz_out = joinpath(workdir, "out.laz")
testfile = File{format"LAZ_"}(testfile_str)

header, pointdata_all = load(testfile)
@testset "Header" begin
    @test header.version_major == 1
    @test header.version_minor == 2
    @test header.data_format_id == 0
    @test header.records_count == 497536
    @test header.x_scale == 0.01
    @test header.y_max == 379999.99
    @test header.header_size == 227
    @test LazIO.boundingbox(header) == (xmin = 1.44e6, ymin = 375000.03, zmin = 832.1800000000001, xmax = 1.44499996e6, ymax = 379999.99, zmax = 972.6700000000001)
end

@testset "Range indexing" begin
    _, pointdata_276 = load(testfile, range=276)
    _, pointdata_array = load(testfile, range=[1,276,277,497536])
    _, pointdata_colon = load(testfile, range=:)

    @test pointdata_all[1] == pointdata_colon[1] == pointdata_array[1]
    @test pointdata_276[1] == pointdata_colon[276] == pointdata_array[2]
    @test pointdata_all[end] == pointdata_colon[end] == pointdata_array[end]
end

@testset "LazDataSet" begin
    ds = LazIO.open(testfile_str)
    @test ds isa LazIO.LazDataset
    @test first(ds) isa LazIO.LazPoint
    @inferred first(ds)
    @test LazIO.boundingbox(ds) == (xmin = 1.44e6, ymin = 375000.03, zmin = 832.1800000000001, xmax = 1.44499996e6, ymax = 379999.99, zmax = 972.6700000000001)
    close(ds)
end

@testset "Table interface" begin
    ds = LazIO.open(testfile_str)
    @test Tables.istable(LazIO.LazDataset)
    @test Tables.rowaccess(LazIO.LazDataset)
    @test first(ds) isa LazIO.LazPoint
    @inferred first(Tables.rows(ds))
    close(ds)
end

@testset "Writing" begin
    # open reader
    reader = Ref{Ptr{Cvoid}}(C_NULL)
    LazIO.@check reader[] LazIO.laszip_create(reader)
    LazIO.@check reader[] LazIO.laszip_open_reader(reader[], testfile_str, Ref(Cint(0)))

    # create writer
    writer = Ref{Ptr{Cvoid}}(C_NULL)
    LazIO.@check writer[] LazIO.laszip_create(writer)

    # copy header from reader to writer (to be updated later)
    header_ptr = Ref{Ptr{LazIO.LazHeader}}(C_NULL)
    LazIO.@check header_ptr[] LazIO.laszip_get_header_pointer(reader[], header_ptr)
    LazIO.@check writer[] LazIO.laszip_set_header(writer[], header_ptr[])

    # open writer
    LazIO.@check writer[] LazIO.laszip_open_writer(writer[], laz_out, Cint(1))

    # read the first point
    point_ptr = Ref{Ptr{LazIO.LazPoint}}(C_NULL)
    LazIO.@check point_ptr[] LazIO.laszip_get_point_pointer(reader[], point_ptr)
    LazIO.@check reader[] LazIO.laszip_read_point(reader[])

    # modify point
    p = unsafe_load(point_ptr[])
    p.classification = LasIO.raw_classification(
        0x05,
        Bool((p.classification & 0b00100000) >> 5),
        Bool((p.classification & 0b01000000) >> 6),
        Bool((p.classification & 0b10000000) >> 7),
    )

    # copy the point
    LazIO.@check writer[] LazIO.laszip_set_point(writer[], Ref(p))

    # write it to the writer
    LazIO.@check writer[] LazIO.laszip_write_point(writer[])

    # update the inventory
    LazIO.@check writer[] LazIO.laszip_update_inventory(writer[])

    # close files and destroy pointers
    LazIO.@check reader[] LazIO.laszip_close_reader(reader[])
    LazIO.@check writer[] LazIO.laszip_close_writer(writer[])
    LazIO.@check reader[] LazIO.laszip_destroy(reader[])
    LazIO.@check writer[] LazIO.laszip_destroy(writer[])

    # check the resulting LAZ file
    ds = LazIO.open(laz_out)
    @test ds.header.number_of_point_records === UInt32(1)
    @test length(ds) === 1
    @test first(ds).classification === 0x05
    close(ds)
    rm(laz_out)
end

@testset "Stream output (dataset)" begin
    ds = LazIO.open(testfile_str)
    laz_stream_out = joinpath(workdir, "out-stream.laz")

    LazIO.write(laz_stream_out, ds) do io
        for p in ds
            LazIO.writepoint(io, p)
        end
    end

    ds_out = LazIO.open(laz_stream_out)

    # TODO write a method such that we can check point equality
    p1 = first(ds)
    p2 = first(ds_out)
    header1 = ds.header
    header2 = ds_out.header
    @test p1.X === p2.X
    @test p1.Y === p2.Y
    @test p1.Z === p2.Z
    @test p1.intensity === p2.intensity
    @test header1.number_of_point_records === header2.number_of_point_records

    close(ds)
    close(ds_out)
    rm(laz_stream_out)
end

@testset "Stream output (header)" begin
    ds = LazIO.open(testfile_str)
    laz_stream_out = joinpath(workdir, "out-stream2.laz")

    LazIO.write(laz_stream_out, ds.header) do io
        for p in ds
            LazIO.writepoint(io, p)
        end
    end

    ds_out = LazIO.open(laz_stream_out)

    # TODO write a method such that we can check point equality
    p1 = first(ds)
    p2 = first(ds_out)
    header1 = ds.header
    header2 = ds_out.header
    @test p1.X === p2.X
    @test p1.Y === p2.Y
    @test p1.Z === p2.Z
    @test p1.intensity === p2.intensity
    @test header1.number_of_point_records === header2.number_of_point_records

    close(ds)
    close(ds_out)
    rm(laz_stream_out)
end


text:
 # Generate GCSPET instance according to Guo

 import Distributions: DiscreteUniform

 """
 generate_location_constrained_jobs(njob, D_l, maxperl)
Sample `njob` job locations from distribution D_l, with the constraint that at most `maxperl` jobs can occur at one location. Returns a jobs-per-location-dict of type `Dict{T,Vector{Int}}`, where `T` is the eltype of `D_l`.
"""
function generate_location_constrained_jobs(njob, Dₗ, maxperl; base_id = 1)
    id = 0
    posjobmap = Dict{eltype(Dₗ), Vector{Int}}()
    while id < njob
        l = rand(Dₗ)
        !haskey(posjobmap, l) && (posjobmap[l] = Int[])
        if length(posjobmap[l]) < maxperl
            push!(posjobmap[l], base_id + id)
            id += 1
        end
    end
    return posjobmap
end

"""
 generate_crane_starting_positions(njob, ncrane)
Generate `ncrane` starting positions at random locations but with fixed intervals between the positions.
"""
function generate_crane_starting_positions(njob, ncrane)
    step = Int(round(njob/ncrane))
    lq_left = rand(1:step-1)
    return [lq_left + (c-1) * step for c in 1:ncrane]
end

"""
 generate_jobs(njobs, ncrane, load)
Generate the data for a set of jobs for the GCSPET, according tot the distributions presented by Guo et al.
"""
function generate_jobdata(njob, ncrane, load)
    ntruck = round(Int, load * njob)
    ntrain = njob - ntruck
    Dₚ= DiscreteUniform(3,10)                                       # Generate processing times
    P = rand(Dₚ, njob)
    Dₐ = DiscreteUniform(1, sum(P) ÷ ncrane)                        # Generate arrival times
    A = vcat(zeros(Int, ntrain), rand(Dₐ, ntruck))
    Dₗ = DiscreteUniform(1,njob)                                    # Generate locations
    JT = vcat(ones(Int, ntrain), zeros(Int, ntruck))
    trainjobs = generate_location_constrained_jobs(ntrain, Dₗ, 2)
    truckjobs = generate_location_constrained_jobs(ntruck, Dₗ, 2; base_id = ntrain+1)
    L = zeros(Int, njob)
    MT = zeros(Int, njob)
    encode_jobtypes!(MT, L, trainjobs)
    encode_jobtypes!(MT, L, truckjobs)
    ID = 1:njob
    return (ID, L, A, P, JT, MT)
end

generate_jobs(njobs, ncranes, load) = map(Job, generate_jobdata(njobs, ncranes, load))

function generate_instance(njobs, ncranes, load; sample_id = 1, safety = 1, speed = 1)
    jobdata = generate_jobdata(njobs, ncranes, load)
    Ω = map(Job, jobdata...)
    lQ = generate_crane_starting_positions(njobs, ncranes)
    Q = map(a -> Crane(a..., speed, zone(first(a), ncranes, njobs, safety), safety), Iterators.zip(1:ncranes, lQ))
    return Instance(toname(Ω, lQ, s_id = sample_id), Ω, Q)
end

function encode_jobtypes!(MT, L, posjobmap)
    for (l, jobs) in posjobmap
        for (i,j) in enumerate(jobs)
            L[j] = l
            MT[j] = i-1
        end
    end
end

toname(njobs::Int, ncranes::Int, load::Float64) = join((njobs, ncranes, load), "_") * ".dat"

toname(Ω, lQ; s_id = 1) = join((length(Ω), length(lQ), count(istruck, Ω)/length(Ω), s_id), "_")


text:
module IntCode

using AdventOfCode2019

function run_program!(data::Array{T, 1}, input::Union{Channel{T},Nothing}, output::Union{Channel{T},Nothing}, done::Union{Channel{Bool},Nothing}; waitingForInput::Union{Channel{Bool},Nothing} = nothing, preserve_program = false, reboot = false) where T <: Integer
    out = Array{T,1}()
    @label start
    extMemory = Dict{T,T}()

    i = T(1)  # instruction pointer
    relativeBase = T(0)  # relative base
    digs = Array{UInt8,1}(undef, 5)
    modes = Array{UInt8,1}(undef, 3)

    params = Array{T,1}(undef, 3)
    while true
        fill!(digs, zero(UInt8))
        digits!(digs, data[i])
        optcode = digs[1] + 10digs[2]

        # modes for parameters 1, 2 and 3
        # mode 0: position mode
        # mode 1: immediate mode
        # mode 2: relative mode
        for j = 1:3
            modes[j] = digs[3+j-1]
        end

        if optcode == 99
            break
        elseif optcode == 1 || optcode == 2  # addition and multiplication
            for j = 1:3
                params[j] = _parameter(data, i + j, relativeBase, modes[j])
            end
            if optcode == 1
                val = _get(data, extMemory, params[1]; preserve = preserve_program) + _get(data, extMemory, params[2]; preserve = preserve_program)
            else
                val = _get(data, extMemory, params[1]; preserve = preserve_program) * _get(data, extMemory, params[2]; preserve = preserve_program)
            end
            _set!(data, extMemory, params[3], val; preserve = preserve_program)
            i += 4
        elseif optcode == 3  # read input
            param = _parameter(data, i + 1, relativeBase, modes[1])
            if input != nothing
                if waitingForInput != nothing
                    put!(waitingForInput, true)
                end
                _set!(data, extMemory, param, take!(input); preserve = preserve_program)
                if waitingForInput != nothing
                    take!(waitingForInput)
                end
            end
            i += 2
        elseif optcode == 4  # output
            param = _parameter(data, i + 1, relativeBase, modes[1])
            value = _get(data, extMemory, param; preserve = preserve_program)
            if output == nothing
                push!(out, value)
            else
                put!(output, value)
            end
            i += 2
        elseif optcode == 5 || optcode == 6 # jump-if-true and jump-if-false
            par1 = _parameter(data, i + 1, relativeBase, modes[1])
            op = optcode == 5 ? eval(!=) : eval(==)
            if op(_get(data, extMemory, par1; preserve = preserve_program), 0)
                par2 = _parameter(data, i + 2, relativeBase, modes[2])
                i = _get(data, extMemory, par2; preserve = preserve_program) + 1
            else
                i += 3
            end
        elseif optcode == 7 || optcode == 8  # less than and equals
            for j = 1:3
                params[j] = _parameter(data, i + j, relativeBase, modes[j])
            end
            if optcode == 7
                value = (_get(data, extMemory, params[1]; preserve = preserve_program) < _get(data, extMemory, params[2]; preserve = preserve_program))
            else
                value = (_get(data, extMemory, params[1]; preserve = preserve_program) == _get(data, extMemory, params[2]; preserve = preserve_program))
            end
            _set!(data, extMemory, params[3], value ? T(1) : T(0); preserve = preserve_program)
            i += 4
        elseif optcode == 9  # adjust relative base
            param = _parameter(data, i + 1, relativeBase, modes[1])
            relativeBase += _get(data, extMemory, param; preserve = preserve_program)
            i += 2
        else
            throw(AssertionError("Invalid optcode: $optcode"))
        end
    end
    if done != nothing
        put!(done, true)
    end
    if reboot
        @goto start
    end
    return out
end

@inline function _get(data::Array{T,1}, memory::Dict{T,T}, index::T; preserve = false) where T <: Integer
    preserve && haskey(memory, index) && return memory[index]
    index > length(data) && return T(0)
    return data[index]
end

@inline function _set!(data::Array{T,1}, memory::Dict{T,T}, index::T, value::T; preserve = false) where T <: Integer
    if preserve
        memory[index] = value
        return
    end
    n = length(data)
    if index > n
        resize!(data, index)
        data[n+1:index] .= T(0)
    end
    data[index] = value
end

@inline function _parameter(data::Array{T,1}, index::T, relativeBase::T, mode::UInt8) where T <: Integer
    if mode == 0
        return data[index] + 1
    elseif mode == 1
        return index
    elseif mode == 2
        return data[index] + relativeBase + 1
    end
    throw(AssertionError("Invalid mode $mode"))
end

function run_program!(data::Array{T,1}, input::Array{T,1}) where T <: Integer
    c = Channel{Int}(length(input)+1)
    for value in input
        put!(c, value)
    end
    out = run_program!(data, c, nothing, nothing)
    close(c)
    return out
end

end # module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewayRedirectConfigurationPropertiesFormat <: SwaggerModel
    redirectType::Any # spec type: Union{ Nothing, RedirectTypeEnum } # spec name: redirectType
    targetListener::Any # spec type: Union{ Nothing, SubResource } # spec name: targetListener
    targetUrl::Any # spec type: Union{ Nothing, String } # spec name: targetUrl
    includePath::Any # spec type: Union{ Nothing, Bool } # spec name: includePath
    includeQueryString::Any # spec type: Union{ Nothing, Bool } # spec name: includeQueryString
    requestRoutingRules::Any # spec type: Union{ Nothing, Vector{SubResource} } # spec name: requestRoutingRules
    urlPathMaps::Any # spec type: Union{ Nothing, Vector{SubResource} } # spec name: urlPathMaps
    pathRules::Any # spec type: Union{ Nothing, Vector{SubResource} } # spec name: pathRules

    function ApplicationGatewayRedirectConfigurationPropertiesFormat(;redirectType=nothing, targetListener=nothing, targetUrl=nothing, includePath=nothing, includeQueryString=nothing, requestRoutingRules=nothing, urlPathMaps=nothing, pathRules=nothing)
        o = new()
        validate_property(ApplicationGatewayRedirectConfigurationPropertiesFormat, Symbol("redirectType"), redirectType)
        setfield!(o, Symbol("redirectType"), redirectType)
        validate_property(ApplicationGatewayRedirectConfigurationPropertiesFormat, Symbol("targetListener"), targetListener)
        setfield!(o, Symbol("targetListener"), targetListener)
        validate_property(ApplicationGatewayRedirectConfigurationPropertiesFormat, Symbol("targetUrl"), targetUrl)
        setfield!(o, Symbol("targetUrl"), targetUrl)
        validate_property(ApplicationGatewayRedirectConfigurationPropertiesFormat, Symbol("includePath"), includePath)
        setfield!(o, Symbol("includePath"), includePath)
        validate_property(ApplicationGatewayRedirectConfigurationPropertiesFormat, Symbol("includeQueryString"), includeQueryString)
        setfield!(o, Symbol("includeQueryString"), includeQueryString)
        validate_property(ApplicationGatewayRedirectConfigurationPropertiesFormat, Symbol("requestRoutingRules"), requestRoutingRules)
        setfield!(o, Symbol("requestRoutingRules"), requestRoutingRules)
        validate_property(ApplicationGatewayRedirectConfigurationPropertiesFormat, Symbol("urlPathMaps"), urlPathMaps)
        setfield!(o, Symbol("urlPathMaps"), urlPathMaps)
        validate_property(ApplicationGatewayRedirectConfigurationPropertiesFormat, Symbol("pathRules"), pathRules)
        setfield!(o, Symbol("pathRules"), pathRules)
        o
    end
end # type ApplicationGatewayRedirectConfigurationPropertiesFormat

const _property_map_ApplicationGatewayRedirectConfigurationPropertiesFormat = Dict{Symbol,Symbol}(Symbol("redirectType")=>Symbol("redirectType"), Symbol("targetListener")=>Symbol("targetListener"), Symbol("targetUrl")=>Symbol("targetUrl"), Symbol("includePath")=>Symbol("includePath"), Symbol("includeQueryString")=>Symbol("includeQueryString"), Symbol("requestRoutingRules")=>Symbol("requestRoutingRules"), Symbol("urlPathMaps")=>Symbol("urlPathMaps"), Symbol("pathRules")=>Symbol("pathRules"))
const _property_types_ApplicationGatewayRedirectConfigurationPropertiesFormat = Dict{Symbol,String}(Symbol("redirectType")=>"RedirectTypeEnum", Symbol("targetListener")=>"SubResource", Symbol("targetUrl")=>"String", Symbol("includePath")=>"Bool", Symbol("includeQueryString")=>"Bool", Symbol("requestRoutingRules")=>"Vector{SubResource}", Symbol("urlPathMaps")=>"Vector{SubResource}", Symbol("pathRules")=>"Vector{SubResource}")
Base.propertynames(::Type{ ApplicationGatewayRedirectConfigurationPropertiesFormat }) = collect(keys(_property_map_ApplicationGatewayRedirectConfigurationPropertiesFormat))
Swagger.property_type(::Type{ ApplicationGatewayRedirectConfigurationPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayRedirectConfigurationPropertiesFormat[name]))}
Swagger.field_name(::Type{ ApplicationGatewayRedirectConfigurationPropertiesFormat }, property_name::Symbol) =  _property_map_ApplicationGatewayRedirectConfigurationPropertiesFormat[property_name]

function check_required(o::ApplicationGatewayRedirectConfigurationPropertiesFormat)
    true
end

function validate_property(::Type{ ApplicationGatewayRedirectConfigurationPropertiesFormat }, name::Symbol, val)
end


text:
@testset "917.reverse-only-letters.jl" begin
    @test reverse_only_letters("Test1ng-Leet=code-Q!") == "Qedo1ct-eeLg=ntse-T!"
    @test reverse_only_letters("a-bC-dEf-ghIj") == "j-Ih-gfE-dCba"
    @test reverse_only_letters("ab-cd") == "dc-ba"
end


text:
@inline vreduce(::typeof(+), v::VectorizationBase.AbstractSIMDVector) = vsum(v)
@inline vreduce(::typeof(*), v::VectorizationBase.AbstractSIMDVector) = vprod(v)
@inline vreduce(::typeof(max), v::VectorizationBase.AbstractSIMDVector) = vmaximum(v)
@inline vreduce(::typeof(min), v::VectorizationBase.AbstractSIMDVector) = vminimum(v)
@inline vreduce(op, v::VectorizationBase.AbstractSIMDVector) = vec_vreduce(op, v)
@inline vec_reduce(op, v::VectorizationBase.AbstractSIMDVector) = vec_reduce(op, Vec(v))
vec_vreduce(op, v::Vec{1}) = VectorizationBase.extractelement(v, 0)
@inline function vec_vreduce(op, v::Vec{W}) where {W}
    a = op(VectorizationBase.extractelement(v,0), VectorizationBase.extractelement(v, 1))
    for i ∈ 2:W-1
        a = op(a, VectorizationBase.extractelement(v, i))
    end
    a
end

function mapreduce_simple(f::F, op::OP, args::Vararg{AbstractArray,A}) where {F,OP,A}
    ptrargs = ntuple(a -> pointer(args[a]), Val(A))
    N = length(first(args))
    iszero(N) && throw("Length of vector is 0!")
    st = ntuple(a -> VectorizationBase.static_sizeof(eltype(args[a])), Val(A))
    a_0 = f(VectorizationBase.__vload.(ptrargs, False(), register_size())...); i = 1
    while i < N
        a_0 = op(a_0, f(VectorizationBase.__vload.(ptrargs, VectorizationBase.lazymul.(st, i), False(), register_size())...)); i += 1
    end
    a_0
end

"""
    vmapreduce(f, op, A::DenseArray...)

Vectorized version of `mapreduce`. Applies `f` to each element of the arrays `A`, and reduces the result with `op`.
"""
@inline function vmapreduce(f::F, op::OP, arg1::AbstractArray{T}, args::Vararg{AbstractArray{T},A}) where {F,OP,T<:NativeTypes,A}
    if !(check_args(arg1, args...) && all_dense(arg1, args...))
        return mapreduce(f, op, arg1, args...)
    end
    N = length(arg1)
    iszero(A) || @assert all(length.(args) .== N)
    W = VectorizationBase.pick_vector_width(T)
    if N < W
        mapreduce_simple(f, op, arg1, args...)
    else
        _vmapreduce(f, op, W, N, T, arg1, args...)
    end
end
@inline function _vmapreduce(f::F, op::OP, ::StaticInt{W}, N, ::Type{T}, args::Vararg{AbstractArray{<:NativeTypes},A}) where {F,OP,A,W,T}
    ptrargs = VectorizationBase.zero_offsets.(stridedpointer.(args))
    if N ≥ 4W
        index = VectorizationBase.Unroll{1,W,4,1,W,0x0000000000000000}((Zero(),)); i = 4W
        au = f(vload.(ptrargs, index)...)
        while i < N - ((W << 2) - 1)
            index = VectorizationBase.Unroll{1,W,4,1,W,0x0000000000000000}((i,)); i += 4W
            au = op(au, f(vload.(ptrargs, index)...))
        end
        a_0 = VectorizationBase.reduce_to_onevec(op, au)
    else
        a_0 = f(vload.(ptrargs, ((MM{W}(Zero()),),))...); i = W
    end
    while i < N - (W - 1)
        a_0 = op(a_0, f(vload.(ptrargs, ((MM{W}(i),),))...)); i += W
    end
    if i < N
        m = mask(T, N & (W - 1))
        a_0 = ifelse(m, op(a_0, f(vload.(ptrargs, ((MM{W}(i),),))...)), a_0)
    end
    vreduce(op, a_0)
end
@inline vmapreduce(f, op, args...) = mapreduce(f, op, args...)

length_one_axis(::Base.OneTo) = Base.OneTo(1)
length_one_axis(::Any) = 1:1

"""
    vreduce(op, destination, A::DenseArray...)

Vectorized version of `reduce`. Reduces the array `A` using the operator `op`.
"""
@inline vreduce(op, arg) = vmapreduce(identity, op, arg)

for (op, init) in zip((:+, :max, :min), (:zero, :typemin, :typemax))
    @eval @inline function vreduce(::typeof($op), arg; dims = nothing)
        if !(check_args(arg) && all_dense(arg))
            return reduce($op, arg, dims = dims)
        end
        dims === nothing && return _vreduce($op, arg)
        isone(ndims(arg)) && return [_vreduce($op, arg)]
        @assert length(dims) == 1
        axes_arg = axes(arg)
        axes_out = Base.setindex(axes_arg, length_one_axis(axes_arg[dims]), dims)
        out = similar(arg, axes_out)
        # fill!(out, $init(first(arg)))
        # TODO: generated function with Base.Cartesian.@nif to set to ndim(arg)
        Base.Cartesian.@nif 5 d -> (d <= ndims(arg) && dims == d) d -> begin
            Rpre = CartesianIndices(ntuple(i -> axes_arg[i], d-1))
            Rpost = CartesianIndices(ntuple(i -> axes_arg[i+d], ndims(arg) - d))
            _vreduce_dims!(out, $op, Rpre, 1:size(arg, dims), Rpost, arg)
        end d -> begin
            Rpre = CartesianIndices(axes_arg[1:dims-1])
            Rpost = CartesianIndices(axes_arg[dims+1:end])
            _vreduce_dims!(out, $op, Rpre, 1:size(arg, dims), Rpost, arg)
        end
    end

    @eval function _vreduce_dims!(out, ::typeof($op), Rpre, is, Rpost, arg)
        s = $init(first(arg))
        @turbo for Ipost in Rpost, Ipre in Rpre
            accum = s
            for i in is
                accum = $op(accum, arg[Ipre, i, Ipost])
            end
            out[Ipre, 1, Ipost] = accum
        end
        return out
    end

    @eval function _vreduce(::typeof($op), arg)
        s = $init(first(arg))
        @turbo for i in eachindex(arg)
            s = $op(s, arg[i])
        end
        return s
    end
end


text:
module Parsers # Bukdu.Plug

import ...Deps: HTTP, URIParser, Request
import .HTTP.Messages: hasheader, header

# getindex_header

function getindex_header(headers::Vector{Pair{String,String}}, key::String)::Union{String,Nothing}
    if hasheader(headers, key)
        header(headers, key)
    else
        nothing
    end
end

# application/x-www-form-urlencoded

struct UrlEncodedScanner
    data::Vector{UInt8}
end

# https://github.com/JuliaWeb/HttpCommon.jl/blob/v0.2.6/src/HttpCommon.jl#L141
function scan(s::UrlEncodedScanner)::Vector{Pair{String,String}}
    query = String(s.data)
    ps = Vector{Pair{String,String}}()
    isempty(query) && return assoc
    for field in split(query, '&'; keepempty=false)
        (k, v) = split(field, '=')
        key = URIParser.unescape_form(k)
        value = URIParser.unescape_form(v)
        push!(ps, Pair(key, value))
    end
    ps
end

# multipart/form-data

mutable struct FormScanner
    data::Vector{UInt8}
    boundary::String
    FormScanner(data::Vector{UInt8}, boundary::String) = new(data, boundary)
end

const CR = 0x0d
const LF = 0x0a
const ContentDisposition = "Content-Disposition: form-data;"

function rstripcr(vec::Vector{UInt8})::Vector{UInt8}
    if !isempty(vec) && vec[end] == CR
        vec[1:end-1]
    else
        vec
    end
end

function rstripcrlf(vec::Vector{UInt8})::Vector{UInt8}
    if !isempty(vec) && vec[end] == LF
        rstripcr(vec[1:end-1])
    else
        vec
    end
end

function content_disposition_fields(io::IOBuffer)::Dict{String,String}
    fields = Dict{String,String}()
    str = String(readuntil(io, LF))
    h = rstripcr(readuntil(io, LF))
    if isempty(h)
    else
        (k, v) = split(String(h), ": ")
        fields[k] = v
        readuntil(io, LF)
    end
    pat = r" (?P<key>[^\"]*)=\"(?P<val>[^\"]*)\""
    for m in eachmatch(pat, str)
        fields[m[:key]] = m[:val]
    end
    fields
end

function scan(s::FormScanner)::Vector{Pair{String,String}}
    ps = Vector{Pair{String, String}}()
    io = IOBuffer(s.data, read=true, truncate=false)
    while !eof(io)
        chunk = readuntil(io, s.boundary)
        if !isempty(chunk)
            chunkio = IOBuffer(chunk)
            skip(chunkio, CR)
            skip(chunkio, LF)
            contentdisposition = readuntil(chunkio, ContentDisposition)
            if !isempty(contentdisposition)
            contentio = IOBuffer(contentdisposition)
            fields = content_disposition_fields(contentio)
            value = String(rstripcrlf(read(contentio)))
            push!(ps, Pair(fields["name"], value))
            end
        end
    end
    ps
end

function fetch_body_params(req::Request)::Vector{Pair{String,String}}
    if hasheader(req.headers, "Content-Type")
        content_type = header(req.headers, "Content-Type")
        if "application/x-www-form-urlencoded" == content_type
            scanner = UrlEncodedScanner(req.body)
            return scan(scanner)
        elseif startswith(content_type, "multipart/form-data")
            pat = r"boundary=\"?([\W\w]*)\"?"
            m = match(pat, content_type)
            if m isa RegexMatch
                boundary = m[1]
                scanner = FormScanner(req.body, string("--", boundary))
                return scan(scanner)
            end
        end
    end
    Vector{Pair{String,String}}()
end

end # module Bukdu.Plug.Parsers

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
mutable struct GFRegenerateState
    prev_trace::GFTrace
    trace::GFTrace
    selection::AddressSet
    score::Float64
    weight::Float64
    visitor::AddressVisitor
    params::Dict{Symbol,Any}
    args_change::Any
    retchange::Union{Some{Any},Nothing}
    callee_output_changes::HomogenousTrie{Any,Any}
end

function GFRegenerateState(args_change, prev_trace, selection, params)
    visitor = AddressVisitor()
    GFRegenerateState(prev_trace, GFTrace(), selection, 0., 0., visitor,
                    params, args_change, nothing, HomogenousTrie{Any,Any}())
end

get_args_change(state::GFRegenerateState) = state.args_change

function set_ret_change!(state::GFRegenerateState, value)
    if state.retchange === nothing
        state.retchange = value
    else
        error("@retchange! was already used")
    end
end

function get_addr_change(state::GFRegenerateState, addr)
    get_leaf_node(state.callee_output_changes, addr)
end

function addr(state::GFRegenerateState, dist::Distribution{T}, args, addr, args_change) where {T}
    visit!(state.visitor, addr)
    has_previous = has_primitive_call(state.prev_trace, addr)
    if has_internal_node(state.selection, addr)
        error("Got internal node but expected leaf node in selection at $addr")
    end
    in_selection = has_leaf_node(state.selection, addr)
    local retval::T
    if has_previous && !in_selection
        # there was a previous value, and it was not in the selection
        # use the previous value
        # it contributes to the weight
        prev_call = get_primitive_call(state.prev_trace, addr)
        retval = prev_call.retval
        score = logpdf(dist, retval, args...)
        state.weight += score - prev_call.score
        retchange = NoChange()
    elseif has_previous
        # there was a previous value, and it was in the selection
        # simulate a new value
        # it does not contribute to the weight
        # retchange indicates that there was a change and gives the previous value
        retval = random(dist, args...)
        score = logpdf(dist, retval, args...)
        prev_call = get_primitive_call(state.prev_trace, addr)
        retchange = (true, prev_call.retval)
    else
        # there is no previous value
        # simulate a new value
        # it does not contribute to the weight
        # retchange is nothing, because the address is new
        retval = random(dist, args...)
        score = logpdf(dist, retval, args...)
        retchange = nothing
    end
    set_leaf_node!(state.callee_output_changes, addr, retchange)
    call = CallRecord(score, retval, args)
    state.trace = assoc_primitive_call(state.trace, addr, call)
    state.score += score
    retval
end

function addr(state::GFRegenerateState, gen::Generator{T}, args, addr, args_change) where {T}
    visit!(state.visitor, addr)
    has_previous = has_subtrace(state.prev_trace, addr)
    if has_previous
        if addr in state.selection
            error("Entire sub-traces cannot be selected, tried to select $addr")
        end
        prev_trace = get_subtrace(state.prev_trace, addr)
        prev_call = get_call_record(prev_trace)
        if has_internal_node(state.selection, addr)
            selection = get_internal_node(state.selection, addr)
        else
            selection = EmptyAddressSet()
        end
        (trace, weight, retchange) = regenerate(
            gen, args, args_change, prev_trace, selection)
        state.weight += weight
        set_leaf_node!(state.callee_output_changes, addr, retchange)
    else
        trace = simulate(gen, args)
        set_leaf_node!(state.callee_output_changes, addr, nothing)
    end
    call = get_call_record(trace)
    state.score += call.score
    state.trace = assoc_subtrace(state.trace, addr, trace)
    call.retval::T
end

function splice(state::GFRegenerateState, gf::GenFunction, args::Tuple)
    exec(gf, state, args)
end

function regenerate(gen::GenFunction, new_args, args_change, trace, selection)
    state = GFRegenerateState(args_change, trace, selection, gen.params)
    retval = exec(gen, state, new_args)
    new_call = CallRecord(state.score, retval, new_args)
    state.trace.call = new_call
    (state.trace, state.weight, state.retchange)
end

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
"""
Module Event
"""
module Event

    using Nextion.Return
    using Nextion: PageID, ComponentID
    using Nextion: hasend, ensurehasend

    """
    Module storing Touch Event enum

    Usage:
    - `Event.Touch.Press`
    - `Event.Touch.Release`
    """
    module Touch
        @enum TouchEventCode::UInt8 begin
            Press = 0x01
            Release = 0x00
        end
        function code(val::UInt8)
            TouchEventCode(val)
        end
    end

    abstract type  AbstractNexEvent end

    abstract type ExpectedMessageSize end
    struct SizeUnknown <: ExpectedMessageSize end
    struct HasLength <: ExpectedMessageSize end

    ExpectedMessageSize(::Type{<: AbstractNexEvent}) = SizeUnknown()

    function ensurehasexpectedlength(msg::Vector{UInt8}, typ::Type)
        if ExpectedMessageSize(typ) == HasLength()
            n = length(msg)
            expectedlength = length(typ)
            if n != expectedlength
                error("Event message must have $expectedlength bytes not $n")
            end
        end
    end

    """
        TouchEvent(msg)

    0X65+Page ID+Component ID+TouchEvent+End
    """
    struct TouchEvent <: AbstractNexEvent
        code::Return.Code.ReturnCode
        pid::PageID
        cid::ComponentID
        tevts::Event.Touch.TouchEventCode  # touch event state

        function TouchEvent(msg::Vector{UInt8})
            ensurehasend(msg)
            ensurehasexpectedlength(msg, TouchEvent)
            code = Return.code(msg[1])
            pid = PageID(msg[2])
            cid = ComponentID(msg[3])
            tevts = Event.Touch.code(msg[4])
            new(code, pid, cid, tevts)
        end
    end
    Base.length(::Type{TouchEvent}) = 7

    """
        CurrentPageIDHeadEvent(msg)

    0X66+Page ID+End
    """
    struct CurrentPageIDHeadEvent <: AbstractNexEvent
        code::Return.Code.ReturnCode
        pid::PageID

        function CurrentPageIDHeadEvent(msg::Vector{UInt8})
            ensurehasend(msg)
            ensurehasexpectedlength(msg, CurrentPageIDHeadEvent)
            code = Return.code(msg[1])
            pid = PageID(msg[2])
            new(code, pid)
        end
    end
    Base.length(::Type{CurrentPageIDHeadEvent}) = 5

    function Base.UInt16(x1::UInt8, x2::UInt8)::UInt16
        UInt16(x1) << 8 + UInt16(x2)
    end

    """
        PositionHeadEvent(msg)

    0X67++ Coordinate X High-order+Coordinate X Low-order+Coordinate Y High-order+Coordinate Y Low-order+TouchEvent State+End
    """
    struct PositionHeadEvent <: AbstractNexEvent
        code::Return.Code.ReturnCode
        x::UInt16
        y::UInt16
        tevts::Event.Touch.TouchEventCode

        function PositionHeadEvent(msg::Vector{UInt8})
            ensurehasend(msg)
            ensurehasexpectedlength(msg, PositionHeadEvent)
            code = Return.code(msg[1])
            x = UInt16(msg[2], msg[3])
            y = UInt16(msg[4], msg[5])
            tevts = Event.Touch.code(msg[6])
            new(code, x, y, tevts)
        end
    end
    Base.length(::Type{PositionHeadEvent}) = 9

    """
        SleepPositionHeadEvent(msg)

    0X68++Coordinate X High-order+Coordinate X Low-order+Coordinate Y High-order+Coordinate Y Low-order+TouchEvent State+End
    """
    struct SleepPositionHeadEvent <: AbstractNexEvent
        code::Return.Code.ReturnCode
        x::UInt16
        y::UInt16
        tevts::Event.Touch.TouchEventCode

        function SleepPositionHeadEvent(msg::Vector{UInt8})
            ensurehasend(msg)
            ensurehasexpectedlength(msg, SleepPositionHeadEvent)
            code = Return.code(msg[1])
            x = UInt16(msg[2], msg[3])
            y = UInt16(msg[4], msg[5])
            tevts = Event.Touch.code(msg[6])
            new(code, x, y, tevts)
        end
    end
    Base.length(::Type{SleepPositionHeadEvent}) = 9

    """
        StringHeadEvent(msg)

    0X70+Variable Content in ASCII code+End
    """
    struct StringHeadEvent <: AbstractNexEvent
        code::Return.Code.ReturnCode
        value::String

        function StringHeadEvent(msg::Vector{UInt8})
            ensurehasend(msg)
            ensurehasexpectedlength(msg, StringHeadEvent)
            code = Return.code(msg[1])
            @assert code == first(StringHeadEvent)
            value = String(msg[2:end - 3])
            new(code, value)
        end
    end
    first(::Type{StringHeadEvent}) = Return.Code.STRING_HEAD
    #Base.length(::Type{StringHeadEvent}) = ???  # variable length

    """
        NumberHeadEvent(msg)

    0X71+variable binary data(4 bytes little endian mode, low in front)+End
    """
    struct NumberHeadEvent <: AbstractNexEvent
        code::Return.Code.ReturnCode
        value::UInt32
        signedvalue::Int32

        function NumberHeadEvent(msg::Vector{UInt8})
            ensurehasend(msg)
            ensurehasexpectedlength(msg, NumberHeadEvent)
            code = Return.code(msg[1])
            @assert code == first(NumberHeadEvent)
            value = UInt32(msg[2]) + UInt32(msg[3]) << 8 + UInt32(msg[4]) << 16 + UInt32(msg[5]) << 24
            signedvalue = reinterpret(Int32, value)
            new(code, value, signedvalue)
        end
    end
    first(::Type{NumberHeadEvent}) = Return.Code.NUMBER_HEAD
    Base.length(::Type{NumberHeadEvent}) = 8

end


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

using Markdown
using InteractiveUtils

# ╔═╡ f054ac12-47b3-11eb-0d47-43990d3906a8
using Revise, FourierTools, Plots, FFTW

# ╔═╡ c5555382-4924-11eb-15d1-17cb887ca876
md"#### Preparing data

Please note, that we always exclude `x_max` from the x-positions since the FFT based interpolation is basically a high resolution evaluation of the Fourier series which assumes a periodic signal.
"

# ╔═╡ 0511e9da-47b4-11eb-3bde-c78c192be0cf
begin
	N_low = 64
	x_min = 0.0
	x_max = 8*2π

	xs_low = range(x_min, x_max, length=N_low+1)[1:N_low]
	xs_high = range(x_min, x_max, length=5000)[1:end-1]
	f(x) = sin(0.5*x) + cos(x) + cos(2 * x) + sin(0.25*x)
	arr_low = f.(xs_low)
	arr_high = f.(xs_high)
end

# ╔═╡ 5e5cc26a-49eb-11eb-01a7-150d675fee58
md"#### Calculate the interpolation"

# ╔═╡ b8fe21d8-47b5-11eb-1854-09254d3cbba1
begin
	N = 1000
	xs_interp = range(x_min, x_max, length=N+1)[1:N]
	arr_interp = resample(arr_low, N)
end

# ╔═╡ bf11f160-4924-11eb-3812-777bb547f499
md"#### Plot the final results"

# ╔═╡ 0511d792-47b4-11eb-316a-3d0148d68406
begin
	scatter(xs_low, arr_low, legend=:bottomleft, markersize=2, label="Low sampling")
	plot!(xs_interp, arr_interp, label="FFT based sinc interpolation", linestyle=:dash)
	plot!(xs_high, arr_high, linestyle=:dashdotdot, label="High sampling")
end

# ╔═╡ 20f9491e-511b-11eb-1a6e-c3af6f5e11a5
md"#### Downsampling
32 samples in the downsampled signal should be sufficient for Nyquist sampling.
And as we can see, the downsampled signal still matches the original one.
"

# ╔═╡ 2e7e1800-511b-11eb-3334-4ddf5076143e
begin
	N_ds = 32
	xs_ds = range(x_min, x_max, length=N_ds+1)[1:N_ds]
	arr_ds = real(resample(arr_high, (N_ds,)))
end

# ╔═╡ 29db9534-511b-11eb-1bdf-37ea0594681f
begin
	scatter(xs_low, arr_low, legend=:bottomleft, markersize=2, label="Low sampling")
	plot!(xs_interp, arr_interp, label="FFT based sinc interpolation", linestyle=:dash)
	plot!(xs_ds, arr_ds, label="downsampled array", linestyle=:dot)
end

# ╔═╡ Cell order:
# ╠═f054ac12-47b3-11eb-0d47-43990d3906a8
# ╟─c5555382-4924-11eb-15d1-17cb887ca876
# ╠═0511e9da-47b4-11eb-3bde-c78c192be0cf
# ╟─5e5cc26a-49eb-11eb-01a7-150d675fee58
# ╠═b8fe21d8-47b5-11eb-1854-09254d3cbba1
# ╟─bf11f160-4924-11eb-3812-777bb547f499
# ╠═0511d792-47b4-11eb-316a-3d0148d68406
# ╠═20f9491e-511b-11eb-1a6e-c3af6f5e11a5
# ╠═2e7e1800-511b-11eb-3334-4ddf5076143e
# ╠═29db9534-511b-11eb-1bdf-37ea0594681f


text:
export NFFTOp
import Base.adjoint

mutable struct NFFTOp{T} <: AbstractLinearOperator{T}
  nrow :: Int
  ncol :: Int
  symmetric :: Bool
  hermitian :: Bool
  prod :: Function
  tprod :: Nothing
  ctprod :: Function
  nprod :: Int
  ntprod :: Int
  nctprod :: Int
  plan
  toeplitz :: Bool
end

"""
    NFFTOp(shape::Tuple, tr::Trajectory; nodes=nothing, kargs...)

generates a `NFFTOp` which evaluates the MRI Fourier signal encoding operator using the NFFT.

# Arguments:
* `shape::NTuple{D,Int64}`  - size of image to encode/reconstruct
* `tr::Trajectory`          - Trajectory with the kspace nodes to sample
* (`nodes=nothing`)         - Array containg the trajectory nodes (redundant)
* (`kargs`)                   - additional keyword arguments
"""
function NFFTOp(shape::Tuple, tr::Trajectory; nodes=nothing, toeplitz=false,
                oversamplingFactor=1.25, kernelSize=3, kargs...)
  nodes==nothing ? nodes=kspaceNodes(tr) : nothing
  # plan = NFFTPlan(nodes, shape, kernelSize, oversamplingFactor, precompute=NFFT.FULL)
  plan = plan_nfft(nodes, shape, kernelSize, oversamplingFactor, precompute=NFFT.FULL)

  return NFFTOp{ComplexF64}(size(nodes,2), prod(shape), false, false
            , x->produ(plan,x)
            , nothing
            , y->ctprodu(plan,y), 0, 0, 0, plan, toeplitz)
end

function produ(plan::NFFT.NFFTPlan, x::Vector{T}) where T<:Union{Real,Complex}
  y = nfft(plan,reshape(x[:],plan.N))
  return vec(y)
end

function ctprodu(plan::NFFT.NFFTPlan, y::Vector{T}) where T<:Union{Real,Complex}
  x = nfft_adjoint(plan, y[:])
  return vec(x)
end

function Base.copy(S::NFFTOp)
  plan = copy(S.plan)
  return NFFTOp{ComplexF64}(size(plan.x,2), prod(plan.N), false, false
         , x->produ(plan,x)
         , nothing
         , y->ctprodu(plan,y), 0, 0, 0, plan, S.toeplitz)
end

### Normal Matrix Code ###

struct NFFTNormalOp{S,D}
  shape::S
  weights::D
  fftplan
  ifftplan
  λ::Array{ComplexF64}
  xL::Matrix{ComplexF64}
end

function Base.copy(S::NFFTNormalOp)
  fftplan = plan_fft(zeros(ComplexF64, Tuple(2*collect(S.shape)));flags=FFTW.MEASURE)
  ifftplan = plan_ifft(zeros(ComplexF64, Tuple(2*collect(S.shape)));flags=FFTW.MEASURE)
  return NFFTNormalOp(S.shape, S.weights, fftplan, ifftplan, S.λ, S.xL)
end

function NFFTNormalOp(S::NFFTOp, W)
  shape = S.plan.N
  weights = W*ones(S.nrow)

  λ, ft, ift = diagonalizeOp(S.plan, weights)
  xL = zeros(ComplexF64,2*shape[1], 2*shape[2])

  return NFFTNormalOp(shape, W, ft, ift, λ, xL)
end

function SparsityOperators.normalOperator(S::NFFTOp, W=I)
  if S.toeplitz
    return NFFTNormalOp(S,W)
  else
    return NormalOp(S,W)
  end
end

function Base.:*(S::NFFTNormalOp, x::AbstractVector{T}) where T
  shape = S.shape

  # xL = zeros(T,Tuple(2*collect(shape)))
  # xL[1:shape[1],1:shape[2]] = x
  # λ = reshape(S.λ,Tuple(2*collect(shape)))

  # y = (S.ifftplan*( λ.*(S.fftplan*xL)))[1:shape[1],1:shape[2]]

  S.xL .= 0
  S.xL[1:shape[1],1:shape[2]] .= reshape(x,shape)
  λ = reshape(S.λ,Tuple(2*collect(shape)))

  y = (S.ifftplan*( λ.*(S.fftplan*S.xL)))[1:shape[1],1:shape[2]]

  return vec(y)
end

function diagonalizeOp(p::NFFT.NFFTPlan, weights=nothing)
  shape = p.N
  nodes = p.x

  # plan the FFTs
  fftplan = plan_fft(zeros(ComplexF64, Tuple(2*collect(shape)));flags=FFTW.MEASURE)
  ifftplan = plan_ifft(zeros(ComplexF64, Tuple(2*collect(shape)));flags=FFTW.MEASURE)

  # calculate first column of block Toeplitz matrix
   e1 = zeros(ComplexF64,shape)
   e1[1] = 1.
   firstCol = nfft_adjoint(p, nfft(p,e1)[:].*weights)
   firstCol = reshape(firstCol, shape)

  # calculate first rows of the leftmost Toeplitz blocks
  p2 = plan_nfft([-1,1].*nodes, shape, 3, 1.25)
  firstRow = nfft_adjoint(p2, nfft(p2,e1)[:].*weights)
  firstRow = reshape(firstRow, shape)

  # construct FT matrix of the eigentvalues
  # here a row and column in the center is filled with zeros.
  # this is done to ensure that the FFT can operate
  # on arrays with even size in each dimension
  eigMat = zeros(ComplexF64, Tuple(2*collect(shape)))
  eigMat[1:shape[1], 1:shape[2]] .= firstCol
  eigMat[shape[1]+2:end, 1:shape[2]] .= firstRow[end:-1:2,:]
  eigMat[1:shape[1], shape[2]+2:end] .= conj.(firstRow[:,end:-1:2])
  eigMat[shape[1]+2:end, shape[2]+2:end] .=  conj.(firstCol[end:-1:2,end:-1:2])

  return vec(fftplan*eigMat), fftplan, ifftplan
end

#
# calculate the matrix element A_{j,k} explicitely
#
#function getMatrixElement(j::Int, k::Int, shape::Tuple, nodes::Matrix; weights=nothing)
#  elem=0.
#  x = mod(k-1,shape[1])-mod(j-1,shape[1])
#  y = div(k-1,shape[1])-div(j-1,shape[1])
#  if weights != nothing
#    for i=1:size(nodes,2)
#      elem += exp( -2*pi*1im*(nodes[1,i]*x + nodes[2,i]*y) )*weights[i]
#    end
#  else
#    for i=1:size(nodes,2)
#      elem += exp( -2*pi*1im*(nodes[1,i]*x + nodes[2,i]*y) )
#    end
#  end
#
#  return elem
#end


text:
# Autogenerated wrapper script for GAP_pkg_json_jll for aarch64-apple-darwin
export json

using GAP_jll
using GAP_lib_jll
JLLWrappers.@generate_wrapper_header("GAP_pkg_json")
JLLWrappers.@declare_file_product(json)
function __init__()
    JLLWrappers.@generate_init_header(GAP_jll, GAP_lib_jll)
    JLLWrappers.@init_file_product(
        json,
        "lib/gap/json.so",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module UdKaDyn

include("RigidBody.jl")
include("Joint.jl")
include("OrientationConversion.jl")
include("extF.jl")
include("simulate.jl")
include("Force.jl")

include("linearize.jl")
include("nRotor.jl")
include("plotSol.jl")
# include("PointMass.jl")
include("trim_FinDiff.jl")
include("trim_kron.jl")
include("trim_kronLazy.jl")

export gen_nRotor, simulate, RigidBody, InertialFrameAsRB, Joint, checkRevJointIn, plotErrNorm, plotPos, plotQuat, plotVel, plotAngVel, linearize, initialiseRigidBody!
end # module

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
using Parameters
using StaticArrays
using LinearAlgebra
using Rotations
using BenchmarkTools
using MeshCat
using CoordinateTransformations
using GeometryTypes
using TrajOptPlots
const TO = TrajectoryOptimization
import TrajectoryOptimization: Rotation

# Set up visualizer
vis = Visualizer(); open(vis);
dims = [1,1,3]*0.5
robot = HyperRectangle(Vec((-1 .* dims ./2)...), Vec(dims...));
setobject!(vis["robot"], robot);

# Define Spinner Model
@with_kw struct Spinner{R<:Rotation} <: AbstractModel
    mass::Float64 = 1.0
    J::Diagonal{Float64, SVector{3, Float64}} = Diagonal(@SVector ones(3))
end
Base.size(::Spinner{<:UnitQuaternion}) = (7,3)
Base.size(::Spinner) = (7,3)
Dynamics.orientation(::Spinner{UnitQuaternion{T,D}}, x::SVector{N,T2}) where {T,T2,N,D} =
    UnitQuaternion{T2,D}(x[4], x[5], x[6], x[7])
Dynamics.orientation(::Spinner, x::SVector) = R(x[4], x[5], x[6])

function TrajectoryOptimization.dynamics(model::Spinner, x::SVector, u::SVector)
    J = model.J
    ω = @SVector [x[1], x[2], x[3]]
    q = Dynamics.orientation(model, x)
    ωdot = J\(u - ω × (J*ω))
    qdot = kinematics(q,ω)
    return [ωdot; qdot]
end

function TrajectoryOptimization.state_diff_jacobian!(G, model::Spinner, Z::Traj)
    for k in eachindex(Z)
        G[k] = TrajectoryOptimization.state_diff_jacobian(model, state(Z[k]))
    end
end

function TrajectoryOptimization.state_diff_jacobian(::Spinner{<:UnitQuaternion}, x::SVector)
    q = UnitQuaternion(x[4], x[5], x[6], x[7])
    G = TrajectoryOptimization.∇differential(q)
    @SMatrix [1 0 0 0 0 0;
    0 1 0 0 0 0;
    0 0 1 0 0 0;
    0 0 0 G[1] G[5] G[ 9];
    0 0 0 G[2] G[6] G[10];
    0 0 0 G[3] G[7] G[11];
    0 0 0 G[4] G[8] G[12]]
end

function TrajectoryOptimization.state_diff_jacobian(
        ::Spinner{UnitQuaternion{T,IdentityMap}}, x::SVector) where T
    return I
end

function TrajectoryOptimization.state_diff(model::Spinner, x, x0)
    δω = @SVector [x[1]-x0[1], x[2]-x0[2], x[3]-x0[3]]
    q  = Dynamics.orientation(model, x)
    q0 = Dynamics.orientation(model, x0)
    inv(q0)*q
    δq = VectorPart(inv(q0)*q)
    return @SVector [δω[1], δω[2], δω[3], δq[1], δq[2], δq[3]]
end

TrajectoryOptimization.state_diff(::Spinner{UnitQuaternion{T,IdentityMap}}, x, x0) where T = x - x0

TrajectoryOptimization.state_diff_size(model::Spinner) = 6
TrajectoryOptimization.state_diff_size(model::Spinner{UnitQuaternion{T,IdentityMap}}) where T = 7

function TrajOptPlots.visualize!(vis, model::Spinner, Z::Traj)
    X = states(Z)
    fps = Int(round(length(Z)/Z[end].t))
    anim = MeshCat.Animation(fps)
    for k in eachindex(Z)
        atframe(anim, k) do
            x = @SVector zeros(3)
            q = UnitQuaternion(Dynamics.orientation(model, X[k]))
            settransform!(vis["robot"], compose(Translation(x), LinearMap(Quat(q))))
        end
    end
    setanimation!(vis, anim)
    return anim
end

model = Spinner{UnitQuaternion{Float64,MRPMap}}()
Dynamics.rotation_type(::Spinner{R}) where R = R

# Params
N = 101
tf = 5.

# Initial and Final conditions
x0 = @SVector [0,0,0, 1,0,0,0.]
xf = @SVector [0,0,0, 1,0,0,0.]

# Intermediate frames
qs = [UnitQuaternion(Quat(RotZ(ang))) for ang in range(0,2pi,length=N)]
X_guess = map(1:N) do k
    om = @SVector [x0[1], x0[2], x0[3]]
    if k > 76
        q = -qs[k]
    else
        q = qs[k]
    end
    [om; SVector(q)]
end

# Build Objective
Q_diag = [(@SVector fill(1e-3,3)); (@SVector fill(1e-1,4))]
R_diag = @SVector fill(1e-4, 3)
no_quat = SVector{6}(deleteat!(collect(1:7), 4))

costs = map(1:N) do k
    x = X_guess[k]
    if Dynamics.rotation_type(model) <: UnitQuaternion{T,IdentityMap} where T
        Q = Diagonal(Q_diag)
    else
        Q = Diagonal(Q_diag[no_quat])
        Gf = TO.state_diff_jacobian(model, x)
        Q = Gf*Q*Gf'
        # Q = Diagonal(Q_diag)
    end
    R = Diagonal(R_diag)

    if k == N
        Q *= 100
    end

    LQRCost(Q,R,x, checks=false)
    # QuadraticCost(Q,R, checks=false)
    # Q = Diagonal(Q_diag[no_quat])
    # RBCost(model, Q, R, x)
    # QuatLQRCost(Q,R,x)
end
obj = Objective(costs)
obj[51].Q

# Problem
model = Spinner{UnitQuaternion{Float64,MRPMap}}()
model = Spinner{UnitQuaternion{Float64,IdentityMap}}()
prob = Problem(model, obj, xf, tf, x0=x0)

# Solve
ilqr = iLQRSolver(prob)
ilqr.opts.verbose = true
U0 = [@SVector zeros(3) for k = 1:N-1]
initial_controls!(ilqr, U0)
solve!(ilqr)
cost(ilqr)
ilqr.K[1]

q = states(ilqr)[N][4:7]
q'xf[4:7]
visualize!(vis, model, ilqr.Z)

# With Quaternion
model_quat = Spinner{UnitQuaternion{VectorPart}}()
Gf = TO.state_diff_jacobian(model_quat, xf)
Q_quat = Gf*Diagonal(diag(Q)[no_quat])*Gf'
no_quat = @SVector [1,2,3,5,6,7]
Qf_quat = Gf*Diagonal(diag(Qf)[no_quat])*Gf'
costfun = QuadraticCost(Q_quat,R, checks=false)
costfun_term = QuadraticCost(Qf_quat, R*0, checks=false)

obj_quat = Objective(costfun, costfun_term, N)

prob_quat = Problem(model_quat, obj_quat, xf, tf, x0=x0)
# Solve
ilqr_quat = iLQRSolver(prob_quat)
ilqr_quat.opts.verbose = true
initial_controls!(ilqr_quat, U0)
solve!(ilqr_quat)
cost(ilqr_quat)
xf[4:7]'states(ilqr_quat)[end][4:7]
visualize!(vis, model, ilqr.Z)

ilqr_quat.opts.verbose = false
@btime begin
    initial_controls!($ilqr_quat, $U0)
    solve!($ilqr_quat)
end

ilqr.opts.verbose = false
@btime begin
    initial_controls!($ilqr, $U0)
    solve!($ilqr)
end

# Step through the solve
solver = iLQRSolver(prob_quat)
initial_controls!(solver, U0)
rollout!(solver)
states(solver)[end]
cost(solver)

Z = get_trajectory(solver)
TO.state_diff_jacobian!(solver.G, solver.model, Z)
TO.discrete_jacobian!(solver.∇F, solver.model, Z)
TO.cost_expansion!(solver.Q, solver.obj, solver.Z)
solver.∇F[50]
solver.Q.x[50]
solver.Q.u[50]
solver.G[50]
ΔV = backwardpass!(solver)
forwardpass!(solver, ΔV, cost(solver))

for k = 1:N
    solver.Z[k].z = solver.Z̄[k].z
end
states(solver)[end]


text:
using AFL
using Test

@testset "AFL" begin
    @testset "run_target" begin
        target = AFL.init_target(joinpath(@__DIR__, "../deps/afl-test"))
        @test AFL.run_target(target, b"1234678") == AFL.Ok
    end

    @testset "run_target_qemu" begin
        target = AFL.init_target(joinpath(@__DIR__, "../deps/afl-shim"), qemu=true)
        @test AFL.run_target(target, b"") == AFL.Ok
    end

    @testset "close" begin
        target = AFL.init_target(joinpath(@__DIR__, "../deps/afl-test"))
        close(target)
    end

    @testset "has_new_bits_pos" begin
        virgin_bits = AFL.InvCoverageMap()
        trace_bits = AFL.CoverageMap()

        trace_bits.edges[3] = 0x1
        @test AFL.has_new_coverage(virgin_bits, trace_bits) == true

        trace_bits.edges[3] = 0x8
        @test AFL.has_new_coverage(virgin_bits, trace_bits) == true
    end

    @testset "has_new_bits_neg" begin
        virgin_bits = AFL.InvCoverageMap()
        trace_bits = AFL.CoverageMap()

        @test AFL.has_new_coverage(virgin_bits, trace_bits) == false

        virgin_bits.edges[2] = 0xFE
        @test AFL.has_new_coverage(virgin_bits, trace_bits) == false

        trace_bits.edges[2] = 0x1
        @test AFL.has_new_coverage(virgin_bits, trace_bits) == false
    end

    @testset "has_new_bits_multiple" begin
        virgin_bits = AFL.InvCoverageMap()
        trace_bits = AFL.CoverageMap()

        virgin_bits.edges[3] = 0b11110101
        trace_bits.edges[3]  = 0b00001101

        @test AFL.has_new_coverage(virgin_bits, trace_bits) == true
    end

    @testset "mark_covered" begin
        virgin_bits = AFL.InvCoverageMap()
        trace_bits = AFL.CoverageMap()

        trace_bits.edges[3] = 0x1

        @test AFL.has_new_coverage(virgin_bits, trace_bits) == true

        AFL.mark_covered!(virgin_bits, trace_bits)

        @test AFL.has_new_coverage(virgin_bits, trace_bits) == false

        backup = copy(virgin_bits)
        AFL.mark_covered!(virgin_bits, trace_bits)

        @test virgin_bits == backup
    end

    @testset "deterministic" begin
        target = AFL.init_target(joinpath(@__DIR__, "../deps/afl-test"))

        AFL.run_target(target, b"12345678")
        cov1 = AFL.classify_coverage(target)

        AFL.run_target(target, b"12345678")
        cov2 = AFL.classify_coverage(target)

        @test cov1 == cov2
    end

    @testset "deterministic_qemu" begin
        target = AFL.init_target(joinpath(@__DIR__, "../deps/afl-shim"), qemu=true)

        AFL.run_target(target, b"12345678")
        cov1 = AFL.classify_coverage(target)

        AFL.run_target(target, b"12345678")
        cov2 = AFL.classify_coverage(target)

        @test cov1 == cov2
    end
end


text:
@testset "Test 2d mesh generators" begin
using jFEMTools
import Tensors
jF = jFEMTools
#Test generated mesh
mesh = rectangle_mesh(TriangleCell, (2,2), Tensors.Vec{2}((0.0,0.0)), Tensors.Vec{2}((1.0,1.0)))
@test getncells(mesh) == 8
@test getnvertices(mesh) == 9
for cell_idx in 1:getncells(mesh)
    @test cell_volume(mesh, cell_idx) ≈ 1/8
end
@test get_vertices_matrix(mesh) == [0.0 0.0; 0.5 0.0; 1.0 0.0; 0.0 0.5; 0.5 0.5; 1.0 0.5; 0.0 1.0; 0.5 1.0; 1.0 1.0]
@test get_cell_connectivity_list(mesh) == [(1, 2, 4), (2, 5, 4), (2, 3, 5), (3, 6, 5), (4, 5, 7), (5, 8, 7), (5, 6, 8), (6, 9, 8)]
# Check expected data for cell 1
@test mesh.cells[1].vertices == (1,2,4)
#@test mesh.cells[1].faces == (1,2,3)
#@test [face_orientation(mesh,1,i) for i in 1:3] == [true,false,true]
@test cell_diameter(mesh,1) == sqrt(2)/2
@test getnedges(mesh.cells[1]) == 3

# Quad mesh
mesh = rectangle_mesh(RectangleCell, (2,2), Tensors.Vec{2}((0.0,0.0)), Tensors.Vec{2}((1.0,1.0)))

# Hex mesh
mesh = unitSquareMesh(HexagonCell,(2,2))
@test getncells(mesh) == 10
@test getnvertices(mesh) == 19
# Check expected data for cell 1
@test cell_diameter(mesh,1) == sqrt(2)/4
@test getnedges(mesh.cells[1]) == 4

# Mixed mesh
cells = [
	TriangleCell((1,2,4)),
	RectangleCell((2,3,6,5)),
	TriangleCell((5,4,2)),
    RectangleCell((4, 5, 8, 7)),
	RectangleCell((5,6,9,8))
    ];
mixmesh = PolytopalMesh(cells, mesh.vertices)

# Mesh2.0
mesh = rectangle_mesh2(TriangleCell, (2,2), Tensors.Vec{2}((0.0,0.0)), Tensors.Vec{2}((1.0,1.0)))
@test getncells(mesh) == 8
@test getnvertices(mesh) == 9
@test cell_diameter(mesh,1) == sqrt(2)/2
@test getnfacets(mesh,1) == 3
for cell_idx in 1:getncells(mesh)
    @test cell_volume(mesh, cell_idx) ≈ 1/8
end
mesh = rectangle_mesh2(RectangleCell, (2,2), Tensors.Vec{2}((0.0,0.0)), Tensors.Vec{2}((1.0,1.0)))
@test getncells(mesh) == 4
@test getnvertices(mesh) == 9
@test cell_diameter(mesh,1) == sqrt(2)/2
@test getnfacets(mesh,1) == 4

mesh = rectangle_mesh2(HexagonCell, (2,2), Tensors.Vec{2}((0.0,0.0)), Tensors.Vec{2}((1.0,1.0)))
@test getncells(mesh) == 10
@test getnvertices(mesh) == 19
@test cell_diameter(mesh,1) == sqrt(2)/4
@test getnfacets(mesh,1) == 4

end

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
using Documenter
using EffectSizes
using Test

DocMeta.setdocmeta!(
    EffectSizes,
    :DocTestSetup,
    :(using EffectSizes);
    recursive=true
)

# Only test one Julia version to avoid differences due to changes in printing.
if v"1.6" ≤ VERSION
    EffectSizes._update_module_doc()
    doctest(EffectSizes)
else
    @warn "Skipping doctests"
end

@testset "EffectSizes.jl" begin
    include("test_confint.jl")
    include("test_effectsize.jl")
end


text:
println("Begin importing packages for `runtests.jl`...")
import AircraftSystems
AS = AircraftSystems
CC = AS.CC
FM = AS.FLOWMath
LS = AS.LaTeXStrings
LA = AS.LinearAlgebra
plt = AS.PyPlot
# plt.pygui(false)
T = AS.Test
VL = AS.VortexLattice
println("Finished.")

plot_directory = joinpath(AS.topdirectory, "data", "plots", AS.TODAY)
if !isdir(plot_directory); mkpath(plot_directory); end

const save_figs = false #! Set true if you need to save the test figures.

# prepare Epema rotor
include("EpemaData.jl")
contour_directory = joinpath(AS.topdirectory, "test", "data", "airfoil", "contours")
polar_directory = joinpath(AS.topdirectory, "test", "data", "airfoil", "polars")
rs_desired = [0.207, 0.3, 0.4, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0]
# EpemaData.scans_to_interpolated_contours(contour_directory, plot_directory, rs_desired)

# Fetch PROWIM data.
include("PROWIMData.jl")

# T.@testset "rotors" begin

include("epema_rotor_sweep.jl")

include("rotor_wakefunction.jl") #currently doesn't work for me, I think the parameters need adjusting first but I'll let you handle it

# end # test rotors

# T.@testset "wings" begin

include("prowim_propsoff_clalpha.jl")

include("prowim_propsoff_liftdist.jl")

# include("prowim_test_cl_CL.jl")

# end

# @testset "Epema Wing Validation" begin

# include("epema_liftdist.jl")

# end # Epema Blown Wing Validation

# T.@testset "Epema Blown Wing Validation" begin #* still a work in progress. No errors that keep it from running, just not great results yet.

    # include("epema_blown_wing_cl.jl")
    include("epema_blown_wing_cl_epema_polar.jl") # uses the Epema polars from above

# end

# T.@testset "Veldhuis Blown Wing Validation" begin

include("prowim_propson_liftdist.jl")

# end

# benchmarks

include("benchmark_ccblade.jl")

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


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

# TODO: remove this type stub between 0.7 and 1.0

"""
    ConjArray(array)

A lazy-view wrapper of an `AbstractArray`, taking the elementwise complex conjugate. This
type is usually constructed (and unwrapped) via the [`conj`](@ref) function (or related
[`adjoint`](@ref)), but currently this is the default behavior for `RowVector` only. For
other arrays, the `ConjArray` constructor can be used directly.

# Examples
```jldoctest
julia> ConjArray([1+im 0; 0 1-im])
2×2 ConjArray{Complex{Int64},2,Array{Complex{Int64},2}}:
 1-1im  0+0im
 0+0im  1+1im
```
"""
struct ConjArray{T,N,A<:AbstractArray} <: AbstractArray{T,N}
    parent::A
end
const ConjVector{T,V<:AbstractVector} = ConjArray{T,1,V}
const ConjMatrix{T,M<:AbstractMatrix} = ConjArray{T,2,M}

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
###############################################################
# SymTensorValue Type
###############################################################

"""
Type representing a symmetric fourth-order tensor
"""
struct SymFourthOrderTensorValue{D,T,L} <: MultiValue{Tuple{D,D,D,D},T,4,L}
  data::NTuple{L,T}
  function SymFourthOrderTensorValue{D,T}(data::NTuple{L,T}) where {D,T,L}
    @assert L == (D*(D+1)÷2)^2
    new{D,T,L}(data)
  end
end

###############################################################
# Constructors (SymTensorValue)
###############################################################

# Empty SymTensorValue constructor

SymFourthOrderTensorValue()                   = SymFourthOrderTensorValue{0,Int}(NTuple{0,Int}())
SymFourthOrderTensorValue{0}() where {T}      = SymFourthOrderTensorValue{0,Int}(NTuple{0,Int}())
SymFourthOrderTensorValue{0,T}() where {T}    = SymFourthOrderTensorValue{0,T}(NTuple{0,T}())
SymFourthOrderTensorValue(data::NTuple{0})    = SymFourthOrderTensorValue{0,Int}(data)
SymFourthOrderTensorValue{0}(data::NTuple{0}) = SymFourthOrderTensorValue{0,Int}(data)

# SymTensorValue single NTuple argument constructor

@generated function SymFourthOrderTensorValue(data::NTuple{L,T}) where {L,T}
  msg = "Invalid number of scalar arguments in SymFourthOrderTensorValue constructor"
  V = (sqrt(1+8*sqrt(L))-1)/2
  @assert floor(Int,V) == ceil(Int,V) msg
  D = Int(V)
  quote
    SymFourthOrderTensorValue{$D,T}(data)
  end
end
SymFourthOrderTensorValue{D}(data::NTuple{L,T}) where {D,L,T}           = SymFourthOrderTensorValue{D,T}(data)
SymFourthOrderTensorValue{D,T1}(data::NTuple{L,T2}) where {D,L,T1,T2}   = SymFourthOrderTensorValue{D,T1}(NTuple{L,T1}(data))
SymFourthOrderTensorValue{D,T1,L}(data::NTuple{L,T2}) where {D,L,T1,T2} = SymFourthOrderTensorValue{D,T1}(NTuple{L,T1}(data))

# SymTensorValue single Tuple argument constructor

SymFourthOrderTensorValue(data::Tuple) = SymFourthOrderTensorValue(promote(data...))
SymFourthOrderTensorValue{D}(data::Tuple) where {D} = SymFourthOrderTensorValue{D}(promote(data...))
SymFourthOrderTensorValue{D,T1}(data::Tuple) where {D,T1} = SymFourthOrderTensorValue{D,T1}(NTuple{length(data),T1}(data))

# SymTensorValue Vararg constructor

SymFourthOrderTensorValue(data...) = SymFourthOrderTensorValue(data)
SymFourthOrderTensorValue{D}(data...) where {D} = SymFourthOrderTensorValue{D}(data)
SymFourthOrderTensorValue{D,T1}(data...) where {D,T1} = SymFourthOrderTensorValue{D,T1}(data)

###############################################################
# Conversions (SymTensorValue)
###############################################################

# Direct conversion
convert(::Type{<:SymFourthOrderTensorValue{D,T}}, arg::Tuple) where {D,T} = SymFourthOrderTensorValue{D,T}(arg)

# Inverse conversion
convert(::Type{<:NTuple{L,T}}, arg::SymFourthOrderTensorValue) where {L,T} = NTuple{L,T}(Tuple(arg))

# Internal conversion
convert(::Type{<:SymFourthOrderTensorValue{D,T}}, arg::SymFourthOrderTensorValue{D}) where {D,T} = SymFourthOrderTensorValue{D,T}(Tuple(arg))
convert(::Type{<:SymFourthOrderTensorValue{D,T}}, arg::SymFourthOrderTensorValue{D,T}) where {D,T} = arg

###############################################################
# Other constructors and conversions (SymTensorValue)
###############################################################

@generated function zero(::Type{<:SymFourthOrderTensorValue{D,T}}) where {D,T}
  L=(D*(D+1)÷2)^2
  quote
    SymFourthOrderTensorValue{D,T}(tfill(zero(T),Val{$L}()))
  end
end
zero(::Type{<:SymFourthOrderTensorValue{D,T,L}}) where {D,T,L} = SymFourthOrderTensorValue{D,T}(tfill(zero(T),Val{L}()))
zero(::SymFourthOrderTensorValue{D,T,L}) where {D,T,L} = zero(SymFourthOrderTensorValue{D,T,L})

# This is in fact the "symmetrized" 4th order identity
@generated function one(::Type{<:SymFourthOrderTensorValue{D,T}}) where {D,T}
  S = typeof(one(T)/2)
  str = join(["($i==$k && $j==$l) ?  ( $i==$j ? one($S) :  one(T)/2) : zero($S), " for i in 1:D for j in i:D for k in 1:D for l in k:D])
  Meta.parse("SymFourthOrderTensorValue{D,$S}(($str))")
end
one(::SymFourthOrderTensorValue{D,T}) where {D,T} = one(SymFourthOrderTensorValue{D,T})

change_eltype(::Type{SymFourthOrderTensorValue{D,T1,L}},::Type{T2}) where {D,T1,T2,L} = SymFourthOrderTensorValue{D,T2,L}
change_eltype(::SymFourthOrderTensorValue{D,T1,L},::Type{T2}) where {D,T1,T2,L} = change_eltype(SymFourthOrderTensorValue{D,T1,L},T2)

###############################################################
# Introspection (SymTensorValue)
###############################################################

eltype(::Type{<:SymFourthOrderTensorValue{D,T}}) where {D,T} = T
eltype(::SymFourthOrderTensorValue{D,T}) where {D,T} = eltype(SymFourthOrderTensorValue{D,T})

size(::Type{<:SymFourthOrderTensorValue{D}}) where {D} = (D,D,D,D)
size(::SymFourthOrderTensorValue{D}) where {D} = size(SymFourthOrderTensorValue{D})

length(::Type{<:SymFourthOrderTensorValue{D}}) where {D} = D*D*D*D
length(::SymFourthOrderTensorValue{D}) where {D} = length(SymFourthOrderTensorValue{D})

num_components(::Type{<:SymFourthOrderTensorValue{D}}) where {D} = length(SymFourthOrderTensorValue{D})
num_components(::SymFourthOrderTensorValue{D}) where {D} = num_components(SymFourthOrderTensorValue{D})