text:
using CompileBot

snoop_bot(
    BotConfig(
        "Plots",
        yml_path= "SnoopCompile.yml",
        else_os = "linux",
        else_version = "1.4",
    ),
    joinpath(@__DIR__, "precompile_script.jl"),
)


text:
include("common.jl")

Δt = 0.01
mmesh = Mesh(30, 30, 1/30)
bdnode = bcnode(mmesh)
bdnode = [bdnode; bdnode .+ mmesh.ndof]
nu = constant(0.0001ones(get_ngauss(mmesh)))
nodes = fem_nodes(mmesh)
u = @. sin(2π * nodes[:,1])
v = @. sin(2π * nodes[:,2])
u0 = [u;v]
u0[bdnode] .= 0.0
us = solve_burgers(u0, 10, nu)

sess = Session(); init(sess)
U = run(sess, us)

matwrite("fenics/fwd1.mat", Dict("U"=>U))
close("all")
figure(figsize=(10,3))
subplot(121)
title("u displacement")
visualize_scalar_on_fem_points(U[end,1:mmesh.nnode], mmesh)
subplot(122)
title("v displacement")
visualize_scalar_on_fem_points(U[end,mmesh.ndof + 1:mmesh.ndof + mmesh.nnode], mmesh)
savefig("fenics/fwd1.png")

close("all")
visualize_vector_on_fem_points(U[end,1:mmesh.nnode], U[end,mmesh.ndof + 1:mmesh.ndof + mmesh.nnode], mmesh)
savefig("fenics/fwd1_quiver.png")

ext: jl


text:
export copy, size, sort
export dropmonomorphic, dropmonomorphic!
export dropmultiallelic, dropmultiallelic!
export convertcoord

#=
General use utilities
=#

function Base.copy(data::PopData)
    PopData(deepcopy(data.metadata), copy(data.genodata))
end

function Base.size(data::PopData)
    return (samples = data.metadata.samples, loci = data.metadata.loci)
end

"""
    Base.sort(x::NTuple{N,T}) where N where T <: Signed
Sort the integers within a Tuple and return the sorted Tuple.
"""
function Base.sort(x::NTuple{N,T}) where N where T <: Signed
    Tuple(sort(SVector(x)))
end

"""
    convertcoord(coordinate::String)
Takes non-decimal-degree format as a `String` and returns it as a decimal degree
`Float32`. Can be broadcasted over an array of coordinate strings to convert them.
## Formatting requirements
- Coordinates as a `String` separated by spaces (`"11 43 41"`) or colons (`"11:43:41"`)
- Must use negative sign (`"-11 43.52"`) or single-letter cardinal direction (`"11 43.52W"`)
- Missing data should be coded as `missing` or the string `"missing"`
- Can mix colons and spaces (although it's bad practice)
### Example
```
julia> convertcoord("-41 31.52")
-41.5253f0
julia> convertcoord.(["-41 31.52", "25 11:54S"])
2-element Array{Float32,1}:
-41.5253
-25.1983
```
"""
function convertcoord(coordinate::String)
    lowercase(coordinate) == "missing" && return missing
    coord_strip = replace(uppercase(coordinate), r"[NSEW]" => "")
    split_coord = parse.(Float32, split(coord_strip, r"\s|:"))
    split_coord[2] /= 60.0
    if length(split_coord) == 3
        split_coord[3] /= 3600.0
    end
    conv = mapreduce(abs, +, split_coord)
    # N + E are positive | S + W are negative
    if split_coord[1] < 0 || occursin(r"[SW]", uppercase(coordinate))
        # negative
        return round(conv * -1, digits = 4)
    else
        # positive
        return round(conv, digits = 4)
    end
end

convertcoord(coordinate::Missing) = missing

"""
    dropmonomorphic(data::PopData; silent::Bool = false)
Return a `PopData` object omitting any monomorphic loci. Will inform you which
loci were removed.
"""
function dropmonomorphic(data::PopData; silent::Bool = false)
    all_loci = loci(data)
    prt = Base.Iterators.partition(data.genodata.genotype, length(samplenames(data)))
    monomorphs = string.(all_loci[[length(unique(i))==1 for i in prt]])
    if length(monomorphs) == 0
        #!silent && println("No monomorphic loci found.\n")
        return data
    elseif !silent
        if length(monomorphs) == 1
            println("Removing monomorphic locus: " * monomorphs[1])
            println()
        else
            println("Removing $(length(monomorphs)) monomorphic loci:" * "\n $monomorphs")
            println()
        end
    end
    exclude(data, locus = loci_to_rm)
end

"""
    dropmonomorphic!(data::PopData; silent::Bool = false)
Edit a `PopData` object in place by omitting any monomorphic loci. Will inform you which
loci were removed.
"""
function dropmonomorphic!(data::PopData; silent::Bool = false)
    all_loci = loci(data)
    prt = Base.Iterators.partition(data.genodata.genotype, length(samplenames(data)))
    monomorphs = string.(all_loci[[length(unique(i))==1 for i in prt]])
    if length(monomorphs) == 0
        #!silent && println("No monomorphic loci found.\n")
        return data
    elseif !silent
        if length(monomorphs) == 1
            println("Removing monomorphic locus: " * monomorphs[1])
            println()
        else
            println("Removing $(length(monomorphs)) monomorphic loci:" * "\n $monomorphs")
            println()
        end
    end
    exclude!(data, locus = monomorphs)
end

"""
    dropmultiallelic(data::PopData)
Return a `PopData` object omitting loci that are not biallelic.
"""
function dropmultiallelic(data::PopData)
    if data.metadata.biallelic == true
        prinln("PopData object is already biallelic")
        return data
    end
    all_loci = loci(data)
    prt = Base.Iterators.partition(data.genodata.genotype, length(samplenames(data)))
    nonbi = string.(all_loci[[!isbiallelic(i) for i in prt]])
    if length(nonbi) == 1
        @info "Removing 1 multiallelic locus"
        println()
    else
        @info "Removing $(length(nonbi)) multialleic loci"
        println()
    end
    _out = exclude(data, locus = nonbi)
    _out.metadata.biallelic = true
    return _out
end

"""
    dropmultiallelic!(data::PopData)
Edit a `PopData` object in place, removing loci that are not biallelic.
"""
function dropmultiallelic!(data::PopData)
    if data.metadata.biallelic == true
        println("PopData object is already biallelic\n")
        return data
    end
    all_loci = loci(data)
    prt = Base.Iterators.partition(data.genodata.genotype, length(samplenames(data)))
    nonbi = string.(all_loci[[!isbiallelic(i) for i in prt]])
    if length(nonbi) == 1
        @info "Removing 1 multiallelic locus"
        println()
    else
        @info "Removing $(length(nonbi)) multialleic loci"
        println()
    end
    exclude!(data, locus = nonbi)
    data.metadata.biallelic = true
    return data
end

function truncatepath(text::String)
    width = displaysize(stdout)[2]
    if length(text) > width
        separated = split(text, "/")
        newtxt = join(separated[[1,2,3]], "/") * "/…/" * join(separated[[end-1, end]],"/")
        return newtxt
    else
        return text
    end
end


text:
export glscreen

let screen_list = Screen[]
    global current_screen, add_screen, get_screens, empty_screens!
    current_screen() = last(screen_list)
    add_screen(screen) = push!(screen_list, screen)
    get_screens() = copy(screen_list)
    empty_screens!() = empty!(screen_list)
end

function cleanup_old_screens()
    for screen in get_screens()
        destroy!(screen)
    end
    empty!(timer_signal_dict)
    empty_screens!()
    GLFW.Terminate()
    GLFW.Init()
    reset_texture_atlas!()
    GLAbstraction.empty_shader_cache!()
end

function glscreen(name="GLVisualize";
        resolution = GLWindow.standard_screen_resolution(),
        debugging = false,
        background = RGBA(1,1,1,1)
    )

    cleanup_old_screens()

    screen = Screen(name, resolution=resolution, debugging=debugging, color=background)
    add_screen(screen)

    GLWindow.add_complex_signals!(screen) #add the drag events and such
    add_oit_fxaa_postprocessing!(screen) # add postprocessing

    screen
end

const timer_signal_dict = Dict{Int, WeakRef}()
"""
Creates a timer signal with `updates_per_second` while `window` is open.
It's reusing timer signals with the same update rate and registering the updates
with GLFW.
"""
function get_timer_signal(updates_per_second, window=current_screen())
    signal = get!(timer_signal_dict, updates_per_second) do
        # because this is a function, it'll only get executed if needed
        WeakRef(fpswhen(window.inputs[:window_open], updates_per_second))
    end.value
    signal
end

function fold_loop(v0, _)
    val, range, index = v0
    val = range[index]
    index += 1
    index>length(range) && (index = 1)
    (val, range, index)
end

function loop(range::Range, rate=60)
    t = get_timer_signal(rate)
    map(first, foldp(fold_loop, (first(range), range, 1), t))
end

function fold_bounce(v0, _)
    val, range, index, direction = v0
    val = range[index]
    index += direction
    if index in (length(range)+1, 0)
        direction = -direction
        index += 2direction
    end
    (val, range, index, direction)
end

function bounce{T}(range::Range{T}, rate=60)
    t = get_timer_signal(rate)
    map(first, foldp(fold_bounce, (first(range), range, 1, 1), t))
end

export bounce, loop

"""
Materializes the `Visualization` on the an OpenGL `Screen`, transfering everything
to the GPU and adding it to the windows renderlist.
"""
function push!(window::Screen, v::Visualization)
    visual_specific!(v)
    glp = get(v, :gl_primitive, GL_TRIANGLES)
    pre = get(v, :prerender, GLAbstraction.EmptyPrerender())
    # assert opacity is defined
    get!(v, :is_fully_opaque) do
        isopaque(get_color(v))
    end
    if haskey(v, :instances)
        robj = instanced_renderobject(
            v.parameters, shader,
            bb, glp, v[:instances], pre=pre
        )
    else
        robj = std_renderobject(v.parameters, shader, bb, glp, pre=pre)
    end
    if haskey(v, :postrender)
        tmp = robj.postrenderfunction # needs to be always executed
        pr = v[:postrender] # for cleaning up!
        robj.postrenderfunction = () -> (tmp(); pr())
    end
    for child in v.children
        push!(window, child, camera)
    end
    push!(window, robj, v.camera)
end


text:
# # Multi-thread `sum_gather`

include("plots.jl");

# ## Throughput, L1, and last-level (LL) cache misses
#
# (Note: the LL cache miss data may not be available in some machines.)
#
# The same set of data is plotted with two different metrics for the input size
# (x-axis).

# ### x-axis: input size per thread

# The curves for L1 cache miss are virtually identical across benchmarks that
# are run with different number of threads (color, `nthreads`).

plt_throughput_cache_miss_vs_per_thread_size
#-

# ### x-axis: total input size

# The rising part of the curves for LL cache miss for when the total input size
# is large coincide across benchmarks with different number of threads (color,
# `nthreads`). This is where the total input size crosses the LL cache size.

plt_throughput_cache_miss_vs_total_size
#-


text:
using DynamicGrids, DimensionalData, Test, Dates, Unitful,
      CUDAKernels, FileIO, FixedPointNumbers, Colors
using DynamicGrids: Extent, SimData, gridview

if CUDAKernels.CUDA.has_cuda_gpu()
    CUDAKernels.CUDA.allowscalar(false)
    hardware = (SingleCPU(), ThreadedCPU(), CPUGPU(), CuGPU())
else
    hardware = (SingleCPU(), ThreadedCPU(), CPUGPU())
    hardware = (SingleCPU(), ThreadedCPU(), )
end
opts = (NoOpt(), SparseOpt())

# proc = CPUGPU()
proc = SingleCPU()
# opt = SparseOpt()
opt = NoOpt()

# life glider sims

# Test all cycled variants of the array
cyclei!(arrays) = begin
    for A in arrays
        v = A[1, :]
        @inbounds copyto!(A, CartesianIndices((1:size(A, 1)-1, 1:size(A, 2))),
                          A, CartesianIndices((2:size(A, 1), 1:size(A, 2))))
        A[end, :] = v
    end
end

cyclej!(arrays) = begin
    for A in arrays
        v = A[:, 1]
        @inbounds copyto!(A, CartesianIndices((1:size(A, 1), 1:size(A, 2)-1)),
                          A, CartesianIndices((1:size(A, 1), 2:size(A, 2))))
        A[:, end] = v
    end
end

test6_7 = (
    init =  Bool[
             0 0 0 0 0 0 0
             0 0 0 0 1 1 1
             0 0 0 0 0 0 1
             0 0 0 0 0 1 0
             0 0 0 0 0 0 0
             0 0 0 0 0 0 0
            ],
    test2 = Bool[
             0 0 0 0 0 1 0
             0 0 0 0 0 1 1
             0 0 0 0 1 0 1
             0 0 0 0 0 0 0
             0 0 0 0 0 0 0
             0 0 0 0 0 0 0
            ],
    test3 = Bool[
             0 0 0 0 0 1 1
             0 0 0 0 1 0 1
             0 0 0 0 0 0 1
             0 0 0 0 0 0 0
             0 0 0 0 0 0 0
             0 0 0 0 0 0 0
            ],
    test4 = Bool[
             0 0 0 0 0 1 1
             1 0 0 0 0 0 1
             0 0 0 0 0 1 0
             0 0 0 0 0 0 0
             0 0 0 0 0 0 0
             0 0 0 0 0 0 0
            ],
    test5 = Bool[
             1 0 0 0 0 1 1
             1 0 0 0 0 0 0
             0 0 0 0 0 0 1
             0 0 0 0 0 0 0
             0 0 0 0 0 0 0
             0 0 0 0 0 0 0
            ],
    test7 = Bool[
             1 0 0 0 0 1 0
             1 0 0 0 0 0 0
             0 0 0 0 0 0 0
             0 0 0 0 0 0 0
             0 0 0 0 0 0 0
             1 0 0 0 0 0 1
            ]
)

test5_6 = (
    init =  DimArray(Bool[
             0 0 0 0 0 0
             0 0 0 1 1 1
             0 0 0 0 0 1
             0 0 0 0 1 0
             0 0 0 0 0 0
            ], (Y, X)),
    test2 = Bool[
             0 0 0 0 1 0
             0 0 0 0 1 1
             0 0 0 1 0 1
             0 0 0 0 0 0
             0 0 0 0 0 0
            ],
    test3 = Bool[
             0 0 0 0 1 1
             0 0 0 1 0 1
             0 0 0 0 0 1
             0 0 0 0 0 0
             0 0 0 0 0 0
            ],
    test4 = Bool[
             0 0 0 0 1 1
             1 0 0 0 0 1
             0 0 0 0 1 0
             0 0 0 0 0 0
             0 0 0 0 0 0
            ],
    test5 = Bool[
             1 0 0 0 1 1
             1 0 0 0 0 0
             0 0 0 0 0 1
             0 0 0 0 0 0
             0 0 0 0 0 0
            ],
    test7 = Bool[
             1 0 0 0 1 0
             1 0 0 0 0 0
             0 0 0 0 0 0
             0 0 0 0 0 0
             1 0 0 0 0 1
            ]
)

test = test5_6

@testset "Life simulation Wrap" begin
    # Test on two sizes to test half blocks on both axes
    # Loop over shifing init arrays to make sure they all work
    for test in (test5_6, test6_7), i in 1:size(test[:init], 1)
        for j in 1:size(test[:init], 2)
            for proc in hardware, opt in opts
                tspan = Date(2001, 1, 1):Day(2):Date(2001, 1, 14)
                ruleset = Ruleset(;
                    rules=(Life(),),
                    timestep=Day(2),
                    boundary=Wrap(),
                    proc=proc,
                    opt=opt,
                )
                @testset "$(nameof(typeof(proc))) $(nameof(typeof(opt))) results match glider behaviour" begin
                    output = ArrayOutput(test[:init], tspan=tspan)
                    sim!(output, ruleset)
                    @test output[2] == test[:test2]
                    # || (println(2); display(output[2]); display(test[:test2]))
                    @test output[3] == test[:test3] # || (println(3); display(output[3]); display(test[:test3]))
                    @test output[4] == test[:test4] # || (println(4); display(output[4]); display(test[:test4]))
                    @test output[5] == test[:test5] # || (println(5); display(output[5]); display(test[:test5]))
                    @test output[7] == test[:test7] # || (println(7); display(output[7]); display(test[:test7]))
                end
                @testset "$(nameof(typeof(proc))) $(nameof(typeof(opt))) using step!" begin
                    simdata = DynamicGrids._proc_setup(SimData(Extent(; init=test[:init], tspan=tspan), ruleset))
                    # Need Array here to copy from GPU to CPU
                    @test Array(gridview(first(simdata))) == test[:init]
                    simdata = step!(simdata)
                    Array(gridview(first(simdata)))
                    @test Array(gridview(first(simdata))) == test[:test2] || (println("s2"); display(Array(gridview(first(simdata)))); display(test[:test2]))
                    simdata = step!(simdata)
                    @test Array(gridview(first(simdata))) == test[:test3] || (println("s3"); display(Array(gridview(first(simdata)))); display(test[:test3]))
                    simdata = step!(simdata)
                    @test Array(gridview(first(simdata))) == test[:test4] || (println("s4"); display(Array(gridview(first(simdata)))); display(test[:test4]))
                    simdata = step!(simdata)
                    @test Array(gridview(first(simdata))) == test[:test5] || (println("s5"); display(Array(gridview(first(simdata)))); display(test[:test5]))
                    simdata = step!(simdata)
                    simdata = step!(simdata)
                    @test Array(gridview(first(simdata))) == test[:test7] || (println("s7"); display(Array(gridview(first(simdata)))); display(test[:test7]))
                end
            end
            cyclej!(test)
        end
        cyclei!(test)
    end
    nothing
end

@testset "Life simulation with Remove boudary" begin
    init_ =     DimArray(Bool[
                 0 0 0 0 0 0 0
                 0 0 0 0 1 1 1
                 0 0 0 0 0 0 1
                 0 0 0 0 0 1 0
                 0 0 0 0 0 0 0
                 0 0 0 0 0 0 0
                ], (X, Y))
    test2_rem = Bool[
                 0 0 0 0 0 1 0
                 0 0 0 0 0 1 1
                 0 0 0 0 1 0 1
                 0 0 0 0 0 0 0
                 0 0 0 0 0 0 0
                 0 0 0 0 0 0 0
                ]
    test3_rem = Bool[
                 0 0 0 0 0 1 1
                 0 0 0 0 1 0 1
                 0 0 0 0 0 0 1
                 0 0 0 0 0 0 0
                 0 0 0 0 0 0 0
                 0 0 0 0 0 0 0
                ]
    test4_rem = Bool[
                 0 0 0 0 0 1 1
                 0 0 0 0 0 0 1
                 0 0 0 0 0 1 0
                 0 0 0 0 0 0 0
                 0 0 0 0 0 0 0
                 0 0 0 0 0 0 0
                ]
    test5_rem = Bool[
                 0 0 0 0 0 1 1
                 0 0 0 0 0 0 1
                 0 0 0 0 0 0 0
                 0 0 0 0 0 0 0
                 0 0 0 0 0 0 0
                 0 0 0 0 0 0 0
                ]
    test7_rem = Bool[
                 0 0 0 0 0 1 1
                 0 0 0 0 0 1 1
                 0 0 0 0 0 0 0
                 0 0 0 0 0 0 0
                 0 0 0 0 0 0 0
                 0 0 0 0 0 0 0
                ]
    rule = Life{:a,:a}(neighborhood=Moore(1))

    @testset "Wrong timestep throws an error" begin
        rs = Ruleset(rule; timestep=Day(2), boundary=Remove(), opt=NoOpt())
        output = ArrayOutput((a=init_,); tspan=1:7)
        @test_throws ArgumentError sim!(output, rs; tspan=Date(2001, 1, 1):Month(1):Date(2001, 3, 1))
    end

    @testset "Results match glider behaviour" begin
        output = ArrayOutput((a=init_,); tspan=1:7)
        for proc in hardware, opt in opts
            sim!(output, rule; boundary=Remove(), proc=proc, opt=opt)
            output[2][:a]
            output[3][:a]
            output[4][:a]
            output[5][:a]
            output[7][:a]
            @test output[2][:a] == test2_rem
            @test output[3][:a] == test3_rem
            @test output[4][:a] == test4_rem
            @test output[5][:a] == test5_rem
            @test output[7][:a] == test7_rem
        end
    end

    @testset "Combinatoric comparisons in a larger Life sim" begin
        rule = Life(neighborhood=Moore(1))
        init_ = rand(Bool, 100, 99)
        mask_ = ones(Bool, size(init_)...)
        mask_[1:50, 1:50] .= false
        wrap_rs_ref = Ruleset(rule; boundary=Wrap())
        remove_rs_ref = Ruleset(rule; boundary=Remove())
        wrap_output_ref = ArrayOutput(init_; tspan=1:100, mask=mask_)
        remove_output_ref = ArrayOutput(init_; tspan=1:100, mask=mask_)
        sim!(remove_output_ref, remove_rs_ref)
        sim!(wrap_output_ref, wrap_rs_ref)
        for proc in hardware, opt in opts

            @testset "$(nameof(typeof(opt))) $(nameof(typeof(proc)))" begin
                @testset "Wrap" begin
                    wrap_rs = Ruleset(rule; boundary=Wrap(), proc=proc, opt=opt)
                    wrap_output = ArrayOutput(init_; tspan=1:100, mask=mask_)
                    sim!(wrap_output, wrap_rs)
                    wrap_output_ref[2] .- wrap_output[2]
                    @test wrap_output_ref[2] == wrap_output[2]
                    wrap_output_ref[3] .- wrap_output[3]
                    @test wrap_output_ref[3] == wrap_output[3]
                    @test wrap_output_ref[10] == wrap_output[10]
                    @test wrap_output_ref[100] == wrap_output[100]
                end
                @testset "Remove" begin
                    remove_rs = Ruleset(rule; boundary=Remove(), proc=proc, opt=opt)
                    remove_output = ArrayOutput(init_; tspan=1:100, mask=mask_)
                    sim!(remove_output, remove_rs);
                    @test remove_output_ref[2] == remove_output[2]
                    @test remove_output_ref[3] == remove_output[3]
                    remove_output_ref[3] .- remove_output[3]
                    @test remove_output_ref[10] == remove_output[10]
                    @test remove_output_ref[100] == remove_output[100]
                end
            end

        end
    end

end

@testset "sim! with other outputs" begin
    for proc in hardware, opt in opts
        @testset "$(nameof(typeof(opt))) $(nameof(typeof(proc)))" begin
            @testset "Transformed output" begin
                ruleset = Ruleset(Life();
                    timestep=Month(1),
                    boundary=Wrap(),
                    proc=proc,
                    opt=opt,
                )
                tspan_ = Date(2010, 4):Month(1):Date(2010, 7)
                output = TransformedOutput(sum, test6_7[:init]; tspan=tspan_)
                sim!(output, ruleset)
                @test output[1] == sum(test6_7[:init])
                @test output[2] == sum(test6_7[:test2])
                @test output[3] == sum(test6_7[:test3])
                @test output[4] == sum(test6_7[:test4])
            end
            @testset "REPLOutput block works, in Unitful.jl seconds" begin
                ruleset = Ruleset(;
                    rules=(Life(),),
                    timestep=5u"s",
                    boundary=Wrap(),
                    proc=proc,
                    opt=opt,
                )
                output = REPLOutput(test6_7[:init];
                    tspan=0u"s":5u"s":6u"s", style=Block(), fps=1000, store=true
                )
                @test DynamicGrids.isstored(output) == true
                sim!(output, ruleset)
                resume!(output, ruleset; tstop=30u"s")
                @test output[At(5u"s")] == test6_7[:test2]
                @test output[At(10u"s")] == test6_7[:test3]
                @test output[At(20u"s")] == test6_7[:test5]
                @test output[At(30u"s")] == test6_7[:test7]
            end
            @testset "REPLOutput braile works, in Months" begin
                ruleset = Ruleset(Life();
                    timestep=Month(1),
                    boundary=Wrap(),
                    proc=proc,
                    opt=opt,
                )
                tspan_ = Date(2010, 4):Month(1):Date(2010, 7)
                output = REPLOutput(test6_7[:init]; tspan=tspan_, style=Braile(), fps=1000, store=false)
                sim!(output, ruleset)
                @test output[At(Date(2010, 7))] == test6_7[:test4]
                @test DynamicGrids.tspan(output) == Date(2010, 4):Month(1):Date(2010, 7)
                resume!(output, ruleset; tstop=Date(2010, 10))
                @test DynamicGrids.tspan(output) == Date(2010, 4):Month(1):Date(2010, 10)
                @test output[1] == test6_7[:test7]
            end
        end
    end
end

@testset "GifOutput saves" begin
    @testset "Image generator" begin
        # TODO fix on CUDA: cell_to_rgb indexes a CuArray
        ruleset = Ruleset(;
            rules=(Life(),),
            boundary=Wrap(),
            timestep=5u"s",
            opt=NoOpt(),
        )
        output = GifOutput(test6_7[:init];
            filename="test_gifoutput.gif", text=nothing,
            tspan=0u"s":5u"s":30u"s", fps=10, store=true,
        )
        @test output.imageconfig.renderer isa Image
        @test output.imageconfig.textconfig == nothing
        @test DynamicGrids.isstored(output) == true
        sim!(output, ruleset)
        @test output[At(5u"s")] == test6_7[:test2]
        @test output[At(10u"s")] == test6_7[:test3]
        @test output[At(20u"s")] == test6_7[:test5]
        @test output[At(30u"s")] == test6_7[:test7]
        gif = load("test_gifoutput.gif")
        @test gif == RGB.(output.gif)
        rm("test_gifoutput.gif")
    end
    @testset "Layout" begin
        # TODO fix on CUDA
        zeroed = test6_7[:init]
        ruleset = Ruleset(Life{:a}(); boundary=Wrap())
        output = GifOutput((a=test6_7[:init], b=zeroed);
            filename="test_gifoutput2.gif", text=nothing,
            tspan=0u"s":5u"s":30u"s", fps=10, store=true
        )
        @test DynamicGrids.isstored(output) == true
        @test output.imageconfig.renderer isa Layout
        @test output.imageconfig.textconfig == nothing
        sim!(output, ruleset)
        @test all(map(==, output[At(5u"s")], (a=test6_7[:test2], b=zeroed)))
        @test all(map(==, output[At(10u"s")], (a=test6_7[:test3], b=zeroed)))
        @test all(map(==, output[At(20u"s")], (a=test6_7[:test5], b=zeroed)))
        @test all(map(==, output[At(30u"s")], (a=test6_7[:test7], b=zeroed)))
        gif = load("test_gifoutput2.gif")
        @test gif == RGB.(output.gif)
        @test gif[:, 1, 7] == RGB{N0f8}.([1.0, 1.0, 0.298, 0.298, 0.298, 1.0])
        rm("test_gifoutput2.gif")
    end
end

@testset "SparseOpt rules run everywhere with non zero values" begin
    set_hood = SetNeighbors() do data, hood, val, I
        for p in positions(hood, I)
            data[p...] = 2
        end
    end
    clearcell = Cell() do data, val, I
        zero(val)
    end
    output = ArrayOutput(ones(10, 11); tspan=1:3)
    sim!(output, set_hood; opt=SparseOpt())
    @test all(output[3] .=== 2.0)
    sim!(output, set_hood, clearcell; opt=SparseOpt())
    @test all(output[3] .=== 0.0)
end

@testset "Single dimension rules" begin
    init = zeros(Int, 7)
    init[6] = true
    rule110 = Neighbors(Moore(1; ndims=1)) do data, hood, c, I
        l, r = neighbors(hood)
        (c + r + c * r + l * c * r) % 2
    end
    REPLOutput(init; store=false, tspan=1:5, fps=100)
    output = REPLOutput(init; tspan=1:5, fps=100)
    @test DynamicGrids.isstored(output)
    sim!(output, rule110)

    tests = (
        [0, 0, 0, 0, 0, 1, 0,],
        [0, 0, 0, 0, 1, 1, 0,],
        [0, 0, 0, 1, 1, 1, 0,],
        [0, 0, 1, 1, 0, 1, 0,],
        [0, 1, 1, 1, 1, 1, 0,],
    )
    map(tests, output) do t, o
        @test t == o
    end

    init = zeros(Int, 50)
    init[50] = true
    filename = "test_1d_gifoutput.gif"
    output = GifOutput(init; filename=filename, text=nothing, tspan=1:50)
    sim!(output, rule110)
    simgif = load(filename)
    @test simgif == RGB.(output.gif)
    rm(filename)
end


text:
# function RTMimaging(pathsrc, pathrec)
#     fidsrc = open(pathsrc, "r");
#     fidrec = open(pathrec, "r");
#     seek(fidsrc, sizeof(Float64)*2)
#     nz = read(fidsrc, Int64); nx    = read(fidsrc, Int64);
#     ext= read(fidsrc, Int64); iflag = read(fidsrc, Int64);
#     if iflag == 1
#        Nz = nz + 2*ext
#     elseif iflag == 2
#        Nz = nz +   ext
#     end
#     Nx = nz + 2*ext
#     head = sizeof(Int64)*8; sptsize = sizeof(Float64)*Nz*Nx*4
#     nt = round(Int64, (filesize(fidsrc)-head)/sptsize)
#     Ip = zeros(Nz, Nx)
#     head = sizeof(Int64)*8 + sizeof(Float64)*Nz*Nx*2
#     for it = 1 : nt
#         position = head + (it-1)*sptsize
#         seek(fidsrc, position); seek(fidrec, position)
#         psrc = reshape(read(fidsrc, Float64, Nz*Nx), Nz, Nx) + reshape(read(fidsrc, Float64, Nz*Nx), Nz, Nx)
#         prec = reshape(read(fidrec, Float64, Nz*Nx), Nz, Nx) + reshape(read(fidrec, Float64, Nz*Nx), Nz, Nx)
#         Ip = Ip + psrc .* prec
#     end
#     close(fidsrc); close(fidrec);
#     return Ip
# end

function RTMimaging(pathsrc, pathrec; image_type="p")
    fidsrc = open(pathsrc, "r");
    fidrec = open(pathrec, "r");
    seek(fidsrc, sizeof(Float64)*2)
    nz = read(fidsrc, Int64); nx    = read(fidsrc, Int64);
    ext= read(fidsrc, Int64); iflag = read(fidsrc, Int64);
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
    end
    Nx = nz + 2*ext
    head = sizeof(Int64)*8; sptsize = sizeof(Float64)*Nz*Nx*4
    nt = round(Int64, (filesize(fidsrc)-head)/sptsize)
    if image_type == "p"
       Ip = zeros(Nz, Nx);
       head = sizeof(Int64)*8 + sizeof(Float64)*Nz*Nx*2
       for it = 1 : nt
           position = head + (it-1)*sptsize
           seek(fidsrc, position); seek(fidrec, position)
           psrc = reshape(read(fidsrc, Float64, Nz*Nx), Nz, Nx) + reshape(read(fidsrc, Float64, Nz*Nx), Nz, Nx)
           prec = reshape(read(fidrec, Float64, Nz*Nx), Nz, Nx) + reshape(read(fidrec, Float64, Nz*Nx), Nz, Nx)
           Ip = Ip + psrc .* prec
       end
       close(fidsrc); close(fidrec);
       return Ip
    elseif image_type == "all"
       Ivz = zeros(Nz, Nx); Ivx = zeros(Nz, Nx);
       Ipz = zeros(Nz, Nx); Ipx = zeros(Nz, Nx);
       for it = 1 : nt
           position = head + (it-1)*sptsize
           seek(fidsrc, position); seek(fidrec, position)
           Svz = reshape(read(fidsrc, Float64, Nz*Nx), Nz, Nx)
           Svx = reshape(read(fidsrc, Float64, Nz*Nx), Nz, Nx)
           Spz = reshape(read(fidsrc, Float64, Nz*Nx), Nz, Nx)
           Spx = reshape(read(fidsrc, Float64, Nz*Nx), Nz, Nx)
           Rvz = reshape(read(fidrec, Float64, Nz*Nx), Nz, Nx)
           Rvx = reshape(read(fidrec, Float64, Nz*Nx), Nz, Nx)
           Rpz = reshape(read(fidrec, Float64, Nz*Nx), Nz, Nx)
           Rpx = reshape(read(fidrec, Float64, Nz*Nx), Nz, Nx)
           Ivz = Ivz + Svz .* Rvz
           Ivx = Ivx + Svx .* Rvx
           Ipz = Ipz + Spz .* Rpz
           Ipx = Ipx + Spx .* Rpx
       end
       close(fidsrc); close(fidrec);
       return Ivz, Ivx, Ipz, Ipx
    end
end

max_stars_repo_name: WneleiGao/AcousticWaveLab

max_issues_repo_name: WneleiGao/AcousticWaveLab

max_forks_repo_name: WneleiGao/AcousticWaveLab


text:
using BenchmarkConfigSweeps

if Sys.CPU_THREADS < 8
    threads = [1, Sys.CPU_THREADS]
else
    threads = [1, 4, 8]
end

BenchmarkConfigSweeps.run(
    joinpath(@__DIR__, "build"),
    joinpath(@__DIR__, "benchmarks.jl"),
    BenchmarkConfigSweeps.nthreads.(threads),
)


text:
# This is Julia code for discretizing AR(1) process by using Rouwenhorst method (1995)
# written on 2017/10/17
#
# Model:
# z' = rho*z+e  e~N(0,(sig_e)^2)
#
# INPUTS
# N:     Number of grids
# rho:   Coeffecient of AR(1) process z'=mew+rho*z+e, e~N(0,sigmasq)
# sig_e: Std of innovation, e, in AR(1) process
#
# OUTPUT
# vZ:    grids
# mPI:   transition prob. matrix

#########################################################################
## 1. Parameters
#########################################################################

# Ex. Standard quarterly calibration for productivity process
# of US economy

const N = 5;            # Number of grids, N
const rho = 0.95;       # AR(1) coefficient
const sig_e = 0.007;    # std of innovation

#########################################################################
## 2. Generating grids && Parameterization
#########################################################################

sig_z = sqrt(sig_e^2/(1-rho^2));    # std of AR(1) r.v.
z_max = sig_z*sqrt(N-1);            # max grid
z_min = -z_max;                     # min grid
p = (1+rho)/2;                      # initial prob
q = p;

vZ = linspace(z_min, z_max, N)      # Grids

#########################################################################
## 3. Transition matrix
#########################################################################

mPI_old = [p 1-p; 1-q q]; # Initial PI
mPI_new = zeros(3,3);

for i in 3:N
    mPI_new  = zeros(i,i);
    mPI_new1 = zeros(i,i);  mPI_new1[1:i-1,1:i-1]    = mPI_old;
    mPI_new2 = zeros(i,i);  mPI_new2[1:i-1,2:i]      = mPI_old;
    mPI_new3 = zeros(i,i);  mPI_new3[2:i,1:i-1]      = mPI_old;
    mPI_new4 = zeros(i,i);  mPI_new4[2:i,2:i]        = mPI_old;

    mPI_new = p*mPI_new1 + (1-p)*mPI_new2 + (1-q)*mPI_new3 + q*mPI_new4;
    mPI_new[2:i-1,:] = 0.5*mPI_new[2:i-1,:]; # devide by 2

    mPI_old = mPI_new;
end

###########################################################
## 4. Display result with 4 decimals
###########################################################

print("Grid points are \n")
for i in 1:N
    @printf "%6.4f" vZ[i]
    print(" ")
end
print("\n")
print("Transition prob. matrix is\n")
for i in 1:N
    for j in 1:N
        @printf "%6.4f" mPI_new[i,j]
        print(" ")
    end
    print("\n")
end


text:
using Base.Test
# fourierlayer.jl
import WaveletScattering: FourierLayer

fourier = FourierLayer(
        bottoms = [:signal],
        tops = [:fourier],
        pathkeys = [PathKey(:time)])

@test !Mocha.can_do_bp(fourier)
@test !Mocha.has_neuron(fourier)
@test !Mocha.has_param(fourier)


text:
function hypotour_inequalities!(
    cut_manager::CutManager,
    demand::Vector{Int},
    capacity::Int,
    edge_tail::Vector{Int},
    edge_head::Vector{Int},
    edge_x::Vector{Float64};
    max_n_cuts = 1000::Cint
)
    n_customers, _demand, n_edges, _edge_tail, _edge_head, _edge_x = input_conversion(demand, edge_tail, edge_head, edge_x)

    new_cmp = init_CMP()

    _max_violation = Ref{Cdouble}()

    ccall(
        (:HTOURSEP_SeparateHTours, LIB_CVRPSEP),
        Cvoid,
        (
            Cint, Ptr{Cint}, Cint,
            Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
            CnstrMgrPointer, Cint,
            Ref{Cdouble}, CnstrMgrPointer
        ),
        Cint(n_customers), Cint.(_demand), Cint(capacity),
        Cint(n_edges), Cint.(_edge_tail), Cint.(_edge_head), Cdouble.(_edge_x),
        cut_manager.cmp, Cint(max_n_cuts),
        _max_violation, new_cmp
    )

    cut_manager.new_cuts = convert_CMP(new_cmp)
    # jl_CMGR_MoveCnstr!(new_cmp, cut_manager.cmp)
    cut_manager.violation = Float64(_max_violation[])

    n_new_cuts =  length(cut_manager.new_cuts)
    Tail = Vector{Vector{Int}}(undef, n_new_cuts)
    Head = Vector{Vector{Int}}(undef, n_new_cuts)
    Coeff = Vector{Vector{Float64}}(undef, n_new_cuts)
    RHS = Vector{Float64}(undef, n_new_cuts)

    for j in 1:n_new_cuts
        cr = cut_manager.new_cuts[j]
        @assert cr.Ctype == CMGR_CT_TWOEDGES_HYPOTOUR
        Tail[j] = cr.IntList .+ 1
        Head[j] = cr.ExtList .+ 1
        Coeff[j] = cr.CoeffList
        RHS[j] = cr.RHS

        @show Tail[j], Head[j], Coeff[j], RHS[j]
    end

    return Tail, Head, Coeff, RHS
end


text:
"""
Two-dimensional NEP
"""
module Spectral2D

using ..Common
using SparseArrays
using NonlinearEigenproblems

import ..HelmholtzNEP
import ..MaxwellNEP

"""
    HelmholtzNEP(sim) -> nep

`sim` is a `Simulation`
"""
function HelmholtzNEP(
            sim::Simulation{1,C,T};
            kwargs...
            ) where {C,T}

    ka = 0
    M = Helmholtz(sim)
    Σ0,Σ1,Σ2 = M.Σs
    f = M.sim.self_energy.f
    # f = map(f->(ω->f(ω)),M.sim.self_energy.f)
    Fs, fχs = _compute_Fsfχs(sim,1)
    As = vcat([Σ0,Σ1,Σ2,sim.laplacian.l0,M.αε],Fs)
    fs = [f[1],f[2],one,one,ω->ω^2,map(ϝ->(ω->ω^2*ϝ(ω)),fχs)...]
    nep = SPMF_NEP(As,fs;kwargs...)
    return HelmholtzNEP(M,nep,6:length(As))
end

"""
    MaxwellNEP(::Simulation{1}; ky=0, kz=0) -> nep
"""
function MaxwellNEP(
            sim::Simulation{1,C,T};
            ky::Number = 0,
            kz::Number = 0,
            kwargs...
            ) where {C,T}

    ka = 0
    M = Maxwell(sim; ky=ky, kz=kz)
    f = map(f->(ω->f(ω,ka,ky,kz)),M.sim.self_energy.f)
    Σ1,Σ2,Σ3 = M.Σs
    Fs, fχs = _compute_Fsfχs(sim,3)
    As = vcat([Σ1,Σ2,Σ3,sim.curlcurl(ky,kz),M.αε],Fs)
    fs = [f[1],f[2],one,one,ω->-ω^2,map(ϝ->(ω->-ω^2*ϝ(ω)),fχs)...]
    nep = SPMF_NEP(As,fs;kwargs...)
    return MaxwellNEP(M,nep,6:length(As))
end

################################################################################

function _compute_Fsfχs(sim::Simulation{1},m::Integer)
    χs = map(d->d.χ,sim.dispersive_domains)
    fχs = map(x->(ω->susceptability(x,ω)),χs)
    Fs = Vector{SparseMatrixCSC{ComplexF64,Int}}(undef,length(sim.dispersive_domains))
    for d ∈ eachindex(Fs) Fs[d] = spdiagm(0=>repeat(sim.Fs[d],m)) end
    return Fs,fχs
end

end # module

using .Spectral2D


text:
function symbol_and_type(ex::Expr)
	if ex.head == :(::) && isa(ex.args[1],Symbol)
		return ex.args[1], ex.args[2]
	else
		error("Could not recognize symbol with type: $ex")
	end
end
symbol_and_type(s::Symbol) = s, Float64

function is_new(s::Symbol,arr::Vector)
	reduce(&,arr.!=s)
end

function parse_argument_data(ex::Expr)
	if ex.head != :tuple
		error("Expression of unsupported format: $ex")
	end
	new_var_row = reshape([nothing,-Inf,Inf,false,false,nothing],(1,6))
	vars = new_var_row
	for p in ex.args[1:end-1]
		if length(p.args) == 5 ## a<x<b or a>x>b
			sym, typ = symbol_and_type(p.args[3])
			if is_new(sym, vars[2:end,1])
				vars = vcat(vars, new_var_row)
				ind = size(vars,1)
			else
				ind = find(vars[:,1],sym)[1]
			end
			vars[ind,1] = sym
			vars[ind,4] = true
			vars[ind,5] = true
			vars[ind,6] = typ
			if p.args[2] == :< && p.args[4] == :<     ## a<x<b
				vars[ind,2] = p.args[1]
				vars[ind,3] = p.args[5]
			elseif p.args[2] == :> && p.args[4] == :> ## a>x>b
				vars[ind,2] = p.args[5]
				vars[ind,3] = p.args[1]
			else
				error("Could not parse argument data $p")
			end

		elseif length(p.args) == 3 ## a<<x<b, a>>x>b, a<x<<b or a>x>>b
			if p.args[1] == :<

				if isa(p.args[2], Expr) &&		## a<<x<b
				   length(p.args[2].args) == 3 && p.args[2].args[1] == :<<
					sym, typ = symbol_and_type(p.args[2].args[3])
					if is_new(sym, vars[2:end,1])
						vars = vcat(vars, new_var_row)
						ind = size(vars,1)
					else
						ind = find(vars[:,1],sym)[1]
					end
					vars[ind,1] = sym
					vars[ind,2] = p.args[2].args[2]
					vars[ind,3] = p.args[3]
					vars[ind,4] = false
					vars[ind,5] = true
					vars[ind,6] = typ

				elseif isa(p.args[3], Expr) &&		## a<x<<b
				   length(p.args[3].args) == 3 && p.args[3].args[1] == :<<
					sym, typ = symbol_and_type(p.args[3].args[2])
					if is_new(sym, vars[2:end,1])
						vars = vcat(vars, new_var_row)
						ind = size(vars,1)
					else
						ind = find(vars[:,1],sym)[1]
					end
					vars[ind,1] = sym
					vars[ind,2] = p.args[2]
					vars[ind,3] = p.args[3].args[3]
					vars[ind,4] = true
					vars[ind,5] = false
					vars[ind,6] = typ
				else
					error("Could not parse argument data $p")
				end

			elseif p.args[1] == :>
				if isa(p.args[2], Expr) &&		## a>>x>b
				   length(p.args[2].args) == 3 && p.args[2].args[1] == :>>
					sym, typ = symbol_and_type(p.args[2].args[3])
					if is_new(sym, vars[2:end,1])
						vars = vcat(vars, new_var_row)
						ind = size(vars,1)
					else
						ind = find(vars[:,1],sym)[1]
					end
					vars[ind,1] = sym
					vars[ind,2] = p.args[3]
					vars[ind,3] = p.args[2].args[2]
					vars[ind,4] = true
					vars[ind,5] = false
					vars[ind,6] = typ
				elseif isa(p.args[3], Expr) &&	## a>x>>b
				   length(p.args[3].args) == 3 && p.args[3].args[1] == :>>
					sym, typ = symbol_and_type(p.args[3].args[2])
					if is_new(sym, vars[2:end,1])
						vars = vcat(vars, new_var_row)
						ind = size(vars,1)
					else
						ind = find(vars[:,1],sym)[1]
					end
					vars[ind,1] = sym
					vars[ind,2] = p.args[3].args[3]
					vars[ind,3] = p.args[2]
					vars[ind,4] = false
					vars[ind,5] = true
					vars[ind,6] = typ
				else
					error("Could not parse argument data $p")
				end
			else
				error("Could not parse argument data $p")
			end
		else
			error("Could not parse argument data $p")
		end
	end
	vars = vars[2:end,:]
	vars
end


text:
# ===== VSB Verifications: RunTime_Function Building =====

println(" === BEGINNING Verifications3.jl === ")

# === Imports ===
println("   IMPORTING PACKAGES...")
using VSB
using SimpleVPM
using LinearAlgebra
using Printf

cd("/Users/Damyn/Documents/BYU/FLOW Lab/VSB/verifications")

# === Boundary Geometry ===
println("   GENERATING GEOMETRY...")
# NPTS = 20                               # Number of Points on Boundary
# del_theta = (2*pi)/(NPTS    )           # Step in theta around body
# R = 1.0                                 # Cylinder Radius
# k_hat = [0.0, 0.0, 1.0]                 # Z direction unit vector
# X_coor = zeros(NPTS)                    # X-coordinates
# Y_coor = zeros(NPTS)                    # Y-coordinates
# X_tHat = zeros(NPTS)                    # X Tangent Vector
# Y_tHat = zeros(NPTS)                    # X Tangent Vector
# X_nHat = zeros(NPTS)                    # Y Normal Vector
# Y_nHat = zeros(NPTS)                    # Y Normal Vector
# overlap = 2.0                           # Particle overlap
#
# phase = 0.0 * del_theta
# for i=1:NPTS
#     X_coor[i] = R*cos((i-1)*del_theta + phase)
#     Y_coor[i] = R*sin((i-1)*del_theta + phase)
#
#     X_tHat[i] = sin((i-1)*del_theta + phase)
#     Y_tHat[i] = -cos((i-1)*del_theta + phase)
#
#     X_nHat[i] = cos((i-1)*del_theta + phase)
#     Y_nHat[i] = sin((i-1)*del_theta + phase)
# end
# body_pts = [[X_coor[i], Y_coor[i], 0.0] for i in 1:NPTS]
# t_hats = [[X_tHat[i], Y_tHat[i], 0.0] for i in 1:NPTS]
# n_hats = [[X_nHat[i], Y_nHat[i], 0.0] for i in 1:NPTS]

NPTS = 20
body = VSB.NACA4("1200", NPTS)

sigma = overlap * del_theta             # Guassian Spreading

# === Parameters ===
# Run 5
nu = .1                                     # Kinematic Viscosity
magU_inf = 1.0                             # Magnitude of Free-Stream Velocity
del_t = 0.5                                 # Time step
NSTEPS = 300                                # Number of time steps

U_inf(t) = magU_inf * [1.0, 0.0, 0.0]       # Free-Stream velocity

# === Set up SimpleVPM ===
println("   CREATING PARTICLE FIELD...")
pfield = SimpleVPM.ParticleField(nu,
                                 zeta_gauserf,                # Basis function
                                 g_gauserf;                   # Regularizing function of zeta
                                 transposed=true,             # Transposed scheme
                                 relax=!true,                 # Pedrizzeti relaxation scheme
                                 rlxf=0.3,                    # Relaxation factor
                                 integration="euler",         # Time integration scheme
                                 Uinf=U_inf)                  # Free-Stream velocity function

# === Velocity function ===
function U_slip(X)
    if length(pfield.particles) == 0
        return U_inf(0)
    else
        return SimpleVPM.U(pfield, X) + U_inf(0)
    end
end

# === Initial vorticity field ===
params = VSB.Parameters(body, U_slip, verbose = false, sigma=sigma)
alphas = params.alphas
betas = VSB.CalcVSDiffusionCoefs(body, alphas, del_t, pfield.nu, sigma=sigma)
betas[1] = 0.0
betas[NPTS] = 0.0

offset = 0.01             # Normal direction offset for new particle locations
for i in 1:NPTS
    loc = body_pts[i] + offset .* n_hats[i]
    Gamma = (1/sigma^3) * betas[i] .* k_hat
    vol = 1.0
    particle = SimpleVPM.Particle(loc, Gamma, sigma, vol)
    SimpleVPM.addparticle(pfield, particle)
end

# === RunTimeFunction ===
function RunTimeFunction(pfield, t, del_t)

    if (t != 0)

        for p in pfield.particles
            p.Gamma[1] = 0.0
            p.Gamma[2] = 0.0
        end

        # === Calculate parameters ===
        params = VSB.Parameters(body, U_slip, verbose = false)
        alphas = params.alphas
        betas = VSB.CalcVSDiffusionCoefs(body, alphas, del_t, pfield.nu)

        # === Add particles on surface ===
        offset = 0             # Normal direction offset for new particle locations
        for i in 1:NPTS
            loc = body_pts[i] + offset .* n_hats[i]
            Gamma = (1/sigma^3) * betas[i] .* k_hat
            vol = 1.0
            particle = SimpleVPM.Particle(loc, Gamma, sigma, vol)
            SimpleVPM.addparticle(pfield, particle)
        end

    end

    # === Returns ===
    return false

end

# === Simulation Parameters ===
println("   RUNNING SIMULATION...")

NSTEPS_relax = -1
sigma0 = sigma
file_path = splitdir(@__FILE__)[1]
save_path = joinpath(file_path, "sims/verification3_07/")
run_name = "verifications3_07"
paraview = true

SimpleVPM.run_vpm!(pfield, del_t, NSTEPS,
                   runtime_function = RunTimeFunction,
                   nsteps_relax = NSTEPS_relax,
                   rbf_itmax = 500,
                   sgm0 = sigma0,
                   save_path = save_path,
                   run_name = run_name)

if paraview
    println("   CALLING PARAVIEW...")
    strn = run_name * "...xmf;"
    run(`paraview --data="$save_path$strn"`)
end

println(" === END OF Verifications3.jl === ")


text:
"""
    asquantiles!(x, n)
"""
function asquantiles!(x::AbstractVector, n::Int)
    quants = nquantile(x, n)[1:end-1]
    for (i,j) in enumerate(x)
        x[i] = searchsortedlast(quants, j)
    end
end

"""
    asquantiles(x, n)
"""
function asquantiles(x::AbstractVector, n::Int)
    quants = nquantile(x, n)[1:end-1]
    ret = similar(x)
    for (i,j) in enumerate(x)
        ret[i] = searchsortedlast(quants, j)
    end
    ret
end

eltypest(df) = map(x->Base.nonmissingtype(x), eltype.(eachcol(df)))
eltypet(x) = Base.nonmissingtype(eltype(x))


text:
struct KalmanFilter <: Transducer
    lssm::LinearStateSpaceModel
    lobs::LinearObservation
end

function predict(l::LinearStateSpaceModel,m,P)
    l.A*m, l.A*P*l.A'+l.Q
end

function update(l::LinearObservation,y,m,P)
    v = y-l.H*m
    S = l.H*P*l.H' + l.R
    K = P*l.H'*inv(S)
    mu = m + K*v
    Pu = P - K*S*K'
    return mu, Pu, v, S
end

function kalman(lssm::LinearStateSpaceModel, lo::LinearObservation,y,m,P)
    m_i,P_i = predict(lssm,m,P)
    update(lo,y,m_i,P_i)
end

function Transducers.start(rf::R_{KalmanFilter}, result)
    private_state = (zeros(size(xform(rf).lssm.A,1)), one(xform(rf).lssm.A))
    return wrap(rf, private_state, start(inner(rf), result))
end

function Transducers.next(rf::R_{KalmanFilter}, result, input)
    wrapping(rf, result) do st, iresult
        (m, P) = st
        mu, Pu, v, S = kalman(xform(rf).lssm, xform(rf).lobs, input, m, P)
        iresult = next(inner(rf), iresult, (mu, Pu, v, S))
        return (mu, Pu), iresult
    end
end


text:
# using Base: Int64, Float64
# include("../src/meta-icvi/rocket.jl")
using MetaICVI

n_kernels = 5
input_length = 100

rocket = RocketModule(input_length, n_kernels)

data = rand(input_length)
apply_kernel(rocket.kernels[1], data)
n_example = 5
t = zeros(n_example, 2)
for i = 1:n_example
    t[i, :] = apply_kernel(rocket.kernels[i], data)
end

features = apply_kernels(rocket, data)

@info "Sequential:" typeof(t) size(t) t
@info "Batch:" typeof(features) size(features) features

@test features == t


text:
function test_cgls()
  cgls_tol = 1.0e-5;

  for npower = 1 : 4
    (b, A, D, HY, HZ, Acond, rnorm) = test(40, 40, 4, npower, 0);  # No regularization.

    (x, stats) = cgls(A, b);
    resid = norm(A' * (A*x - b)) / norm(b)
    @printf("CGLS: Relative residual: %8.1e\n", resid);
    @test(resid <= cgls_tol);
    @test(stats.solved);

    λ = 1.0e-3;
    (x, stats) = cgls(A, b, λ=λ);
    resid = norm(A' * (A*x - b) + λ * x) / norm(b)
    @printf("CGLS: Relative residual: %8.1e\n", resid);
    @test(resid <= cgls_tol);
    @test(stats.solved);
  end

  # Test with preconditioning.
  A = rand(10, 6); b = rand(10);
  M = InverseLBFGSOperator(10, 4);
  for _ = 1 : 6
    s = rand(10);
    y = rand(10);
    push!(M, s, y);
  end

  (x, stats) = cgls(A, b, M=M);
  resid = norm(A' * M * (A * x - b)) / sqrt(dot(b, M * b));
  @printf("CGLS: Preconditioned residual: %8.1e\n", resid);
  @test resid <= cgls_tol;

  # test trust-region constraint
  (x, stats) = cgls(A, b)

  radius = 0.75 * norm(x)
  (x, stats) = cgls(A, b, radius=radius)
  @test(stats.solved);
  @test(abs(radius - norm(x)) <= cgls_tol * radius)

  opA = LinearOperator(A)
  (xop, statsop) = cgls(opA, b, radius=radius)
  @test(abs(radius - norm(xop)) <= cgls_tol * radius)

  # Code coverage.
  (b, A, D, HY, HZ, Acond, rnorm) = test(40, 40, 4, 3, 0);
  (x, stats) = cgls(Matrix(A), b);
  (x, stats) = cgls(sparse(Matrix(A)), b);
  show(stats);

  # Test b == 0
  (x, stats) = cgls(A, zeros(size(A,1)))
  @test x == zeros(size(A,1))
  @test stats.status == "x = 0 is a zero-residual solution"

  # Test integer values
  A = [I; rand(1:10, 2, 3)]
  b = A * ones(Int, 3)
  (x, stats) = cgls(A, b)
  @test stats.solved
end

test_cgls()


text:
module LibYAML

using Libdl

# Load `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("LibYAML not installed properly, run Pkg.build(\"LibYAML\"), restart Julia, and try again")
end
include(depsjl_path)

# Module initialization function
function __init__()
    # Always check your dependencies from `deps.jl`
    check_deps()
end

# Version Information

"""
Get the library version as a string.

The function returns the pointer to a static string of the form
`X.Y.Z`, where `X` is the major version number, `Y` is a minor version
number, and `Z` is the patch version number.
"""
function get_version_string()::String
    unsafe_string(ccall((:yaml_get_version_string, libyaml), Cstring, ()))
end

"""
Get the library version numbers.

Returns a tuple `(major, minor, patch)`.
"""
function get_version()::NTuple{3, Int}
    major = Ref{Cint}(0)
    minor = Ref{Cint}(0)
    patch = Ref{Cint}(0)
    ccall((:yaml_get_version, libyaml),
          Cvoid, (Ref{Cint}, Ref{Cint}, Ref{Cint}), major, minor, patch)
    major[], minor[], patch[]
end

# Basic Types

"""The version directive data."""
struct VersionDirective
    """The major version number."""
    major::Int
    """The minor version number."""
    minor::Int
end

"""The tag directive data."""
struct TagDirective
    """The tag handle."""
    handle::Cstring
    """The tag prefix."""
    prefix::Cstring
end

# Tokens

# Events

# Nodes

# NOTE: It turns out that wrapping LibYAML directly in Julia is not
# feasible since LibYAML uses unions in its types. It is unclear to me
# how to create respective Julia types that have the correct size. It
# might be necessary to create C wrappers that allocate or free
# LibYAML's objects, and provide Julia accessors that don't use C
# unions.

"""The document nodes."""
struct DocumentNodes
#     """The beginning of the stack."""
#     start::Ref{Node}
#     """The end of the stack."""
#     end_::Ref{Node}
#     """The top of the stack."""
#     top::Ref{Node}
end

"""The document structure."""
struct Document
    """The document nodes."""
    nodes::DocumentNodes

    """The version directive."""
    version_directive::Ref{VersionDirective}

#     /** The list of tag directives. */
#     struct {
#         /** The beginning of the tag directives list. */
#         yaml_tag_directive_t *start;
#         /** The end of the tag directives list. */
#         yaml_tag_directive_t *end;
#     } tag_directives;
#
#     /** Is the document start indicator implicit? */
#     int start_implicit;
#     /** Is the document end indicator implicit? */
#     int end_implicit;
#
#     /** The beginning of the document. */
#     yaml_mark_t start_mark;
#     /** The end of the document. */
#     yaml_mark_t end_mark;
end

# Parser definitions

struct Parser
    parser::Ptr{Cvoid}
end

"""
Parse the input stream and produce the next YAML document.

Call this function subsequently to produce a sequence of documents
constituting the input stream.

If the produced document has no root node, it means that the document
end has been reached.

An application is responsible for freeing any data associated with the
produced document object using the `document_delete()` function.

An application must not alternate the calls of `parser_load()` with the
calls of `parser_scan()` or `parser_parse()`. Doing this will break
the parser.

@param[in,out]   parser      A parser object.
@param[out]      document    An empty document object.

@return @c 1 if the function succeeded, @c 0 on error.
"""
function parser_load(parser::Parser)::Document
    doc = Ref{Document}(Document())
    ccall((:yaml_parser_load, libyaml),
          Cint, (Ptr{Cvoid}, Ref{Document}), parser.parser, doc)
    doc[]
end

# Emitter definitions

end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to the conversion from small Euler angles to DCM.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export smallangle_to_dcm

"""
    smallangle_to_dcm(θx::Number, θy::Number, θz::Number; normalize = true)

Create a direction cosine matrix from three small rotations of angles `θx`,
`θy`, and `θz` [rad] about the axes X, Y, and Z, respectively.

If the keyword `normalize` is `true`, then the matrix will be normalized using
the function `orthonormalize`.

# Example

```jldoctest
julia> smallangle_to_dcm(+0.01, -0.01, +0.01)
DCM{Float64}:
  0.9999     0.00989903  0.010098
 -0.009999   0.999901    0.00989802
 -0.009999  -0.009998    0.9999

julia> smallangle_to_dcm(+0.01, -0.01, +0.01; normalize = false)
DCM{Float64}:
  1.0    0.01  0.01
 -0.01   1.0   0.01
 -0.01  -0.01  1.0
```
"""
@inline function smallangle_to_dcm(
    θx::T1,
    θy::T2,
    θz::T3;
    normalize = true
) where {T1<:Number, T2<:Number, T3<:Number}
    # Since we might orthonormalize `D`, we need to get the float to avoid type
    # instabilities.
    T = float(promote_type(T1, T2, T3))

    D = DCM{T}(
          1, +θz, -θy,
        -θz,   1, +θx,
        +θy, -θx,   1
    )'

    if normalize
        return orthonormalize(D)
    else
        return D
    end
end


text:
using Distributions
using StaticArrays
import Printf.@sprintf
import Distributions.RealInterval

abstract type AbstractUncertainValue end

Broadcast.broadcastable(d::AbstractUncertainValue) = Ref(d)

##############################
# Support of an observation
##############################
import Distributions.support

"""
    support(o::AbstractUncertainValue) -> IntervalArithmetic.Interval

Find the support of the uncertain observation.
"""
function support(o::AbstractUncertainValue)
    supp = support(o.distribution)
    lowerbound = supp.lb
    upperbound = supp.ub
    @interval(lowerbound, upperbound)
end

export support

##############################
# Intersection of the supports
##############################
"""
    support_overlap(o1::AbstractUncertainValue, o2::AbstractUncertainValue)

Compute the overlap in the supports of two uncertain observations.
"""
function support_overlap(uval1::AbstractUncertainValue, uval2::AbstractUncertainValue)
    s1 = support(uval1)
    s2 = support(uval2)

    if s1 isa RealInterval
        s1 = interval(s1...)
    end

    if s2 isa RealInterval
        s2 = interval(s2...)
    end
    s1 ∩ s2
end

export support_overlap

##################################################################
# Intersection of two UncertainValues as a mixture model
##################################################################
import Base.intersect, Base.∩

"""
    intersect(o1::AbstractUncertainValue, o2::AbstractUncertainValue)

Compute the intersection between two uncertain observations probabilistically.
The intersection is represented as a mixture model of the distributions
furnishing the observations.
"""
function intersect(o1::AbstractUncertainValue, o2::AbstractUncertainValue)
    # Create a mixture model representing the intersection
    if support_overlap(o1, o2) == ∅
        throw(DomainError((o1, o2), "intersect(o1, o2) == ∅. Cannot compute mixture model."))
    end

    MixtureModel([o1.distribution, o2.distribution])
end

export intersect

###################
# Pretty printing
###################

function summarise(o::AbstractUncertainValue)
    _type = typeof(o)
    "$_type"
end
Base.show(io::IO, q::AbstractUncertainValue) = print(io, summarise(q))

dimension(usv::AbstractUncertainValue) = 1

Base.max(uv::AbstractUncertainValue) = maximum(uv.distribution)
Base.min(uv::AbstractUncertainValue) = minimum(uv.distribution)
Base.minimum(uv::AbstractUncertainValue) = minimum(uv.distribution)
Base.maximum(uv::AbstractUncertainValue) = maximum(uv.distribution)

Base.maximum(uvs::AbstractVector{AbstractUncertainValue}) =
    maximum([maximum(uv) for uv in uvs])
Base.minimum(uvs::AbstractVector{AbstractUncertainValue}) =
    minimum([minimum(uv) for uv in uvs])

export
AbstractUncertainValue


text:
# ==================
# Genetic algorithms
# ==================

"""
    identity(v1, v2)

Returns the same parameter individuals `v1` and `v2` as an offspring pair.
"""
identity(v1::T, v2::T; kwargs...) where {T<:AbstractVector} = (v1, v2)
identity(v1::T, v2::T; kwargs...) where {T<:AbstractIndividual} = (v1, v2)

# Binary crossovers
# -----------------

"""
    SPX(v1, v2)

Single point crossover between `v1` and `v2` individuals.
"""
function SPX(
    v1::T,
    v2::T;
    rng::AbstractRNG = Random.default_rng(),
) where {T<:AbstractVector}
    l = length(v1)
    c1 = copy(v1)
    c2 = copy(v2)
    pos = rand(rng, 1:l)
    for i = pos:l
        vswap!(c1, c2, i)
    end
    return c1, c2
end
SPX(v1::T, v2::T; rng::AbstractRNG = Random.default_rng()) where {T<:AbstractIndividual} =
    map(T, SPX(variables(v1), variables(v2), rng=rng))

"""
    TPX(v1, v2)

Two point crossover between `v1` and `v2` individuals.
"""
function TPX(
    v1::T,
    v2::T;
    rng::AbstractRNG = Random.default_rng(),
) where {T<:AbstractVector}
    l = length(v1)
    c1 = copy(v1)
    c2 = copy(v2)
    from, to = rand(rng, 1:l, 2)
    from, to = from > to ? (to, from) : (from, to)
    for i = from:to
        vswap!(c1, c2, i)
    end
    return c1, c2
end
TPX(v1::T, v2::T; rng::AbstractRNG = Random.default_rng()) where {T<:AbstractIndividual} =
    map(T, TPX(variables(v1), variables(v2), rng=rng))

"""
    SHFX(v1, v2)

Shuffle crossover between the parents `v1` and `v2` that performs recombination similar to [`SPX`](@ref) preliminary shuffling these parents.
"""
function SHFX(
    v1::T,
    v2::T;
    rng::AbstractRNG = Random.default_rng(),
) where {T<:AbstractVector}
    l = length(v1)
    # shuffle and perform 1-point XO
    prm = randperm(rng, l)
    tmp1, tmp2 = SPX(view(v1, prm), view(v2, prm), rng = rng)
    # unshuffle offspring
    c1, c2 = similar(tmp1), similar(tmp2)
    for (i, j) in enumerate(prm)
        c1[j] = tmp1[i]
        c2[j] = tmp2[i]
    end
    return c1, c2
end
SHFX(v1::T, v2::T; rng::AbstractRNG = Random.default_rng()) where {T<:AbstractIndividual} =
    map(T, SHFX(variables(v1), variables(v2), rng=rng))

"""
    UX(v1, v2)

Uniform crossover between `v1` and `v2` individuals.
"""
function UX(v1::T, v2::T; rng::AbstractRNG = Random.default_rng()) where {T<:AbstractVector}
    l = length(v1)
    c1 = copy(v1)
    c2 = copy(v2)
    xch = rand(rng, Bool, l)
    for i = 1:l
        if xch[i]
            vswap!(c1, c2, i)
        end
    end
    return c1, c2
end
UX(v1::T, v2::T; rng::AbstractRNG = Random.default_rng()) where {T<:AbstractIndividual} =
    map(T, UX(variables(v1), variables(v2), rng=rng))

"""
    BINX(Cr::Real=0.5)

Returns a uniform (binomial) crossover function, see [Recombination Interface](@ref), function with the probability `Cr` [^2].

The crossover probability value must be in unit interval, ``Cr \\in [0,1]``.
"""
function BINX(Cr::Real = 0.5)
    function binxvr(
        v1::T,
        v2::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractVector}
        l = length(v1)
        c1 = copy(v1)
        c2 = copy(v2)
        j = rand(rng, 1:l)
        for i in (((1:l) .+ j .- 2) .% l) .+ 1
            if rand(rng) <= Cr
                vswap!(c1, c2, i)
            end
        end
        return c1, c2
    end
    function binxvr(
        v1::T,
        v2::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractIndividual}
        c1, c2 = binxvr(variables(v1), variables(v2), rng=rng)
        return T(c1), T(c2)
    end
    return binxvr
end

"""
    EXPX(Cr::Real=0.5)

Returns an exponential crossover function, see [Recombination Interface](@ref), function with the probability `Cr` [^2].

The crossover probability value must be in unit interval, ``Cr \\in [0,1]``.
"""
function EXPX(Cr::Real = 0.5)
    function expxvr(
        v1::T,
        v2::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractVector}
        l = length(v1)
        c1 = copy(v1)
        c2 = copy(v2)
        j = rand(rng, 1:l)
        switch = true
        for i in (((1:l) .+ j .- 2) .% l) .+ 1
            i == j && continue
            if switch && rand(rng) <= Cr
                c1[i] = v1[i]
            else
                switch = false
                c1[i] = v2[i]
            end
        end
        return c1, copy(c1)
    end
    function expxvr(
        v1::T,
        v2::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractIndividual}
        map(T, expxvr(variables(v1), variables(v2), rng=rng))
    end
    return expxvr
end

"""
    BSX(k::Int)

Binary Subset Crossover[^7]. Produces an offspring by first pooling the unique
items of the two parents, and then creating each offspring by sampling without
replacement at most `k` elements from the pool of items.
"""
function BSX(k::Int)
    function BSX(
        v1::T,
        v2::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractVector{Bool}}
        l = length(v1) # get number of available elements
        pooled = findall(v1 .| v2) # pool parents selections
        K = min(k, length(pooled)) # cannot sample more than the items in pool
        c1 = falses(l) # init child 1
        c2 = falses(l) # init child 2
        c1[shuffle(rng, pooled)[1:K]] .= true # fill child 1 with sample from pool w/o replacement
        c2[shuffle(rng, pooled)[1:K]] .= true # fill child 2 with sample from pool w/o replacement
        return c1, c2
    end
    function BSX(
        v1::T,
        v2::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractIndividual}
        map(T, BSX(variables(v1), variables(v2), rng=rng))
    end
    return BSX
end

# Real valued crossovers
# ----------------------

"""
    DC(v1, v2)

Returns a randomly assembled offspring and its inverse from the elements of parents `v1` and `v2`.
"""
function DC(v1::T, v2::T; rng::AbstractRNG = Random.default_rng()) where {T<:AbstractVector}
    l = length(v1)
    c1 = similar(v1)
    c2 = similar(v2)
    sltc = rand(rng, Bool, 2, l)
    for i = 1:l
        c1[i] = sltc[1, i] ? v1[i] : v2[i]
        c2[i] = sltc[2, i] ? v2[i] : v1[i]
    end
    return c1, c2
end
DC(v1::T, v2::T; rng::AbstractRNG = Random.default_rng()) where {T<:AbstractIndividual} =
    map(T, DC(variables(v1), variables(v2), rng=rng))

"""
    WAX(w::Vector{<:Real})(v1, v2)

Returns a weighted average recombination operation, see [Recombination Interface](@ref), which generate an offspring as weighted average of the parents `v1` and `v2` with the weights `w`.
"""
function WAX(w::Vector{<:Real})
    function wavexvr(
        v1::T,
        v2::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractVector}
        c1 = (v1 + v2) ./ w
        return c1, copy(c1)
    end
    function wavexvr(
        v1::T,
        v2::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractIndividual}
        map(T, wavexvr(variables(v1), variables(v2), rng=rng))
    end
    return wavexvr
end

"""
    AX(v1, v2)

Average crossover generates an offspring by taking average of the parents `v1` and `v2`.
"""
function AX(v1::T, v2::T; rng::AbstractRNG = Random.default_rng()) where {T<:AbstractVector}
    c1 = (v1 + v2) ./ 2
    return c1, copy(c1)
end
AX(v1::T, v2::T; rng::AbstractRNG = Random.default_rng()) where {T<:AbstractIndividual} =
    map(T, AX(variables(v1), variables(v2), rng=rng))

"""
    IC(d::Real=0.0)

Returns an extended intermediate recombination operation, see [Recombination Interface](@ref), which generates offspring `u` and `v` as

- ``u_i = x_i + \\alpha_i (y_i - x_i)``
- ``v_i = y_i + \\alpha_i (x_i - y_i)``

where ``\\alpha_i`` is chosen uniform randomly in the interval ``[-d;d+1]``.

"""
function IC(d::Real = 0.0)
    function intermxvr(
        v1::T,
        v2::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractVector}
        l = length(v1)
        α = (1.0 + 2d) * rand(rng, l) .- d
        c1 = v2 .+ α .* (v1 - v2)
        α = (1.0 + 2d) * rand(rng, l) .- d
        c2 = v1 .+ α .* (v2 - v1)
        return c1, c2
    end
    function intermxvr(
        v1::T,
        v2::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractIndividual}
        map(T, intermxvr(variables(v1), variables(v2), rng=rng))
    end
    return intermxvr
end

"""
    LC(d::Real=0.0)

Returns a extended line recombination operation, see [Recombination Interface](@ref), which generates offspring `u` and `v` as

- ``u_i = x_i + \\alpha (y_i - x_i)``
- ``v_i = y_i + \\alpha (x_i - y_i)``

where ``\\alpha`` is chosen uniform randomly in the interval ``[-d;d+1]``.

"""
function LC(d::Real = 0.0)
    function linexvr(
        v1::T,
        v2::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractVector}
        α1, α2 = (1.0 + 2d) * rand(rng, 2) .- d
        c1 = v2 .+ α2 * (v1 - v2)
        c2 = v1 .+ α1 * (v2 - v1)
        return c1, c2
    end
    function linexvr(
        v1::T,
        v2::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractIndividual}
        map(T, linexvr(variables(v1), variables(v2), rng=rng))
    end
    return linexvr
end

"""
    HX(x, y)

Heuristic crossover (HX) recombination operation[^3] generates offspring `u` and `v` as

- ``u = x + r (x - y)``
- ``v = y + r (y - x)``

where ``r`` is chosen uniform randomly in the interval ``[0;1)``.
"""
function HX(v1::T, v2::T; rng::AbstractRNG = Random.default_rng()) where {T<:AbstractVector}
    c1 = v1 .+ rand(rng) * (v1 .- v2)
    c2 = v2 .+ rand(rng) * (v2 .- v1)
    return c1, c2
end
HX(v1::T, v2::T; rng::AbstractRNG = Random.default_rng()) where {T<:AbstractIndividual} =
    map(T, HX(variables(v1), variables(v2), rng=rng))

"""
    LX(μ::Real = 0.0, b::Real = 0.2)

Returns a Laplace crossover (LX) recombination operation[^4], see [Recombination Interface](@ref).
"""
function LX(μ::Real = 0.0, b::Real = 0.2) # location μ, scale b > 0
    function lxxvr(
        v1::T,
        v2::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractVector}
        u = rand(rng)
        β = u > 0.5 ? μ + b * log(u) : μ - b * log(u)
        S = β * abs.(v1 - v2)
        c1 = v1 + S
        c2 = v2 + S
        return c1, c2
    end
    function lxxvr(
        v1::T,
        v2::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractIndividual}
        map(T, lxxvr(variables(v1), variables(v2), rng=rng))
    end
    return lxxvr
end

"""
    MILX(μ::Real = 0.0, b_real::Real = 0.15, b_int::Real = 0.35)

Returns a mixed integer Laplace crossover (MI-LX) recombination operation[^5], see [Recombination Interface](@ref).
"""
function MILX(μ::Real = 0.0, b_real::Real = 0.15, b_int::Real = 0.35) # location μ, scale b > 0
    function milxxvr(
        v1::T,
        v2::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractVector}
        @assert all([typeof(a) == typeof(b) for (a, b) in zip(v1, v2)]) "Types of variables in vectors do not match"
        l = length(v1)
        U, R = rand(rng, l), rand(rng, l)
        B = (isa(x, Integer) ? b_int : b_real for x in v1)
        βs = broadcast((u, r, b) -> r > 0.5 ? μ + b * log.(u) : μ - b * log.(u), U, R, B)
        S = βs .* abs.(v1 - v2)
        c1 = v1 + S
        c2 = v2 + S
        return c1, c2
    end
    function milxxvr(
        v1::T,
        v2::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractIndividual}
        map(T, milxxvr(variables(v1), variables(v2), rng=rng))
    end
    return milxxvr
end

"""
    SBX(pm::Real = 0.5, η::Integer = 2)

Returns a Simulated Binary Crossover (SBX) recombination operation, see [Recombination Interface](@ref),
with the mutation probability `pm` of the recombinant component, and is the crossover distribution index `η`[^6].
"""
function SBX(pm::Real = 0.5, η::Integer = 2)
    function sbxv(
        v1::T,
        v2::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractVector}
        n = length(v1)
        u = rand(rng, n)
        mask = u .<= 0.5
        mask_neg = u .> 0.5
        β = similar(v1)
        β[mask] .= (2 * u[mask]) .^ (1 / (η + 1))
        β[mask_neg] .= (2 * (1 .- u[mask_neg])) .^ (-1 / (η + 1))
        μ = (v1 + v2) ./ 2
        diff = v1 - v2
        c = β .* diff ./ 2
        mask_set = n == 1 ? [1] : rand(rng, n) .<= pm
        c1 = copy(μ) # c2 = μ - c
        c1[mask_set] .-= c[mask_set]
        c2 = copy(μ) # c2 = μ + c
        c2[mask_set] .+= c[mask_set]
        return c1, c2
    end
    function sbxv(
        v1::T,
        v2::T;
        rng::AbstractRNG = Random.default_rng(),
    ) where {T<:AbstractIndividual}
        map(T, sbxv(variables(v1), variables(v2), rng=rng))
    end
    return sbxv
end

# Combinatorial crossovers
# ------------------------

"""
    PMX(v1, v2)

Partially mapped crossover which maps ordering and values information from
the parents `v1` and `v2` to  the offspring. A portion of one parent is mapped onto
a portion of the other parent string and the remaining information is exchanged.
"""
function PMX(
    v1::T,
    v2::T;
    rng::AbstractRNG = Random.default_rng(),
) where {T<:AbstractVector}
    s = length(v1)
    from, to = rand(rng, 1:s, 2)
    from, to = from > to ? (to, from) : (from, to)
    c1 = similar(v1)
    c2 = similar(v2)

    # Swap
    c1[from:to] = v2[from:to]
    c2[from:to] = v1[from:to]

    # Fill in from parents
    for i in vcat(1:from-1, to+1:s)
        # Check conflicting offspring
        in1 = inmap(v1[i], c1, from, to)
        if in1 == 0
            c1[i] = v1[i]
        else
            tmpin = in1
            while tmpin > 0
                tmpin = inmap(c2[in1], c1, from, to)
                in1 = tmpin > 0 ? tmpin : in1
            end
            c1[i] = v1[in1]
        end

        in2 = inmap(v2[i], c2, from, to)
        if in2 == 0
            c2[i] = v2[i]
        else
            tmpin = in2
            while tmpin > 0
                tmpin = inmap(c1[in2], c2, from, to)
                in2 = tmpin > 0 ? tmpin : in2
            end
            c2[i] = v2[in2]
        end
    end
    return c1, c2
end
PMX(v1::T, v2::T; rng::AbstractRNG = Random.default_rng()) where {T<:AbstractIndividual} =
    map(T, PMX(variables(v1), variables(v2), rng=rng))

"""
    OX1(v1, v2)

Order crossover constructs an offspring by choosing a substring of one parent
and preserving the relative order of the elements of the other parent.
"""
function OX1(
    v1::T,
    v2::T;
    rng::AbstractRNG = Random.default_rng(),
) where {T<:AbstractVector}
    s = length(v1)
    from, to = rand(rng, 1:s, 2)
    from, to = from > to ? (to, from) : (from, to)
    c1 = zero(v1)
    c2 = zero(v2)
    # Swap
    c1[from:to] = v2[from:to]
    c2[from:to] = v1[from:to]
    # Fill in from parents
    k = to + 1 > s ? 1 : to + 1 #child1 index
    j = to + 1 > s ? 1 : to + 1 #child2 index
    for i in vcat(to+1:s, 1:from-1)
        while in(v1[k], c1)
            k = k + 1 > s ? 1 : k + 1
        end
        c1[i] = v1[k]
        while in(v2[j], c2)
            j = j + 1 > s ? 1 : j + 1
        end
        c2[i] = v2[j]
    end
    return c1, c2
end
OX1(v1::T, v2::T; rng::AbstractRNG = Random.default_rng()) where {T<:AbstractIndividual} =
    map(T, OX1(variables(v1), variables(v2), rng=rng))

"""
    CX(v1, v2)

Cycle crossover creates an offspring from the parents `v1` and `v2` such that
every position is occupied by a corresponding element from one of the parents.
"""
function CX(v1::T, v2::T; kwargs...) where {T<:AbstractVector}
    s = length(v1)
    c1 = zero(v1)
    c2 = zero(v2)
    c1i = fill(true, length(v1))
    c2i = fill(true, length(v2))
    Z = zero(eltype(T))

    f1 = true #switch
    k = 1
    while k > 0
        idx = k
        if f1
            #cycle from v1
            while c1i[idx]
                c1[idx] = v1[idx]
                c2[idx] = v2[idx]
                c1i[idx] = false
                c2i[idx] = false
                idx = inmap(v2[idx], v1, 1, s)
            end
        else
            #cycle from v2
            while c2i[idx]
                c1[idx] = v2[idx]
                c2[idx] = v1[idx]
                c1i[idx] = false
                c2i[idx] = false
                idx = inmap(v1[idx], v2, 1, s)
            end
        end
        f1 ⊻= true
        k = inmap(true, c2i, 1, s)
    end
    return c1, c2
end
CX(v1::T, v2::T; kwargs...) where {T<:AbstractIndividual} =
    map(T, CX(variables(v1), variables(v2); kwargs...))

"""
    OX2(v1, v2)

Order-based crossover selects at random several positions in the parent `v1`, and
the order of the elements in the selected positions of the parent `v1` is imposed on
the parent `v2`.
"""
function OX2(
    v1::T,
    v2::T;
    rng::AbstractRNG = Random.default_rng(),
) where {T<:AbstractVector}
    s = length(v1)
    c1 = copy(v1)
    c2 = copy(v2)
    Z = zero(eltype(T))

    for i = 1:s
        if rand(rng, Bool)
            idx1 = inmap(v2[i], v1, 1, s)
            idx2 = inmap(v1[i], v2, 1, s)
            c1[idx1] = Z
            c2[idx2] = Z
        end
    end

    for i = 1:s
        if !in(v2[i], c1)
            tmpin = inmap(Z, c1, 1, s)
            c1[tmpin] = v2[i]
        end
        if !in(v1[i], c2)
            tmpin = inmap(Z, c2, 1, s)
            c2[tmpin] = v1[i]
        end
    end
    return c1, c2
end
OX2(v1::T, v2::T; rng::AbstractRNG = Random.default_rng()) where {T<:AbstractIndividual} =
    map(T, OX2(variables(v1), variables(v2), rng=rng))

"""
    POS(v1, v2)

Position-based crossover is a modification of the [`OX1`](@ref) operator.
It selects a random set of positions in the parents `v1` and `v2`, then imposes
the position of the selected elements of one parent on the corresponding elements
of the other parent.
"""
function POS(
    v1::T,
    v2::T;
    rng::AbstractRNG = Random.default_rng(),
) where {T<:AbstractVector}
    s = length(v1)
    c1 = zero(v1)
    c2 = zero(v2)
    Z = zero(eltype(T))

    for i = 1:s
        if rand(rng, Bool)
            c1[i] = v2[i]
            c2[i] = v1[i]
        end
    end

    for i = 1:s
        if !in(v1[i], c1)
            tmpin = inmap(Z, c1, 1, s)
            c1[tmpin] = v1[i]
        end
        if !in(v2[i], c2)
            tmpin = inmap(Z, c2, 1, s)
            c2[tmpin] = v2[i]
        end
    end
    return c1, c2
end
POX(v1::T, v2::T; rng::AbstractRNG = Random.default_rng()) where {T<:AbstractIndividual} =
    map(T, POX(variables(v1), variables(v2), rng=rng))

"""
    SSX(v1, v2)

Subset crossover operator creates new offspring by pooling the unique indices of
the two parent vectors `v1` and `v2`, and then sampling a set of unique indices
from this pool, uniformly at random.
"""
function SSX(
    v1::T,
    v2::T;
    rng::AbstractRNG = Random.default_rng(),
) where {T<:AbstractVector}
    l = length(v1)
    pool = unique([v1; v2])
    p = length(pool)
    c1, c2 = if l <= p
        pool[randperm(p)[1:l]], pool[randperm(p)[1:l]]
    else
        rand(pool, l), rand(pool, l)
    end
    return c1, c2
end
SSX(v1::T, v2::T; rng::AbstractRNG = Random.default_rng()) where {T<:AbstractIndividual} =
    map(T, SSX(variables(v1), variables(v2), rng=rng))


text:
# Autogenerated wrapper script for Cairomm_jll for aarch64-linux-musl-cxx11
export libcariomm

using libsigcpp_jll
using Cairo_jll
JLLWrappers.@generate_wrapper_header("Cairomm")
JLLWrappers.@declare_library_product(libcariomm, "libcairomm-1.16.so.1")
function __init__()
    JLLWrappers.@generate_init_header(libsigcpp_jll, Cairo_jll)
    JLLWrappers.@init_library_product(
        libcariomm,
        "lib/libcairomm-1.16.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# The required / allowed input and output arguments per functional family
# :input refers to input argument and numbers to the derivative order
const ARGUMENTS = Dict(
    :lda => [[:zk], [:vrho], [:v2rho2], [:v3rho3], [:v4rho4], ],
    :gga => [
        [:zk],
        [:vrho, :vsigma],
        [:v2rho2, :v2rhosigma, :v2sigma2],
        [:v3rho3, :v3rho2sigma, :v3rhosigma2, :v3sigma3],
        [:v4rho4, :v4rho3sigma, :v4rho2sigma2, :v4rhosigma3, :v4sigma4],
    ],
    :mgga => [
        [:zk],
        [:vrho, :vsigma, :vlapl, :vtau],
        [:v2rho2, :v2rhosigma, :v2rholapl, :v2rhotau, :v2sigma2, :v2sigmalapl,
         :v2sigmatau, :v2lapl2, :v2lapltau, :v2tau2],
        [:v3rho3, :v3rho2sigma, :v3rho2lapl, :v3rho2tau, :v3rhosigma2, :v3rhosigmalapl,
         :v3rhosigmatau, :v3rholapl2, :v3rholapltau, :v3rhotau2, :v3sigma3, :v3sigma2lapl,
         :v3sigma2tau, :v3sigmalapl2, :v3sigmalapltau, :v3sigmatau2, :v3lapl3,
         :v3lapl2tau, :v3lapltau2, :v3tau3],
        [:v4rho4, :v4rho3sigma, :v4rho3lapl, :v4rho3tau, :v4rho2sigma2, :v4rho2sigmalapl,
         :v4rho2sigmatau, :v4rho2lapl2, :v4rho2lapltau, :v4rho2tau2, :v4rhosigma3,
         :v4rhosigma2lapl, :v4rhosigma2tau, :v4rhosigmalapl2, :v4rhosigmalapltau,
         :v4rhosigmatau2, :v4rholapl3, :v4rholapl2tau, :v4rholapltau2, :v4rhotau3,
         :v4sigma4, :v4sigma3lapl, :v4sigma3tau, :v4sigma2lapl2, :v4sigma2lapltau,
         :v4sigma2tau2, :v4sigmalapl3, :v4sigmalapl2tau, :v4sigmalapltau2, :v4sigmatau3,
         :v4lapl4, :v4lapl3tau, :v4lapl2tau2, :v4lapltau3, :v4tau4],
     ],
)
const INPUT = Dict(
    :lda => [:rho],
    :gga => [:rho, :sigma],
    :mgga => [:rho, :sigma, :lapl, :tau],
)

function derivative_flag(family, argument)
    flags = (:exc, :vxc, :fxc, :kxc, :lxc)  # flags for having 0th to 4th derivative
    for (i, arglist) in enumerate(ARGUMENTS[family])
        (argument in arglist) && return (i - 1, flags[i])
    end
    (-1, :unknown)
end

"""
Evaluate a functional obtaining the energy and / or certain derivatives of it.
What is returned depends on `derivatives`, which should be a range of derivative
orders, e.g. `0:1` (the default) for the energy and potential, `0:0` for just the
energy and so on.
The required input arguments depend on the functional type (`rho` for all functionals,
`sigma` for GGA and mGGA, `tau` and `lapl` for mGGA). Obtained data is returned
as a named tuple.
"""
function evaluate(func::Functional; derivatives=0:1, rho::AbstractArray, kwargs...)
    @assert all(0 .≤ derivatives .≤ 4)

    # Determine the gridshape (i.e. the shape of the grid points without the spin components)
    if ndims(rho) > 1
        if size(rho, 1) != func.spin_dimensions.rho
            error("First axis for multidimensional rho array should be equal " *
                  "to the number of spin components (== $(func.spin_dimensions.rho))")
        end
        gridshape = size(rho)[2:end]
    else
        if mod(length(rho), func.spin_dimensions.rho) != 0
            error("Length of linear rho array should be divisible by number of spin " *
                  "components in rho (== $(func.spin_dimensions.rho)).")
        end
        gridshape = (Int(length(rho) / func.spin_dimensions.rho), )
    end

    # Output arguments, where memory is already allocated
    outargs_allocated = Dict{Symbol, AbstractArray}()
    outargs = Dict{Symbol, AbstractArray}()
    for symbol in vcat(ARGUMENTS[func.family][1 .+ derivatives]...)
        if symbol in keys(kwargs)
            outargs_allocated[symbol] = kwargs[symbol]
        elseif symbol == :zk  # For zk keep just the grid shape
            outargs[symbol] = similar(rho, gridshape)
        else
            n_spin = getfield(func.spin_dimensions, symbol)
            outargs[symbol] = similar(rho, n_spin, gridshape...)
        end
    end

    evaluate!(func; rho=rho, kwargs..., outargs...)
    (; outargs..., outargs_allocated...)
end

"""
Evaluate a functional and store results in passed arrays. If for a particular
quantity no array is passed, it is not computed. Required input arguments
depend on the functional type (`rho` for all functionals, `sigma` for GGA and mGGA,
`tau` and `lapl` for mGGA).
"""
function evaluate!(func::Functional; rho::AbstractArray, kwargs...)
    mod(length(rho), func.spin_dimensions.rho) != 0 && error(
        "Length of rho array should be divisible by number of spin " *
        "components in rho (== $(func.spin_dimensions.rho))."
    )

    n_p = Int(length(rho) / func.spin_dimensions.rho)
    kwargs = Dict(kwargs)
    for argument in vcat(INPUT[func.family]..., ARGUMENTS[func.family]...)
        argument == :rho && continue
        haskey(kwargs, argument) || continue

        # Check dimensionality
        n_spin = getfield(func.spin_dimensions, argument)
        if length(kwargs[argument]) != n_spin * n_p
            throw(DimensionMismatch("$argument should have n_spin * n_p elements " *
                                    "(n_spin=$n_spin, n_p=$n_p"))
        end

        # Check derivative support in functional
        order, flag = derivative_flag(func.family, argument)
        if !(flag in func.flags) && flag != :unknown
            throw(ArgumentError("Functional $(func.identifier) does not support " *
                                "derivatives of order $order (like $argument)."))
        end

    end
    evaluate!(func, Val(func.family), rho; kwargs...)
end

# Array typedefs specifically to the Array type supported in Libxc
# Allows providing special versions of evaluate! for other types of arrays.
# Dispatch can be done on the rho argument, which is not a kwarg in this function.
const LibxcArray = Array{Float64}
const LibxcOptArray = Union{LibxcArray, Ptr{Nothing}}

function evaluate!(func::Functional, ::Val{:lda}, rho::LibxcArray;
                   zk::LibxcOptArray=C_NULL,
                   vrho::LibxcOptArray=C_NULL, v2rho2::LibxcOptArray=C_NULL,
                   v3rho3::LibxcOptArray=C_NULL, v4rho4::LibxcOptArray=C_NULL)
    np = Int(length(rho) / func.spin_dimensions.rho)
    xc_lda(func.pointer_, np, rho, zk, vrho, v2rho2, v3rho3, v4rho4)
end

function evaluate!(func::Functional, ::Val{:gga}, rho::LibxcArray; sigma::LibxcArray,
                   zk::LibxcOptArray=C_NULL,
                   vrho::LibxcOptArray=C_NULL, vsigma::LibxcOptArray=C_NULL,
                   v2rho2::LibxcOptArray=C_NULL, v2rhosigma::LibxcOptArray=C_NULL,
                   v2sigma2::LibxcOptArray=C_NULL,
                   v3rho3::LibxcOptArray=C_NULL, v3rho2sigma::LibxcOptArray=C_NULL,
                   v3rhosigma2::LibxcOptArray=C_NULL, v3sigma3::LibxcOptArray=C_NULL,
                   v4rho4::LibxcOptArray=C_NULL, v4rho3sigma::LibxcOptArray=C_NULL,
                   v4rho2sigma2::LibxcOptArray=C_NULL, v4rhosigma3::LibxcOptArray=C_NULL,
                   v4sigma4::LibxcOptArray=C_NULL)
    np = Int(length(rho) / func.spin_dimensions.rho)
    xc_gga(func.pointer_, np, rho, sigma, zk, vrho, vsigma, v2rho2, v2rhosigma, v2sigma2,
           v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3,
           v4rho4, v4rho3sigma, v4rho2sigma2, v4rhosigma3, v4sigma4)
end

function evaluate!(func::Functional, ::Val{:mgga}, rho::LibxcArray;
                   sigma::LibxcArray, lapl::LibxcArray, tau::LibxcArray,
                   zk::LibxcOptArray=C_NULL,
                   vrho::LibxcOptArray=C_NULL,
                   vsigma::LibxcOptArray=C_NULL,
                   vlapl::LibxcOptArray=C_NULL,
                   vtau::LibxcOptArray=C_NULL,
                   v2rho2::LibxcOptArray=C_NULL,
                   v2rhosigma::LibxcOptArray=C_NULL,
                   v2rholapl::LibxcOptArray=C_NULL,
                   v2rhotau::LibxcOptArray=C_NULL,
                   v2sigma2::LibxcOptArray=C_NULL,
                   v2sigmalapl::LibxcOptArray=C_NULL,
                   v2sigmatau::LibxcOptArray=C_NULL,
                   v2lapl2::LibxcOptArray=C_NULL,
                   v2lapltau::LibxcOptArray=C_NULL,
                   v2tau2::LibxcOptArray=C_NULL,
                   v3rho3::LibxcOptArray=C_NULL,
                   v3rho2sigma::LibxcOptArray=C_NULL,
                   v3rho2lapl::LibxcOptArray=C_NULL,
                   v3rho2tau::LibxcOptArray=C_NULL,
                   v3rhosigma2::LibxcOptArray=C_NULL,
                   v3rhosigmalapl::LibxcOptArray=C_NULL,
                   v3rhosigmatau::LibxcOptArray=C_NULL,
                   v3rholapl2::LibxcOptArray=C_NULL,
                   v3rholapltau::LibxcOptArray=C_NULL,
                   v3rhotau2::LibxcOptArray=C_NULL,
                   v3sigma3::LibxcOptArray=C_NULL,
                   v3sigma2lapl::LibxcOptArray=C_NULL,
                   v3sigma2tau::LibxcOptArray=C_NULL,
                   v3sigmalapl2::LibxcOptArray=C_NULL,
                   v3sigmalapltau::LibxcOptArray=C_NULL,
                   v3sigmatau2::LibxcOptArray=C_NULL,
                   v3lapl3::LibxcOptArray=C_NULL,
                   v3lapl2tau::LibxcOptArray=C_NULL,
                   v3lapltau2::LibxcOptArray=C_NULL,
                   v3tau3::LibxcOptArray=C_NULL,
                   v4rho4::LibxcOptArray=C_NULL,
                   v4rho3sigma::LibxcOptArray=C_NULL,
                   v4rho3lapl::LibxcOptArray=C_NULL,
                   v4rho3tau::LibxcOptArray=C_NULL,
                   v4rho2sigma2::LibxcOptArray=C_NULL,
                   v4rho2sigmalapl::LibxcOptArray=C_NULL,
                   v4rho2sigmatau::LibxcOptArray=C_NULL,
                   v4rho2lapl2::LibxcOptArray=C_NULL,
                   v4rho2lapltau::LibxcOptArray=C_NULL,
                   v4rho2tau2::LibxcOptArray=C_NULL,
                   v4rhosigma3::LibxcOptArray=C_NULL,
                   v4rhosigma2lapl::LibxcOptArray=C_NULL,
                   v4rhosigma2tau::LibxcOptArray=C_NULL,
                   v4rhosigmalapl2::LibxcOptArray=C_NULL,
                   v4rhosigmalapltau::LibxcOptArray=C_NULL,
                   v4rhosigmatau2::LibxcOptArray=C_NULL,
                   v4rholapl3::LibxcOptArray=C_NULL,
                   v4rholapl2tau::LibxcOptArray=C_NULL,
                   v4rholapltau2::LibxcOptArray=C_NULL,
                   v4rhotau3::LibxcOptArray=C_NULL,
                   v4sigma4::LibxcOptArray=C_NULL,
                   v4sigma3lapl::LibxcOptArray=C_NULL,
                   v4sigma3tau::LibxcOptArray=C_NULL,
                   v4sigma2lapl2::LibxcOptArray=C_NULL,
                   v4sigma2lapltau::LibxcOptArray=C_NULL,
                   v4sigma2tau2::LibxcOptArray=C_NULL,
                   v4sigmalapl3::LibxcOptArray=C_NULL,
                   v4sigmalapl2tau::LibxcOptArray=C_NULL,
                   v4sigmalapltau2::LibxcOptArray=C_NULL,
                   v4sigmatau3::LibxcOptArray=C_NULL,
                   v4lapl4::LibxcOptArray=C_NULL,
                   v4lapl3tau::LibxcOptArray=C_NULL,
                   v4lapl2tau2::LibxcOptArray=C_NULL,
                   v4lapltau3::LibxcOptArray=C_NULL,
                   v4tau4::LibxcOptArray=C_NULL)
    np = Int(length(rho) / func.spin_dimensions.rho)
    xc_mgga(func.pointer_, np, rho, sigma, lapl, tau, zk, vrho, vsigma, vlapl, vtau,
            v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, v2sigmalapl,
            v2sigmatau, v2lapl2, v2lapltau, v2tau2, v3rho3, v3rho2sigma,
            v3rho2lapl, v3rho2tau, v3rhosigma2, v3rhosigmalapl, v3rhosigmatau,
            v3rholapl2, v3rholapltau, v3rhotau2, v3sigma3, v3sigma2lapl,
            v3sigma2tau, v3sigmalapl2, v3sigmalapltau, v3sigmatau2, v3lapl3,
            v3lapl2tau, v3lapltau2, v3tau3, v4rho4, v4rho3sigma, v4rho3lapl,
            v4rho3tau, v4rho2sigma2, v4rho2sigmalapl, v4rho2sigmatau,
            v4rho2lapl2, v4rho2lapltau, v4rho2tau2, v4rhosigma3,
            v4rhosigma2lapl, v4rhosigma2tau, v4rhosigmalapl2,
            v4rhosigmalapltau, v4rhosigmatau2, v4rholapl3, v4rholapl2tau,
            v4rholapltau2, v4rhotau3, v4sigma4, v4sigma3lapl, v4sigma3tau,
            v4sigma2lapl2, v4sigma2lapltau, v4sigma2tau2, v4sigmalapl3,
            v4sigmalapl2tau, v4sigmalapltau2, v4sigmatau3, v4lapl4, v4lapl3tau,
            v4lapl2tau2, v4lapltau3, v4tau4)
end