text:
using Knet
using Sloth
import Base: eltype
using Random
using Statistics

function randinit(d...; scale=0.1)
    A = rand(d...)
    A = 2A .* scale .- scale
end

function clip(x)
    x = max.(0, x)
    x = min.(1, x)
end

struct Chain; layers; end
(c::Chain)(x) = (for l in c.layers; x = l(x); end; x)

mutable struct VanillaRNN
    i2h
    h2h
end

function VanillaRNN(xsize::Int, hsize::Int; atype=Sloth._atype, init=xavier)
    return VanillaRNN(Linear(input=xsize, output=hsize, atype=atype),
                      Linear(input=hsize, output=hsize, atype=atype))
end

(m::VanillaRNN)(x, ht) = relu.(m.i2h(x) + m.h2h(ht))

mutable struct Net
    σ
    num_glimpses
    glimpse_net
    controller
    location_net
    softmax_layer
    baseline_net
    locations
    history
end

function Net(patch_size, num_patches, glimpse_scale, num_channels, loc_hidden,
             glimpse_hidden, σ, hidden_size, num_classes, num_glimpses;
             rnnType=:relu, atype=Sloth._atype, init=randinit)
    usegpu = atype <: KnetArray
    etype = eltype(atype)
    return Net(
        σ,
        num_glimpses,
        GlimpseNet(loc_hidden, glimpse_hidden, patch_size,
                   num_patches, glimpse_scale, num_channels;
                   atype=atype, init=init),
        RNN(hidden_size, hidden_size;
            rnnType=rnnType, usegpu=usegpu, dataType=etype, winit=init),
        # VanillaRNN(hidden_size, hidden_size; atype=atype, init=init),
        LocationNet(hidden_size, 2, σ; atype=atype, init=init),
        Linear(input=hidden_size, output=num_classes, atype=atype, init=init),
        BaselineNet(input=hidden_size, output=1, atype=atype, init=init),
        [],
        zeros(6)) # 1-3=>loss,4=>iter,5=>correct,6=># of samples
end

# one step
function (ram::Net)(x, ltprev; deterministic=false)
    gt = ram.glimpse_net(x, ltprev)
    ram.controller(gt)
    ht = ram.controller.h
    H, B = size(ht); ht = reshape(ht, H, B)
    μ, lt = ram.location_net(ht)
    bt = clip(ram.baseline_net(value(ht)))
    F = eltype(gt)
    σ = F(ram.σ)
    σ² = σ .^ 2
    logπ = -(abs.(lt - μ) .^ 2) / 2σ² .- log(σ) .- log(√(F(2)π))
    logπ = sum(logπ, dims=1)

    if deterministic
        lt = μ
    end
    return lt, bt, logπ
end

# all timesteps
function (ram::Net)(x; deterministic=false, centered=false)
    B = size(x)[end]
    lt, ht = initstates(ram, B; centered=centered)
    ram.controller.h = 0
    ram.controller.c = 0
    logπs, baselines, locations = [], [], Any[lt]
    for t = 1:ram.num_glimpses
        lt, bt, logπ = ram(x, lt; deterministic=deterministic)
        push!(locations, lt)
        push!(baselines, bt)
        push!(logπs, logπ)
    end
    baseline, logπ = vcat(baselines...), vcat(logπs...)
    ht = ram.controller.h
    H, B = size(ht); ht = reshape(ht, H, B)
    scores = ram.softmax_layer(mat(ht))
    return scores, baseline, logπ, locations
end

# loss
function (ram::Net)(x, y::Union{Array{Int64}, Array{UInt8}};
                    deterministic=false, centered=false)
    atype = artype(ram)
    etype = eltype(ram)
    scores, baseline, logπ, locations = ram(
        x; deterministic=deterministic, centered=centered)
    grid = exploration(ram, locations)
    rate = reshape(mean(grid .> 0, dims=(1,2)), 1, :);
    ŷ = vec(map(i->i[1], argmax(Array(value(scores)), dims=1)))
    r = ŷ .== y; r = reshape(r, 1, :)
    # r = r .+ rate;
    R = convert(atype{etype}, r)
    # R = zeros(Float32, size(baseline)...); R[end,:] = r
    # R = convert(atype{etype}, R)
    R̂ = R .- value(baseline)
    loss_action = nll(scores, y)
    loss_baseline = sum(abs2, baseline .- R) / length(baseline)
    loss_reinforce = mean(sum(-logπ .* R̂, dims=1))
    return loss_action, loss_baseline, loss_reinforce, sum(r), length(r)
end
loss(ram::Net, x, ygold) = sum(ram(x,ygold)[1:3])
loss(ram::Net, d::Knet.Data) = mean(sum(ram(x,y)[1:3]) for (x,y) in d)

function validate(ram::Net, data; deterministic=false, centered=false)
    losses = zeros(3)
    ncorrect = ninstances = 0
    for (x,y) in data
        ret = ram(x,y; deterministic=deterministic, centered=centered)
        for i = 1:3; losses[i] += ret[i]; end
        ncorrect += ret[4]
        ninstances += ret[5]
    end
    losses = losses / length(data)
    losses = [sum(losses), losses...]
    return losses, ncorrect / ninstances
end

function initstates(ram::Net, B; centered=false)
    etype, atype = eltype(ram), artype(ram)
    hsize = get_hsize(ram)
    lsize = get_lsize(ram)
    l₀ = atype{etype}(2rand(lsize, B) .- 1)
    if centered
        l₀ = atype(zeros(etype, lsize, B))
    end
    h₀ = atype{etype}(zeros(hsize, B))
    return l₀, h₀
end

eltype(ram::Net) = eltype(ram.softmax_layer.w.value)
artype(ram::Net) = ifelse(
    typeof(ram.softmax_layer.w.value) <: KnetArray,
    KnetArray,
    Array)
get_hsize(ram::Net) = get_hsize(ram.controller)
get_hsize(m::RNN) = m.hiddenSize
get_hsize(m::VanillaRNN) = size(m.h2h.w)[1]
get_lsize(ram::Net) = size(ram.location_net.layer.w)[1]

mutable struct GlimpseNet
    fc1
    fc2
    fc3
    fc4
    retina
end

function GlimpseNet(
    loc_hidden, glimpse_hidden, patch_size, num_patches, scale, num_channels;
    atype=Sloth._atype, init=randinit)
    retina = Retina(patch_size, num_patches, scale)
    fc1 = Dense(
        input=patch_size*patch_size*num_patches*num_channels,
        output=glimpse_hidden, atype=atype, init=init)
    fc2 = Dense(input=2, output=loc_hidden, atype=atype, init=init)
    fc3 = Linear(input=glimpse_hidden, output=glimpse_hidden+loc_hidden,
                 atype=atype, init=init)
    fc4 = Linear(input=loc_hidden, output=glimpse_hidden+loc_hidden,
                 atype=atype, init=init)
    return GlimpseNet(fc1, fc2, fc3, fc4, retina)
end

function (m::GlimpseNet)(x, lt)
    phi = m.retina(x, lt)
    lt = mat(lt)
    atype = typeof(m.fc1.w.value) <: KnetArray ? KnetArray : Array
    etype = eltype(m.fc1.w.value)
    phi = atype{etype}(phi)
    yphi = m.fc1(phi)
    yloc = m.fc2(lt)
    y = m.fc3(yphi)
    l = m.fc4(yloc)
    return relu.(y .* l)
end

mutable struct Retina
    patch_size
    num_patches
    glimpse_scale
end

# no need for Retina constructor!

# foveate
function (r::Retina)(x, l)
    phi = []
    sz = r.patch_size
    atype = typeof(l) <: KnetArray ? KnetArray : Array
    etype = eltype(l)

    for k = 1:r.num_patches
        push!(phi, r(x, l, sz))
        sz = Int(r.glimpse_scale * sz)
    end

    for i = 1:length(phi)
        k = div(size(phi[i])[1], r.patch_size)
        phi[i] = pool(phi[i]; window=k, mode=1)
    end

    phi2d = map(mat, phi)
    return atype{etype}(mat(vcat(phi2d...)))
end

# extract patch
# size(x) = W,H,C,B
# size(l) = 2,B
function (r::Retina)(x, l, k)
    W, H, C, B = size(x)
    coords = denormalize(H, l)

    patch_x = coords[1,:] .- div(k,2) .+ 1.
    patch_y = coords[2,:] .- div(k,2) .+ 1.

    patches = []
    atype = typeof(x) <: KnetArray ? KnetArray : Array
    etype = eltype(x)
    x2d = mat(x)
    for i = 1:B
        img = reshape(x2d[:, i:i], W, H, C, 1)
        img = Array(img)
        T = size(img)[1]

        sz = r.patch_size
        from_x, from_y = Int(round(patch_x[i])), Int(round(patch_y[i]))
        to_x, to_y = from_x+sz-1, from_y+sz-1

        if exceeds(from_x, to_x, from_y, to_y, T)
            pd = div(sz,2)+1
            padded = zeros(2pd+W, 2pd+H, C, 1)
            from_x += pd+1; to_x += pd+1; from_y += pd+1; to_y += pd+1
            padded[pd+1:pd+W, pd+1:pd+H, :, :] = img
            img = padded
        end
        push!(patches, img[from_x:to_x,from_y:to_y,:,:])
    end

    glimpses = atype{etype}(cat(patches...; dims=4))
    return glimpses
end

function denormalize(T, coords)
    return convert(Array{Int}, floor.(0.5T * (coords .+ 1.0)))
end

function exceeds(from_x, to_x, from_y, to_y, T)
    return from_x < 1 || from_y < 1 || to_x > T || to_y > T
end

mutable struct LocationNet
    σ
    layer
end

function LocationNet(input_size::Int, output_size::Int, σ;
                     atype=Sloth._atype, init=randinit)
    layer = Dense(input=input_size, output=output_size, f=clip,
                  atype=atype, init=init)
    return LocationNet(σ, layer)
end

function (m::LocationNet)(ht)
    μ = clip(m.layer(value(ht)))
    F = eltype(μ)
    noise = F(m.σ) .* randn!(similar(μ))
    lt = μ + noise
    return μ, clip(lt)
end

BaselineNet = Linear

max_stars_repo_name: ilkerkesen/RAM

max_issues_repo_name: ilkerkesen/RAM

max_forks_repo_name: ilkerkesen/RAM


text:
struct DummyRoot end

BB.get_branching_indices(::DummyRoot) = 1:1

function BB.evaluate_node!(tree::BnBTree{BB.DefaultNode, DummyRoot}, node::BB.DefaultNode)
    lb = -0.9
    ub = sum(BB.get_relaxed_values(tree, node))
    if node.id % 2 == 0
        lb = 1.0
        ub = 1.0
    end
    return lb, ub
end

function BB.get_relaxed_values(::BnBTree{BB.DefaultNode, DummyRoot}, node::BB.DefaultNode)
    return [-node.id/10]
end

function BB.get_branching_nodes_info(tree::BnBTree{BB.DefaultNode, DummyRoot}, node::BB.DefaultNode, vidx::Int)
    node_info = NamedTuple[]
    push!(node_info, NamedTuple())
    push!(node_info, NamedTuple())
    return node_info
end

function dummy_callback(tree, node; node_infeasible=false, worse_than_incumbent=false)
    if node.id % 2 == 0
        @test worse_than_incumbent
    else
        @test !worse_than_incumbent
    end
    @test !node_infeasible
end

@testset "Dummy callback" begin
    bnb_model = BB.initialize(;
        root = DummyRoot(),
        sense = :Min
    )

    BB.set_root!(bnb_model, NamedTuple())

    BB.optimize!(bnb_model; callback=dummy_callback)

    @test -0.9 ≈ BB.get_solution(bnb_model)[1]
    @test BB.get_objective_value(bnb_model) ≈ -0.9

    @test BB.get_num_solutions(bnb_model) > 1
    @test BB.get_objective_value(bnb_model; result=2) >  BB.get_objective_value(bnb_model; result=1)
end

ext: jl


text:
import Base: ==

"""
`BrewPkg`

A simple type to give us some nice ways of representing our packages to the user

It contains important information such as the `name` of the package, the `tap` it
came from, the `version` of the package and whether it was `translated` or not
"""
struct BrewPkg
    # The name of this particular Brew package
    name::String

    # The tap this brew package comes from ("Homebrew/core" in general)
    tap::String

    # The version of this brew package
    version::String
end

function ==(x::BrewPkg, y::BrewPkg)
    return x.name == y.name && x.tap == y.tap && x.version == y.version
end

"""
`fullname(pkg::BrewPkg)`

Return the fully-qualified name for a package, dropping "Homebrew/core"
"""
function fullname(pkg::BrewPkg)
    if pkg.tap == "Homebrew/core"
        return pkg.name
    end
    return joinpath(pkg.tap,pkg.name)
end

"""
`show(io::IO, b::BrewPkg)`

Writes a `BrewPkg` to `io`, showing tap, name and version number
"""
function show(io::IO, b::BrewPkg)
    write(io, "$(fullname(b)): $(b.version)")
end

"""
`StringOrPkg`

A convenience type accepting either an `AbstractString` or a `BrewPkg`
"""
const StringOrPkg = Union{AbstractString, BrewPkg}


text:
using LibGit2, LibGit2_jll, SHA, Scratch

"""
    cached_git_clone(url::String); desired_commit = nothing, verbose = false)

Return the path to a local git clone of the given `url`.  If `desired_commit` is given,
then a cached git repository will not be updated if the commit already exists locally.
"""
function cached_git_clone(url::String;
                          desired_commit::Union{Nothing, String} = nothing,
                          clones_dir::String = @get_scratch!("git_clones"),
                          verbose::Bool = false)
    repo_path = joinpath(clones_dir, string(basename(url), "-", bytes2hex(sha256(url))))
    if isdir(repo_path)
        if verbose
            @info("Using cached git repository", url, repo_path)
        end
        # If we didn't just mercilessly obliterate the cached git repo, use it!
        LibGit2.with(LibGit2.GitRepo(repo_path)) do repo
            # In some cases, we know the hash we're looking for, so only fetch() if
            # this git repository doesn't contain the hash we're seeking.
            # this is not only faster, it avoids race conditions when we have
            # multiple builders on the same machine all fetching at once.
            if desired_commit === nothing || !LibGit2.iscommit(desired_commit, repo)
                LibGit2.fetch(repo)
            end
        end
    else
        if verbose
            @info("Cloning git repository", url, repo_path)
        end
        # If there is no repo_path yet, clone it down into a bare repository
        LibGit2.clone(url, repo_path; isbare=true)
    end
    return repo_path
end

"""
    get_tree_hash(tree::LibGit2.GitTree)

Given a `GitTree`, get the `GitHash` that identifies it.
"""
function get_tree_hash(tree::LibGit2.GitTree)
    oid_ptr = Ref(LibGit2.GitHash())
    oid_ptr = ccall((:git_tree_id, libgit2), Ptr{LibGit2.GitHash}, (Ptr{Cvoid},), tree.ptr)
    oid_ptr == C_NULL && throw(ArgumentError("bad tree ID: $tree"))
    return unsafe_load(oid_ptr)
end

"""
    get_commit_sha(url::String, tree_hash::Base.SHA1; verbose::Bool=false)

Find the latest git commit corresponding to the given git tree SHA1 for the remote
repository with the given `url`.  The repository is cached locally for quicker future
access.  If `verbose` is `true`, print to screen some debugging information.
The return value is the commit SHA as a `String`, if the corresponding revision is found,
`nothing` otherwise.
"""
function get_commit_sha(url::String, tree_hash::Base.SHA1; verbose::Bool=false)
    git_commit_sha = nothing
    dir = cached_git_clone(url; verbose)

    LibGit2.with(LibGit2.GitRepo(dir)) do repo
        LibGit2.with(LibGit2.GitRevWalker(repo)) do walker
            # The repo is cached, so locally it may be checking out an outdated commit.
            # Start the search from HEAD of the tracking upstream repo.
            try
                LibGit2.push!(walker, LibGit2.GitHash(LibGit2.peel(LibGit2.GitCommit, LibGit2.upstream(LibGit2.head(repo)))))
            catch
                @warn("Could not walk from origin branch!")
                LibGit2.push_head!(walker)
            end
            # For each commit in the git repo, check to see if its treehash
            # matches the one we're looking for.
            for oid in walker
                tree = LibGit2.peel(LibGit2.GitTree, LibGit2.GitCommit(repo, oid))
                if all(get_tree_hash(tree).val .== tree_hash.bytes)
                    git_commit_sha = LibGit2.string(oid)
                    break
                end
            end
        end
    end
    return git_commit_sha
end


text:
using Spglib
using Test

@testset "Spglib.jl" begin
    # Write your own tests here.
    include("symmetry.jl")
    include("standardize.jl")
    include("reciprocal.jl")
    include("reduce.jl")
    include("dataset.jl")
end


text:
# conversions
mmHgToPa = 133.32;
MPaToPa = 1e6;
cmTom = 1e-2;
mmTom = 1e-3;
cm3Tom3 = 1e-6;
minTos = 60;
gTokg = 1e-3;

# non-dimensional scalings
us = 200*mmTom;

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
#InspectDR: Draw legends with Cairo layer.
#-------------------------------------------------------------------------------

#==
===============================================================================#

#NOTE: linelength & textoffset are absolute values
function legend_renderitem(ctx::CairoContext, wfrm::DWaveform, x::Float64, y::Float64,
	linelength::Float64, textoffset::Float64)
	linestart = x+textoffset
	lineend = linestart+linelength

Cairo.save(ctx)
	setlinestyle(ctx, LineStyle(wfrm.line))
	drawline(ctx, Point2D(linestart, y), Point2D(lineend, y))
	drawglyph_safe(ctx, wfrm, Point2D(linestart+linelength/2, y))
Cairo.restore(ctx)
	x = lineend+textoffset #Compute new x
	render(ctx, wfrm.id, Point2D(x, y), align=ALIGN_VCENTER|ALIGN_LEFT)
end

function legend_render(canvas::PCanvas2D, plot::Plot2D, istrip::Int)
	ctx = canvas.ctx #WANTCONST
	lyt = plot.layout.values #WANTCONST

	xleft = canvas.bb.xmax - lyt.halloc_legend
	bb = BoundingBox(xleft, canvas.bb.xmax, canvas.graphbb.ymin, canvas.graphbb.ymin)
	Cairo.save(canvas.ctx)
	drawrectangle(ctx, bb, lyt.frame_legend)
	Cairo.restore(canvas.ctx)

	Cairo.save(canvas.ctx)
	setfont(ctx, lyt.font_legend)
	(w, h) = text_dims(Cairo.text_extents(ctx, "M"))
	ypitch = h*(1+lyt.valloc_legenditemsp)

	#Compute absolute values:
	textoffset = lyt.hoffset_legendtext*w

	y = canvas.graphbb.ymin + h/2
	for d in plot.display_data
		if d.strip != istrip; continue; end
		if "" == d.id; continue; end
		legend_renderitem(ctx, d, xleft, y, lyt.halloc_legendlineseg, textoffset)
		y += ypitch
	end
	Cairo.restore(canvas.ctx)
end

#Last line


text:
# Kernels that handle broadcasting with different sized arrays.
# Full broadcasting index computation is complicated.  I am going to handle the
# case where the input arrays are either the same size as the result or vectors
# (have a single dim > 1).  If the result array has size (a1,a2,...) and input
# array b has size (1,...,bn,...,1), we need to find the linear index in b given
# the linear index in a.  If we are at position (i1,i2,...) the answer is just
# i_n=mod(div(i,stride(a,n)),size(a,n)) with 0 indexing.  So we can just pass in
# stride(a,n) and size(a,n) as an argument for each input.

fp = open("cuda12.cu","w")
#using Knet: broadcast_ops

function cuda12src(f, j=f, ex="$f(xi,yi)"; BLK=256, THR=256)
  sprint() do s
    for (T,F) in [("float","$(f)_32"),("double","$(f)_64")]
        print(s,
"""
__global__ void _$(F)_12(int n, $T *x, int sx, int nx, $T *y, int sy, int ny, $T *z) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    $T xi = (nx==n ? x[i] : sx==1 ? x[i%nx] : nx==1 ? x[0] : x[(i/sx)%nx]);
    $T yi = (ny==n ? y[i] : sy==1 ? y[i%ny] : ny==1 ? y[0] : y[(i/sy)%ny]);
    z[i] = $ex;
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  $DLLEXPORT void $(F)_12(int n, $T *x, int sx, int nx, $T *y, int sy, int ny, $T *z) {
    _$(F)_12<<<$BLK,$THR>>>(n,x,sx,nx,y,sy,ny,z);
  }
}
""")
    end
  end
end

for a in broadcast_ops
    if !isa(a,Tuple); a=(a,); end
    print(fp,cuda12src(a...))
end

close(fp)


text:
function coroutine_pagecontrol(ws)
    @info ws
    while isopen(ws)
        data, success = readguarded(ws)
        if success
            @show String(data)
        end
    end
    @info "Closing ", ws
end

max_stars_repo_name: hustf/EquationEditor

max_issues_repo_name: hustf/EquationEditor

max_forks_repo_name: hustf/EquationEditor


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

mutable struct IoK8sApiAutoscalingV2beta2MetricTarget <: SwaggerModel
    averageUtilization::Any # spec type: Union{ Nothing, Int32 } # spec name: averageUtilization
    averageValue::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApiResourceQuantity } # spec name: averageValue
    type::Any # spec type: Union{ Nothing, String } # spec name: type
    value::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApiResourceQuantity } # spec name: value

    function IoK8sApiAutoscalingV2beta2MetricTarget(;averageUtilization=nothing, averageValue=nothing, type=nothing, value=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV2beta2MetricTarget, Symbol("averageUtilization"), averageUtilization)
        setfield!(o, Symbol("averageUtilization"), averageUtilization)
        validate_property(IoK8sApiAutoscalingV2beta2MetricTarget, Symbol("averageValue"), averageValue)
        setfield!(o, Symbol("averageValue"), averageValue)
        validate_property(IoK8sApiAutoscalingV2beta2MetricTarget, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        validate_property(IoK8sApiAutoscalingV2beta2MetricTarget, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        o
    end
end # type IoK8sApiAutoscalingV2beta2MetricTarget

const _property_map_IoK8sApiAutoscalingV2beta2MetricTarget = Dict{Symbol,Symbol}(Symbol("averageUtilization")=>Symbol("averageUtilization"), Symbol("averageValue")=>Symbol("averageValue"), Symbol("type")=>Symbol("type"), Symbol("value")=>Symbol("value"))
const _property_types_IoK8sApiAutoscalingV2beta2MetricTarget = Dict{Symbol,String}(Symbol("averageUtilization")=>"Int32", Symbol("averageValue")=>"IoK8sApimachineryPkgApiResourceQuantity", Symbol("type")=>"String", Symbol("value")=>"IoK8sApimachineryPkgApiResourceQuantity")
Base.propertynames(::Type{ IoK8sApiAutoscalingV2beta2MetricTarget }) = collect(keys(_property_map_IoK8sApiAutoscalingV2beta2MetricTarget))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV2beta2MetricTarget }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAutoscalingV2beta2MetricTarget[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV2beta2MetricTarget }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV2beta2MetricTarget[property_name]

function check_required(o::IoK8sApiAutoscalingV2beta2MetricTarget)
    (getproperty(o, Symbol("type")) === nothing) && (return false)
    true
end

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


text:
"""
    ClimaCore.DataLayouts

Notation:
- `i,j` are horizontal node indices within an element
- `k` is the vertical node index within an element
- `f` is the field index
- `v` is the vertical element index in a stack
- `h` is the element stack index

Data layout is specified by the order in which they appear, e.g. `IJKFVH`
indexes the underlying array as `[i,j,k,f,v,h]`

"""
module DataLayouts

import Base: Base, @propagate_inbounds
import StaticArrays: SOneTo, MArray

import ..slab, ..slab_args, ..column, ..column_args, ..level
export slab, column, IJFH, IJF, IFH, IF, VF, VIJFH, VIFH

include("struct.jl")

abstract type AbstractData{S} end

Base.size(data::AbstractData, i::Integer) = size(data)[i]

function Base.show(io::IO, data::AbstractData)
    indent_width = 2
    (rows, cols) = displaysize(io)
    println(io, summary(data))
    print(io, " "^indent_width)
    print(
        IOContext(
            io,
            :compact => true,
            :limit => true,
            :displaysize => (rows, cols - indent_width),
        ),
        vec(parent(data)),
    )
    return io
end

"""
    DataColumn{S}

Abstract type for data storage for a column. Objects `data` should define a
`data[k,v]`, returning a value of type `S`.
"""
abstract type DataColumn{S} <: AbstractData{S} end

"""
    DataSlab1D{S,Ni}

Abstract type for data storage for a slab of `Ni` values of type `S`.
Objects `data` should define a `data[i]`, returning a value of type `S`.
"""
abstract type DataSlab1D{S, Nij} <: AbstractData{S} end

"""
    DataSlab2D{S,Nij}

Abstract type for data storage for a slab of `Nij × Nij` values of type `S`.
Objects `data` should define a `data[i,j]`, returning a value of type `S`.
"""
abstract type DataSlab2D{S, Nij} <: AbstractData{S} end

"""
    Data1D{S,Ni}

Abstract type for data storage for a 1D field made up of `Ni` values of type `S`.

Objects `data` should define `slab(data, h)` to return a `DataSlab2D{S,Nij}` object.
"""
abstract type Data1D{S, Ni} <: AbstractData{S} end

"""
    Data2D{S,Nij}

Abstract type for data storage for a 2D field made up of `Nij × Nij` values of type `S`.

Objects `data` should define `slab(data, h)` to return a `DataSlab2D{S,Nij}` object.
"""
abstract type Data2D{S, Nij} <: AbstractData{S} end

"""
    Data1DX{S,Ni}

Abstract type for data storage for a 1D field with extruded columns.
The horizontal is made up of `Ni` values of type `S`.

Objects `data` should define `slab(data, v, h)` to return a
`DataSlab1D{S,Ni}` object, and a `column(data, i, h)` to return a `DataColumn`.
"""
abstract type Data1DX{S, Ni} <: AbstractData{S} end

"""
    Data2DX{S,Nij}

Abstract type for data storage for a 2D field with extruded columns.
The horizontal is made  is made up of `Nij × Nij` values of type `S`.

Objects `data` should define `slab(data, v, h)` to return a
`DataSlab2D{S,Nij}` object, and a `column(data, i, j, h)` to return a `DataColumn`.
"""
abstract type Data2DX{S, Nij} <: AbstractData{S} end

"""
    Data3D{S,Nij,Nk}

Abstract type for data storage for a 3D field made up of `Nij × Nij × Nk` values of type `S`.
"""
abstract type Data3D{S, Nij, Nk} <: AbstractData{S} end

# Generic AbstractData methods

Base.eltype(::AbstractData{S}) where {S} = S
@inline function Base.propertynames(::AbstractData{S}) where {S}
    filter(name -> sizeof(fieldtype(S, name)) > 0, fieldnames(S))
end

Base.parent(data::AbstractData) = getfield(data, :array)

Base.similar(data::AbstractData{S}) where {S} = similar(data, S)

function Base.copyto!(dest::D, src::D) where {D <: AbstractData}
    copyto!(parent(dest), parent(src))
    return dest
end

function nfields(data::AbstractData{S}) where {S}
    typesize(eltype(parent(data)), S)
end

@generated function _getproperty(
    data::AbstractData{S},
    ::Val{Name},
) where {S, Name}
    errorstring = "Invalid field name $(Name) for type $(S)"
    i = findfirst(isequal(Name), fieldnames(S))
    if i === nothing
        return :(error($errorstring))
    end
    static_idx = Val{i}()
    return :(Base.@_inline_meta; DataLayouts._property_view(data, $static_idx))
end

@inline function Base.getproperty(data::AbstractData{S}, name::Symbol) where {S}
    _getproperty(data, Val{name}())
end

#= Generic function impl fallback
@noinline function error_invalid_fieldname(@nospecialize(S::Type), name::Symbol)
    error("Invalid field name $(name) for type $(S)")
end

@inline function Base.getproperty(data::AbstractData{S}, name::Symbol) where {S}
   i = findfirst(isequal(name), fieldnames(S))
   i === nothing && error_invalid_fieldname(S, name)
   getproperty(data, i)
end
=#

# ==================
# Data3D DataLayout
# ==================

struct IJKFVH{S, Nij, Nk, A} <: Data3D{S, Nij, Nk}
    array::A
end

function IJKFVH{S, Nij, Nk}(array::AbstractArray{T, 6}) where {S, Nij, Nk, T}
    @assert size(array, 1) == Nij
    @assert size(array, 2) == Nij
    @assert size(array, 3) == Nk
    check_basetype(T, S)
    IJKFVH{S, Nij, Nk, typeof(array)}(array)
end

function replace_basetype(
    data::IJKFVH{S, Nij, Nk},
    ::Type{T},
) where {S, Nij, Nk, T}
    array = parent(data)
    S′ = replace_basetype(eltype(array), T, S)
    return IJKFVH{S′, Nij, Nk}(similar(array, T))
end

@generated function _property_view(
    data::IJKFVH{S, Nij, Nk, A},
    ::Val{Idx},
) where {S, Nij, Nk, A, Idx}
    SS = fieldtype(S, Idx)
    T = eltype(A)
    offset = fieldtypeoffset(T, S, Idx)
    nbytes = typesize(T, SS)
    field_byterange = (offset + 1):(offset + nbytes)
    return :(IJKFVH{$SS, $Nij, $Nk}(
        @inbounds view(parent(data), :, :, :, $field_byterange, :, :)
    ))
end

@inline function Base.getproperty(
    data::IJKFVH{S, Nij, Nk},
    i::Integer,
) where {S, Nij, Nk}
    array = parent(data)
    T = eltype(array)
    SS = fieldtype(S, i)
    offset = fieldtypeoffset(T, S, i)
    nbytes = typesize(T, SS)
    dataview =
        @inbounds view(array, :, :, :, (offset + 1):(offset + nbytes), :, :)
    IJKFVH{SS, Nij, Nk}(dataview)
end

function Base.size(data::IJKFVH{S, Nij, Nk}) where {S, Nij, Nk}
    Nv = size(parent(data), 5)
    Nh = size(parent(data), 6)
    return (Nij, Nij, Nk, Nv, Nh)
end

# ==================
# Data2D DataLayout
# ==================

"""
    IJFH{S, Nij, A} <: Data2D{S, Nij}

Backing `DataLayout` for 2D spectral element slabs.

Element nodal point (I,J) data is contiguous for each datatype `S` struct field (F),
for each 2D mesh element slab (H).
"""
struct IJFH{S, Nij, A} <: Data2D{S, Nij}
    array::A
end

function IJFH{S, Nij}(array::AbstractArray{T, 4}) where {S, Nij, T}
    @assert size(array, 1) == Nij
    @assert size(array, 2) == Nij
    check_basetype(T, S)
    IJFH{S, Nij, typeof(array)}(array)
end

rebuild(data::IJFH{S, Nij}, array) where {S, Nij} = IJFH{S, Nij}(array)

Base.copy(data::IJFH{S, Nij}) where {S, Nij} = IJFH{S, Nij}(copy(parent(data)))

function Base.size(data::IJFH{S, Nij}) where {S, Nij}
    Nv = 1
    Nh = size(parent(data), 4)
    (Nij, Nij, 1, Nv, Nh)
end

"""
    IJFH{S,Nij}(ArrayType, nelements)

Construct a IJFH 2D Spectral DataLayout given the backing `ArrayType`,
quadrature degrees of freedom `Nij × Nij`  ,
and the number of mesh elements `nelements`.
"""
function IJFH{S, Nij}(ArrayType, nelements) where {S, Nij}
    T = eltype(ArrayType)
    IJFH{S, Nij}(ArrayType(undef, Nij, Nij, typesize(T, S), nelements))
end

function replace_basetype(data::IJFH{S, Nij}, ::Type{T}) where {S, Nij, T}
    array = parent(data)
    S′ = replace_basetype(eltype(array), T, S)
    return IJFH{S′, Nij}(similar(array, T))
end

Base.length(data::IJFH) = size(parent(data), 4)

@generated function _property_view(
    data::IJFH{S, Nij, A},
    ::Val{Idx},
) where {S, Nij, A, Idx}
    SS = fieldtype(S, Idx)
    T = eltype(A)
    offset = fieldtypeoffset(T, S, Idx)
    nbytes = typesize(T, SS)
    field_byterange = (offset + 1):(offset + nbytes)
    return :(IJFH{$SS, $Nij}(
        @inbounds view(parent(data), :, :, $field_byterange, :)
    ))
end

@inline function Base.getproperty(data::IJFH{S, Nij}, i::Integer) where {S, Nij}
    array = parent(data)
    T = eltype(array)
    SS = fieldtype(S, i)
    offset = fieldtypeoffset(T, S, i)
    nbytes = typesize(T, SS)
    dataview = @inbounds view(array, :, :, (offset + 1):(offset + nbytes), :)
    IJFH{SS, Nij}(dataview)
end

@inline function slab(data::IJFH{S, Nij}, h::Integer) where {S, Nij}
    @boundscheck (1 <= h <= length(data)) || throw(BoundsError(data, (h,)))
    dataview = @inbounds view(parent(data), :, :, :, h)
    IJF{S, Nij}(dataview)
end

@inline function slab(data::IJFH{S, Nij}, v::Integer, h::Integer) where {S, Nij}
    @boundscheck (v >= 1 && 1 <= h <= length(data)) ||
                 throw(BoundsError(data, (v, h)))
    dataview = @inbounds view(parent(data), :, :, :, h)
    IJF{S, Nij}(dataview)
end

# ==================
# Data1D DataLayout
# ==================

Base.length(data::Data1D) = size(parent(data), 3)

"""
    IFH{S, Ni, A} <: Data1D{S, Ni}

Backing `DataLayout` for 1D spectral element slabs.

Element nodal point (I) data is contiguous for each datatype `S` struct field (F),
for each 1D mesh element (H).
"""
struct IFH{S, Ni, A} <: Data1D{S, Ni}
    array::A
end

function IFH{S, Ni}(array::AbstractArray{T, 3}) where {S, Ni, T}
    @assert size(array, 1) == Ni
    check_basetype(T, S)
    IFH{S, Ni, typeof(array)}(array)
end

function replace_basetype(data::IFH{S, Ni}, ::Type{T}) where {S, Ni, T}
    array = parent(data)
    S′ = replace_basetype(eltype(array), T, S)
    return IFH{S′, Ni}(similar(array, T))
end

"""
    IFH{S,Ni}(ArrayType, nelements)

Construct a IFH 1D Spectral DataLayout given the backing `ArrayType`,
quadrature degrees of freedom `Ni`  ,
and the number of mesh elements `nelements`.
"""
function IFH{S, Ni}(ArrayType, nelements) where {S, Ni}
    T = eltype(ArrayType)
    IFH{S, Ni}(ArrayType(undef, Ni, typesize(T, S), nelements))
end

rebuild(data::IFH{S, Ni}, array::AbstractArray{T, 3}) where {S, Ni, T} =
    IFH{S, Ni}(array)

Base.copy(data::IFH{S, Ni}) where {S, Ni} = IFH{S, Ni}(copy(parent(data)))

function Base.size(data::IFH{S, Ni}) where {S, Ni}
    Nv = 1
    Nh = size(parent(data), 3)
    (Ni, 1, 1, Nv, Nh)
end

@inline function slab(data::IFH{S, Ni}, h::Integer) where {S, Ni}
    @boundscheck (1 <= h <= length(data)) || throw(BoundsError(data, (h,)))
    dataview = @inbounds view(parent(data), :, :, h)
    IF{S, Ni}(dataview)
end
slab(data::IFH, v::Integer, h::Integer) = slab(data, h)

@generated function _property_view(
    data::IFH{S, Ni, A},
    ::Val{Idx},
) where {S, Ni, A, Idx}
    SS = fieldtype(S, Idx)
    T = eltype(A)
    offset = fieldtypeoffset(T, S, Idx)
    nbytes = typesize(T, SS)
    field_byterange = (offset + 1):(offset + nbytes)
    return :(IFH{$SS, $Ni}(
        @inbounds view(parent(data), :, $field_byterange, :)
    ))
end

@inline function Base.getproperty(data::IFH{S, Ni}, f::Integer) where {S, Ni}
    array = parent(data)
    T = eltype(array)
    SS = fieldtype(S, f)
    offset = fieldtypeoffset(T, S, f)
    nbytes = typesize(T, SS)
    dataview = @inbounds view(array, :, (offset + 1):(offset + nbytes), :)
    IFH{SS, Ni}(dataview)
end

# ======================
# DataSlab2D DataLayout
# ======================

@propagate_inbounds function Base.getindex(
    slab::DataSlab2D{S},
    I::CartesianIndex,
) where {S}
    slab[I[1], I[2]]
end

@propagate_inbounds function Base.setindex!(
    slab::DataSlab2D{S},
    val,
    I::CartesianIndex,
) where {S}
    slab[I[1], I[2]] = val
end

Base.size(::DataSlab2D{S, Nij}) where {S, Nij} = (Nij, Nij, 1, 1, 1)
Base.axes(::DataSlab2D{S, Nij}) where {S, Nij} = (SOneTo(Nij), SOneTo(Nij))

@inline function slab(data::DataSlab2D, h)
    @boundscheck (h >= 1) || throw(BoundsError(data, (h,)))
    data
end

@inline function slab(data::DataSlab2D, v, h)
    @boundscheck (v >= 1 && h >= 1) || throw(BoundsError(data, (v, h)))
    data
end

"""
    IJF{S, Nij, A} <: DataSlab2D{S, Nij}

Backing `DataLayout` for 2D spectral element slab data.

Nodal element data (I,J) are contiguous for each `S` datatype struct field (F) for a single element slab.

A `DataSlab2D` view can be returned from other `Data2D` objects by calling `slab(data, idx...)`.
"""
struct IJF{S, Nij, A} <: DataSlab2D{S, Nij}
    array::A
end

function IJF{S, Nij}(array::AbstractArray{T, 3}) where {S, Nij, T}
    @assert size(array, 1) == Nij
    @assert size(array, 2) == Nij
    check_basetype(T, S)
    IJF{S, Nij, typeof(array)}(array)
end

function replace_basetype(data::IJF{S, Nij}, ::Type{T}) where {S, Nij, T}
    array = parent(data)
    S′ = replace_basetype(eltype(array), T, S)
    return IJF{S′, Nij}(similar(array, T))
end

function Base.size(data::IJF{S, Nij}) where {S, Nij}
    return (Nij, Nij, 1, 1, 1)
end

@generated function _property_view(
    data::IJF{S, Nij, A},
    ::Val{Idx},
) where {S, Nij, A, Idx}
    SS = fieldtype(S, Idx)
    T = eltype(A)
    offset = fieldtypeoffset(T, S, Idx)
    nbytes = typesize(T, SS)
    field_byterange = (offset + 1):(offset + nbytes)
    return :(IJF{$SS, $Nij}(
        @inbounds view(parent(data), :, :, $field_byterange)
    ))
end

@inline function Base.getproperty(data::IJF{S, Nij}, i::Integer) where {S, Nij}
    array = parent(data)
    T = eltype(array)
    SS = fieldtype(S, i)
    offset = fieldtypeoffset(T, S, i)
    nbytes = typesize(T, SS)
    dataview = @inbounds view(array, :, :, (offset + 1):(offset + nbytes))
    IJF{SS, Nij}(dataview)
end

@inline function Base.getindex(
    data::IJF{S, Nij},
    i::Integer,
    j::Integer,
) where {S, Nij}
    @boundscheck (1 <= i <= Nij && 1 <= j <= Nij) ||
                 throw(BoundsError(data, (i, j)))
    dataview = @inbounds view(parent(data), i, j, :)
    get_struct(dataview, S)
end

@inline function Base.setindex!(
    data::IJF{S, Nij},
    val,
    i::Integer,
    j::Integer,
) where {S, Nij}
    @boundscheck (1 <= i <= Nij && 1 <= j <= Nij) ||
                 throw(BoundsError(data, (i, j)))
    dataview = @inbounds view(parent(data), i, j, :)
    set_struct!(dataview, convert(S, val))
end

# ======================
# DataSlab1D DataLayout
# ======================

@propagate_inbounds function Base.getindex(slab::DataSlab1D, I::CartesianIndex)
    slab[I[1]]
end

@propagate_inbounds function Base.setindex!(
    slab::DataSlab1D,
    val,
    I::CartesianIndex,
)
    slab[I[1]] = val
end

function Base.size(::DataSlab1D{<:Any, Ni}) where {Ni}
    return (Ni, 1, 1, 1, 1)
end
Base.axes(::DataSlab1D{S, Ni}) where {S, Ni} = (SOneTo(Ni),)
Base.lastindex(::DataSlab1D{S, Ni}) where {S, Ni} = Ni

@inline function slab(data::DataSlab1D, h)
    @boundscheck (h >= 1) || throw(BoundsError(data, (h,)))
    data
end

@inline function slab(data::DataSlab1D, v, h)
    @boundscheck (v >= 1 && h >= 1) || throw(BoundsError(data, (v, h)))
    data
end

"""
    IF{S, Ni, A} <: DataSlab1D{S, Ni}

Backing `DataLayout` for 1D spectral element slab data.

Nodal element data (I) are contiguous for each `S` datatype struct field (F) for a single element slab.

A `DataSlab1D` view can be returned from other `Data1D` objects by calling `slab(data, idx...)`.
"""
struct IF{S, Ni, A} <: DataSlab1D{S, Ni}
    array::A
end

function IF{S, Ni}(array::AbstractArray{T, 2}) where {S, Ni, T}
    @assert size(array, 1) == Ni
    check_basetype(T, S)
    IF{S, Ni, typeof(array)}(array)
end

function replace_basetype(data::IF{S, Ni}, ::Type{T}) where {S, Ni, T}
    array = parent(data)
    S′ = replace_basetype(eltype(array), T, S)
    return IF{S′, Ni}(similar(array, T))
end

@generated function _property_view(
    data::IF{S, Ni, A},
    ::Val{Idx},
) where {S, Ni, A, Idx}
    SS = fieldtype(S, Idx)
    T = eltype(A)
    offset = fieldtypeoffset(T, S, Idx)
    nbytes = typesize(T, SS)
    field_byterange = (offset + 1):(offset + nbytes)
    return :(IF{$SS, $Ni}(@inbounds view(parent(data), :, $field_byterange)))
end

@inline function Base.getproperty(data::IF{S, Ni}, f::Integer) where {S, Ni}
    array = parent(data)
    T = eltype(array)
    SS = fieldtype(S, f)
    offset = fieldtypeoffset(T, S, f)
    len = typesize(T, SS)
    dataview = @inbounds view(array, :, (offset + 1):(offset + len))
    IF{SS, Ni}(dataview)
end

@inline function Base.getindex(data::IF{S, Ni}, i::Integer) where {S, Ni}
    @boundscheck (1 <= i <= Ni) || throw(BoundsError(data, (i,)))
    dataview = @inbounds view(parent(data), i, :)
    get_struct(dataview, S)
end

@inline function Base.setindex!(data::IF{S, Ni}, val, i::Integer) where {S, Ni}
    @boundscheck (1 <= i <= Ni) || throw(BoundsError(data, (i,)))
    dataview = @inbounds view(parent(data), i, :)
    set_struct!(dataview, convert(S, val))
end

# ======================
# DataColumn DataLayout
# ======================

Base.length(data::DataColumn) = size(parent(data), 1)
Base.size(data::DataColumn) = (1, 1, 1, length(data), 1)

"""
    VF{S, A} <: DataColumn{S}

Backing `DataLayout` for 1D FV column data.

Column level data (V) are contiguous for each `S` datatype struct field (F).

A `DataColumn` view can be returned from other `Data1DX`, `Data2DX` objects by calling `column(data, idx...)`.
"""
struct VF{S, A} <: DataColumn{S}
    array::A
end

function VF{S}(array::AbstractArray{T, 2}) where {S, T}
    check_basetype(T, S)
    VF{S, typeof(array)}(array)
end

function VF{S}(array::AbstractVector{T}) where {S, T}
    @assert typesize(T, S) == 1
    VF{S}(reshape(array, (:, 1)))
end

function VF{S}(ArrayType, nelements) where {S}
    T = eltype(ArrayType)
    VF{S}(ArrayType(undef, nelements, typesize(T, S)))
end

function replace_basetype(data::VF{S}, ::Type{T}) where {S, T}
    array = parent(data)
    S′ = replace_basetype(eltype(array), T, S)
    return VF{S′}(similar(array, T))
end

Base.copy(data::VF{S}) where {S} = VF{S}(copy(parent(data)))
Base.lastindex(data::VF) = length(data)

@generated function _property_view(data::VF{S, A}, ::Val{Idx}) where {S, A, Idx}
    SS = fieldtype(S, Idx)
    T = eltype(A)
    offset = fieldtypeoffset(T, S, Idx)
    nbytes = typesize(T, SS)
    field_byterange = (offset + 1):(offset + nbytes)
    return :(VF{$SS}(@inbounds view(parent(data), :, $field_byterange)))
end

@inline function Base.getproperty(data::VF{S}, i::Integer) where {S}
    array = parent(data)
    T = eltype(array)
    SS = fieldtype(S, i)
    offset = fieldtypeoffset(T, S, i)
    nbytes = typesize(T, SS)
    dataview = @inbounds view(array, :, (offset + 1):(offset + nbytes))
    VF{SS}(dataview)
end

@inline function Base.getindex(data::VF{S}, v::Integer) where {S}
    @boundscheck 1 <= v <= size(parent(data), 1) ||
                 throw(BoundsError(data, (v,)))
    dataview = @inbounds view(parent(data), v, :)
    @inbounds get_struct(dataview, S)
end

@propagate_inbounds function Base.getindex(
    col::DataColumn,
    I::CartesianIndex{5},
)
    col[I[4]]
end

@propagate_inbounds function Base.setindex!(
    col::DataColumn,
    val,
    I::CartesianIndex{5},
)
    col[I[4]] = val
end

@inline function Base.setindex!(data::VF{S}, val, v::Integer) where {S}
    @boundscheck (1 <= v <= length(parent(data))) ||
                 throw(BoundsError(data, (v,)))
    dataview = @inbounds view(parent(data), v, :)
    @inbounds set_struct!(dataview, convert(S, val))
end

@inline function column(data::VF, i, h)
    @boundscheck (i >= 1 && h >= 1) || throw(BoundsError(data, (i, h)))
    data
end

@inline function column(data::VF, i, j, h)
    @boundscheck (i >= 1 && j >= 1 && h >= 1) ||
                 throw(BoundsError(data, (i, j, h)))
    data
end

# ======================
# Data2DX DataLayout
# ======================

"""
    VIJFH{S, Nij, A} <: Data2DX{S, Nij}

Backing `DataLayout` for 2D spectral element slab + extruded 1D FV column data.

Column levels (V) are contiguous for every element nodal point (I, J)
for each `S` datatype struct field (F), for each 2D mesh element slab (H).
"""
struct VIJFH{S, Nij, A} <: Data2DX{S, Nij}
    array::A
end

function VIJFH{S, Nij}(array::AbstractArray{T, 5}) where {S, Nij, T}
    @assert size(array, 2) == size(array, 3) == Nij
    VIJFH{S, Nij, typeof(array)}(array)
end

function replace_basetype(data::VIJFH{S, Nij}, ::Type{T}) where {S, Nij, T}
    array = parent(data)
    S′ = replace_basetype(eltype(array), T, S)
    return VIJFH{S′, Nij}(similar(array, T))
end

function Base.copy(data::VIJFH{S, Nij}) where {S, Nij}
    VIJFH{S, Nij}(copy(parent(data)))
end

function Base.size(data::VIJFH{<:Any, Nij}) where {Nij}
    Nv = size(parent(data), 1)
    Nh = size(parent(data), 5)
    return (Nij, Nij, 1, Nv, Nh)
end

function Base.length(data::VIJFH)
    size(parent(data), 1) * size(parent(data), 5)
end

@generated function _property_view(
    data::VIJFH{S, Nij, A},
    ::Val{Idx},
) where {S, Nij, A, Idx}
    SS = fieldtype(S, Idx)
    T = eltype(A)
    offset = fieldtypeoffset(T, S, Idx)
    nbytes = typesize(T, SS)
    field_byterange = (offset + 1):(offset + nbytes)
    return :(VIJFH{$SS, $Nij}(
        @inbounds view(parent(data), :, :, :, $field_byterange, :)
    ))
end

@propagate_inbounds function Base.getproperty(
    data::VIJFH{S, Nij},
    i::Integer,
) where {S, Nij}
    array = parent(data)
    T = eltype(array)
    SS = fieldtype(S, i)
    offset = fieldtypeoffset(T, S, i)
    nbytes = typesize(T, SS)
    dataview = @inbounds view(array, :, :, :, (offset + 1):(offset + nbytes), :)
    VIJFH{SS, Nij}(dataview)
end

# Note: construct the subarray view directly as optimizer fails in Base.to_indices (v1.7)
@inline function slab(data::VIJFH{S, Nij}, v, h) where {S, Nij}
    array = parent(data)
    Nv = size(array, 1)
    Nh = size(array, 5)
    @boundscheck (1 <= v <= Nv && 1 <= h <= Nh) ||
                 throw(BoundsError(data, (v, h)))
    Nf = size(array, 4)
    dataview = @inbounds SubArray(
        array,
        (
            v,
            Base.Slice(Base.OneTo(Nij)),
            Base.Slice(Base.OneTo(Nij)),
            Base.Slice(Base.OneTo(Nf)),
            h,
        ),
    )
    IJF{S, Nij}(dataview)
end

# Note: construct the subarray view directly as optimizer fails in Base.to_indices (v1.7)
@inline function column(data::VIJFH{S, Nij}, i, j, h) where {S, Nij}
    array = parent(data)
    Nh = size(array, 5)
    @boundscheck (1 <= i <= Nij && 1 <= j <= Nij && 1 <= h <= Nh) ||
                 throw(BoundsError(data, (i, j, h)))
    Nv = size(array, 1)
    Nf = size(array, 4)
    dataview = @inbounds SubArray(
        array,
        (Base.Slice(Base.OneTo(Nv)), i, j, Base.Slice(Base.OneTo(Nf)), h),
    )
    VF{S}(dataview)
end

@inline function level(data::VIJFH{S, Nij}, v) where {S, Nij}
    array = parent(data)
    Nv = size(array, 1)
    @boundscheck (1 <= v <= Nv) || throw(BoundsError(data, (v,)))
    dataview = @inbounds view(array, v, :, :, :, :)
    IJFH{S, Nij}(dataview)
end

@propagate_inbounds function Base.getindex(data::VIJFH, I::CartesianIndex{5})
    data[I[1], I[2], I[4]]
end

@propagate_inbounds function Base.setindex!(
    data::VIJFH,
    val,
    I::CartesianIndex{5},
)
    data[I[1], I[2], I[4]] = val
end

# ======================
# Data1DX DataLayout
# ======================

"""
    VIFH{S, Ni, A} <: Data1DX{S, Ni}

Backing `DataLayout` for 1D spectral element slab + extruded 1D FV column data.

Column levels (V) are contiguous for every element nodal point (I)
for each datatype `S` struct field (F), for each 1D mesh element slab (H).
"""
struct VIFH{S, Ni, A} <: Data1DX{S, Ni}
    array::A
end

function VIFH{S, Ni}(array::AbstractArray{T, 4}) where {S, Ni, T}
    @assert size(array, 2) == Ni
    check_basetype(T, S)
    VIFH{S, Ni, typeof(array)}(array)
end

function replace_basetype(data::VIFH{S, Ni}, ::Type{T}) where {S, Ni, T}
    array = parent(data)
    S′ = replace_basetype(eltype(array), T, S)
    return VIFH{S′, Ni}(similar(array, T))
end

Base.copy(data::VIFH{S, Ni}) where {S, Ni} = VIFH{S, Ni}(copy(parent(data)))

function Base.size(data::VIFH{<:Any, Ni}) where {Ni}
    Nv = size(parent(data), 1)
    Nh = size(parent(data), 4)
    return (Ni, 1, 1, Nv, Nh)
end

function Base.length(data::VIFH)
    size(parent(data), 1) * size(parent(data), 4)
end

@generated function _property_view(
    data::VIFH{S, Ni, A},
    ::Val{Idx},
) where {S, Ni, A, Idx}
    SS = fieldtype(S, Idx)
    T = eltype(A)
    offset = fieldtypeoffset(T, S, Idx)
    nbytes = typesize(T, SS)
    field_byterange = (offset + 1):(offset + nbytes)
    return :(VIFH{$SS, $Ni}(
        @inbounds view(parent(data), :, :, $field_byterange, :)
    ))
end

@inline function Base.getproperty(data::VIFH{S, Ni}, i::Integer) where {S, Ni}
    array = parent(data)
    T = eltype(array)
    SS = fieldtype(S, i)
    offset = fieldtypeoffset(T, S, i)
    nbytes = typesize(T, SS)
    dataview = @inbounds view(array, :, :, (offset + 1):(offset + nbytes), :)
    VIFH{SS, Ni}(dataview)
end

# Note: construct the subarray view directly as optimizer fails in Base.to_indices (v1.7)
@inline function slab(data::VIFH{S, Ni}, v, h) where {S, Ni}
    array = parent(data)
    Nv = size(array, 1)
    Nh = size(array, 4)
    @boundscheck (1 <= v <= Nv && 1 <= h <= Nh) ||
                 throw(BoundsError(data, (v, h)))
    Nf = size(array, 3)
    dataview = @inbounds SubArray(
        array,
        (v, Base.Slice(Base.OneTo(Ni)), Base.Slice(Base.OneTo(Nf)), h),
    )
    IF{S, Ni}(dataview)
end

# Note: construct the subarray view directly as optimizer fails in Base.to_indices (v1.7)
@inline function column(data::VIFH{S, Ni}, i, h) where {S, Ni}
    array = parent(data)
    Nh = size(array, 4)
    @boundscheck (1 <= i <= Ni && 1 <= h <= Nh) ||
                 throw(BoundsError(data, (i, h)))
    Nv = size(array, 1)
    Nf = size(array, 3)
    dataview = @inbounds SubArray(
        array,
        (Base.Slice(Base.OneTo(Nv)), i, Base.Slice(Base.OneTo(Nf)), h),
    )
    VF{S}(dataview)
end

@inline function column(data::VIFH{S, Ni}, i, j, h) where {S, Ni}
    array = parent(data)
    Nh = size(array, 4)
    @boundscheck (1 <= i <= Ni && j == 1 && 1 <= h <= Nh) ||
                 throw(BoundsError(data, (i, j, h)))
    Nv = size(array, 1)
    Nf = size(array, 3)
    dataview = @inbounds SubArray(
        array,
        (Base.Slice(Base.OneTo(Nv)), i, Base.Slice(Base.OneTo(Nf)), h),
    )
    VF{S}(dataview)
end

@inline function level(data::VIFH{S, Nij}, v) where {S, Nij}
    array = parent(data)
    Nv = size(array, 1)
    @boundscheck (1 <= v <= Nv) || throw(BoundsError(data, (v,)))
    dataview = @inbounds view(array, v, :, :, :)
    IFH{S, Nij}(dataview)
end

@propagate_inbounds function Base.getindex(data::VIFH, I::CartesianIndex)
    data[I[1], I[4]]
end

@propagate_inbounds function Base.setindex!(data::VIFH, val, I::CartesianIndex)
    data[I[1], I[4]] = val
end

# =========================================
# Special DataLayouts for regular gridding
# =========================================

struct IH1JH2{S, Nij, A} <: Data2D{S, Nij}
    array::A
end

"""
    IH1JH2{S, Nij}(data::AbstractMatrix{S})

Stores a 2D field in a matrix using a column-major format.
The primary use is for interpolation to a regular grid for ex. plotting / field output.
"""
function IH1JH2{S, Nij}(array::AbstractMatrix{S}) where {S, Nij}
    @assert size(array, 1) % Nij == 0
    @assert size(array, 2) % Nij == 0
    IH1JH2{S, Nij, typeof(array)}(array)
end

Base.copy(data::IH1JH2{S, Nij}) where {S, Nij} =
    IH1JH2{S, Nij}(copy(parent(data)))

function Base.size(data::IH1JH2{S, Nij}) where {S, Nij}
    Nv = 1
    Nh = div(length(parent(data)), Nij * Nij)
    (Nij, Nij, 1, Nv, Nh)
end

Base.length(data::IH1JH2{S, Nij}) where {S, Nij} =
    div(length(parent(data)), Nij * Nij)

function Base.similar(
    data::IH1JH2{S, Nij, A},
    ::Type{Eltype},
) where {S, Nij, A, Eltype}
    array = similar(A, Eltype)
    return IH1JH2{Eltype, Nij}(array)
end

@inline function slab(data::IH1JH2{S, Nij}, h::Integer) where {S, Nij}
    N1, N2 = size(parent(data))
    n1 = div(N1, Nij)
    n2 = div(N2, Nij)
    z2, z1 = fldmod(h - 1, n1)
    @boundscheck (1 <= h <= n1 * n2) || throw(BoundsError(data, (h,)))
    dataview =
        @inbounds view(parent(data), Nij * z1 .+ (1:Nij), Nij * z2 .+ (1:Nij))
    return dataview
end

struct IV1JH2{S, Ni, A} <: Data1DX{S, Ni}
    array::A
end

"""
    IV1JH2{S, Ni}(data::AbstractMatrix{S})

Stores values from an extruded 1D spectral field in a matrix using a column-major format.
The primary use is for interpolation to a regular grid for ex. plotting / field output.
"""
function IV1JH2{S, Ni}(array::AbstractMatrix{S}) where {S, Ni}
    @assert size(array, 2) % Ni == 0
    IV1JH2{S, Ni, typeof(array)}(array)
end

Base.copy(data::IV1JH2{S, Ni}) where {S, Ni} = IV1JH2{S, Ni}(copy(parent(data)))

function Base.size(data::IV1JH2{S, Ni}) where {S, Ni}
    Nv = size(parent(data), 1)
    Nh = div(size(parent(data), 2), Ni)
    (Ni, 1, 1, Nv, Nh)
end

Base.length(data::IV1JH2{S, Ni}) where {S, Ni} = div(length(parent(data)), Ni)

function Base.similar(
    data::IV1JH2{S, Ni, A},
    ::Type{Eltype},
) where {S, Ni, A, Eltype}
    array = similar(A, Eltype)
    return IV1JH2{Eltype, Ni}(array)
end

@inline function slab(data::IV1JH2{S, Ni}, v::Integer, h::Integer) where {S, Ni}
    N1, N2 = size(parent(data))
    n1 = N1
    n2 = div(N2, Ni)
    _, z2 = fldmod(h - 1, n2)
    @boundscheck (1 <= v <= n1) && (1 <= h <= n2) ||
                 throw(BoundsError(data, (v, h)))
    dataview = @inbounds view(parent(data), v, Ni * z2 .+ (1:Ni))
    return dataview
end

# broadcast machinery
include("broadcast.jl")

# GPU method specializations
include("cuda.jl")

end # module


text:
@testset "device arrays" begin

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

@testset "constructors" begin
    # inner constructors
    let
        p = Ptr{Int}(C_NULL)
        @on_device CuDeviceArray{Int,1}((1,), $p)
    end

    # outer constructors
    for I in [Int32,Int64]
        a = I(1)
        b = I(2)
        p = Ptr{I}(C_NULL)

        # not parameterized
        @on_device CuDeviceArray($b, $p)
        @on_device CuDeviceArray(($b,), $p)
        @on_device CuDeviceArray(($b,$a), $p)

        # partially parameterized
        @on_device CuDeviceArray{$I}($b, $p)
        @on_device CuDeviceArray{$I}(($b,), $p)
        @on_device CuDeviceArray{$I}(($a,$b), $p)

        # fully parameterized
        @on_device CuDeviceArray{$I,1}($b, $p)
        @on_device CuDeviceArray{$I,1}(($b,), $p)
        @test_throws ErrorException @on_device CuDeviceArray{$I,1}(($a,$b), $p)
        @test_throws ErrorException @on_device CuDeviceArray{$I,2}($b, $p)
        @test_throws ErrorException @on_device CuDeviceArray{$I,2}(($b,), $p)
        @on_device CuDeviceArray{$I,2}(($a,$b), $p)
    end
end

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

@testset "basics" begin     # argument passing, get and setindex, length
    dims = (16, 16)
    len = prod(dims)

    @eval function array_copy(input::CuDeviceArray{Float32}, output::CuDeviceArray{Float32})
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x

        if i <= length(input)
            output[i] = Float64(input[i])   # force conversion upon setindex!
        end

        return nothing
    end

    input = round.(rand(Float32, dims) * 100)

    input_dev = CuArray(input)
    output_dev = similar(input_dev)

    @cuda (1,len) array_copy(input_dev, output_dev)
    output = Array(output_dev)
    @test input ≈ output
end

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

@testset "bounds checking" begin
    @eval function array_oob_1d(array)
        return array[1]
    end

    # NOTE: these tests verify that bounds checking is _disabled_ (see #4)

    ir = sprint(io->CUDAnative.code_llvm(io, array_oob_1d, (CuDeviceArray{Int,1},)))
    @test !contains(ir, "trap")

    @eval function array_oob_2d(array)
        return array[1, 1]
    end

    ir = sprint(io->CUDAnative.code_llvm(io, array_oob_2d, (CuDeviceArray{Int,2},)))
    @test !contains(ir, "trap")
end

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

@testset "views" begin
    @eval function array_view(array)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x

        sub = view(array, 2:length(array)-1)
        if i <= length(sub)
            sub[i] = i
        end

        return nothing
    end

    array = zeros(Int64, 100)
    array_dev = CuArray(array)

    sub = view(array, 2:length(array)-1)
    for i in 1:length(sub)
        sub[i] = i
    end

    @cuda (100, 1) array_view(array_dev)
    @test array == Array(array_dev)
end

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

end


text:
# file xim_haircut.jl
# is a part of the xim module
# contains code to do haircut to images

#*************************************************
function xim_haircut_testsum(x,y)
  c=x+y
  c
end

#*************************************************
#*************************************************
function xim_hotHaircut(im::Array{Float64,2},r1::Int16,r2::Int16,ratio1::Float64,rd::Int16,ratio2::Float64)
    # cuts exceptionally hot pixels (compared to a local average parametrized by r1 r2 rd ratio1 and ratio2)
    # this function uses the following functions contained in this module...
    #    ...   weigthed_average, rshape1, rshape2, pad_im
    #
    # detailed description
    #

    si      = size(im)
    l       = Int16(2*rd+1)
    a       = Array{Float64,2}(undef,l,1); a = LinRange(-rd,rd,l)

    F1      = Array{Float64,2}(undef,2*rd+1,2*rd+1); F1      = rshape1(rd)
    avmap   = Array{Float64,2}(undef,si[1],si[2]);   avmap   = weigthed_average(Array{Float64,2}(deepcopy(im)),F1)
    ratmap  = Array{Float64,2}(undef,si[1],si[2]);   ratmap  = im./avmap
    badpix  = Array{Bool,2}(undef,si[1],si[2]);      badpix  = ratmap .>ratio2
    goodpix = ones(Bool,si[1],si[2]) .- badpix

    F2      = Array{Float64,2}(undef,2*r1+1,2*r2+1); F2      = rshape2(r1,r2)
    weight  = Array{Float64,2}(undef,si[1],si[2]);   weight  = weigthed_average(Array{Float64,2}(goodpix),F2)
    goodim  = Array{Float64,2}(undef,si[1],si[2]);   goodim  =  deepcopy(im) .* goodpix

    am      = Array{Float64,2}(undef,si[1],si[2]);   am      = weigthed_average(goodim,F2)
    bm      = Array{Float64,2}(undef,si[1],si[2]);   bm      = im ./ am .* weight

    hothair = Array{Float64,2}(undef,si[1],si[2]);   hothair = bm .> ratio1
    ok      = ones(Bool,si[1],si[2]) .- hothair

    new_im  = Array{Float64,2}(undef,si[1],si[2])
      new_im  = ok .* deepcopy(im) + hothair .* am ./ weight

    return new_im
  end

#*************************************************
function xim_coldHaircut(im,r1,r2,ratio1,rd,ratio2)

  im[50:55,:] .= 0.2
  im[:,50:55] .= 0.8
  s = subfunction2(40)
  im = subfunction4(im,s)
  return im
end

#*************************************************
#*************************************************
function weigthed_average(im::Array{Float64,2},w::Array{Float64,2})
  si = size(im)
  sw = size(w)
  padded_im = pad_im(  im , Int16((sw[1]-1)/2) , Int16((sw[2]-1)/2))

  am = zeros(Float64,si[1],si[2])
  for lili = 1:sw[1]
    for coco = 1:sw[2]
      am = am .+ w[lili,coco].*padded_im[lili:lili+si[1]-1,coco:coco+si[2]-1];
    end
  end

  return am
end
#*************************************************

#*************************************************
#*************************************************
function rshape1(R::Int16)

  p  = Array{Float64,2}(undef,7,1); p = LinRange(0.0,R,7)
  f  = Array{Float64,2}(undef,7,1); f = [0.0,0.2,1.0,1.0,0.5,0.2,0.0]
  s  = Float64((R/6)^2)
  l  = Int16(2*R+1)
  a  = Array{Float64,2}(undef,l,1); a = LinRange(-R,R,l)
  rs = sqrt.(repeat(a',l,1).^2 + repeat(a,1,l).^2)

  sh = zeros(Float64,l,l)
  sh += (rs.>=p[1]).*(rs.<p[2]) .* (f[1].+(rs.-p[1]*(f[2]-f[1])./(p[2]-p[1]))./s)
  sh += (rs.>=p[2]).*(rs.<p[3]) .* (f[2].+(rs.-p[2]*(f[3]-f[2])./(p[3]-p[2]))./s)
  sh += (rs.>=p[3]).*(rs.<p[4]) .* (f[3].+(rs.-p[3]*(f[4]-f[3])./(p[4]-p[3]))./s)
  sh += (rs.>=p[4]).*(rs.<p[5]) .* (f[4].+(rs.-p[4]*(f[5]-f[4])./(p[5]-p[4]))./s)
  sh += (rs.>=p[5]).*(rs.<p[6]) .* (f[5].+(rs.-p[5]*(f[6]-f[5])./(p[6]-p[5]))./s)
  sh += (rs.>=p[6]).*(rs.<p[7]) .* (f[6].+(rs.-p[6]*(f[7]-f[6])./(p[7]-p[6]))./s)

  sh = sh/sum(sh)
  return sh
end
#*************************************************

#*************************************************
#*************************************************
function rshape2(R1::Int16,R2::Int16)

  p  = Array{Float64,2}(undef,7,1); p = LinRange(R1,R2,7)
  f  = Array{Float64,2}(undef,7,1); f = [0.0,0.2,1.0,1.0,0.5,0.2,0.0]
  s  = Float64(((R2-R1)/6)^2)
  l  = Int16(2*R2+1)
  a  = Array{Float64,2}(undef,l,1); a = LinRange(-R2,R2,l)
  rs = sqrt.(repeat(a',l,1).^2 + repeat(a,1,l).^2)

  sh = zeros(Float64,l,l)
  sh += (rs.>=p[1]).*(rs.<p[2]) .* (f[1].+(rs.-p[1]*(f[2]-f[1])./(p[2]-p[1]))./s)
  sh += (rs.>=p[2]).*(rs.<p[3]) .* (f[2].+(rs.-p[2]*(f[3]-f[2])./(p[3]-p[2]))./s)
  sh += (rs.>=p[3]).*(rs.<p[4]) .* (f[3].+(rs.-p[3]*(f[4]-f[3])./(p[4]-p[3]))./s)
  sh += (rs.>=p[4]).*(rs.<p[5]) .* (f[4].+(rs.-p[4]*(f[5]-f[4])./(p[5]-p[4]))./s)
  sh += (rs.>=p[5]).*(rs.<p[6]) .* (f[5].+(rs.-p[5]*(f[6]-f[5])./(p[6]-p[5]))./s)
  sh += (rs.>=p[6]).*(rs.<p[7]) .* (f[6].+(rs.-p[6]*(f[7]-f[6])./(p[7]-p[6]))./s)

  sh = sh/sum(sh)
  return sh
end
#*************************************************

#*************************************************
#*************************************************
function pad_im(im::Array{Float64,2},n1::Int16,n2::Int16)
  # function to padd image im in all sides by n1 vertically and n3 horizontally
  si  = size(im)
  nim = Array{Float64,2}(undef,si[1]+2*n1,si[2]+2*n2)
  nim[n1+1:n1+si[1],:] = hcat( reverse(im[:,1:n2],dims=2) , im , reverse(im[:,si[2]-n2+1:si[2]],dims=2) )
  nim[1:n1,:] = reverse(deepcopy(nim[n1+1:2*n1,:]),dims=1)
  nim[n1+si[1]+1:2*n1+si[1],:] = reverse(nim[si[1]+1:n1+si[1],:],dims=1)
  return nim
end
#*************************************************

max_stars_repo_name: NelsonMCLopes/xim

max_issues_repo_name: NelsonMCLopes/xim

max_forks_repo_name: NelsonMCLopes/xim


text:
int main() {
   int =;
   return 1;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
function ga_v2(
        nth_param_set::Int64,
        max_generation::Int64,
        n_population::Int64,
        n_children::Int64,
        n_gene::Int64,
        allowable_error::Float64)::Tuple{Array{Float64, 1}, Float64}
    if n_population < n_gene + 2
        error("n_population must be larger than $(n_gene+2)")
    end

    N_iter::Int64 = 1
    N0::Vector{Float64} = zeros(3*n_population)

    population = get_initial_population(n_population, n_gene)
    N0[1] = population[1, end]
    println(
        @sprintf(
            "Generation%d: Best Fitness = %.6e", 1, population[1, end]
        )
    )
    flush(stdout)

    best_indiv = decode_gene2val(population[1, 1:n_gene])
    best_fitness = population[1, end]

    f = open("./fitparam/$nth_param_set/fit_param1.dat", "w")
    for val in best_indiv
        write(f, @sprintf("%.6e\n", val))
    end
    close(f)

    open("./fitparam/$nth_param_set/generation.dat", "w") do f
        write(f, @sprintf("%d", 1))
    end

    open("./fitparam/$nth_param_set/best_fitness.dat", "w") do f
        write(f, @sprintf("%.6e", best_fitness))
    end

    if population[1, end] <= allowable_error
        best_indiv = decode_gene2val(population[1, 1:n_gene])
        best_fitness = population[1, end]

        return best_indiv, best_fitness
    end

    generation::Int64 = 2
    while generation <= max_generation
        ip = randperm(n_population)[1:n_gene+2]
        population = converging!(ip, population, n_population, n_gene)
        population = local_search!(ip, population, n_population, n_children, n_gene)
        if N_iter > 1
            for _ in 1:N_iter
                ip = randperm(n_population)[1:n_gene+2]
                population = converging!(ip, population, n_population, n_gene)
            end
        end

        if generation%length(N0) == 0
            N0[end] = population[1, end]
            if N0[1] == N0[end]
                N_iter *= 2
            else
                N_iter = 1
            end
        else
            N0[generation%length(N0)] = population[1, end]
        end

        println(
            @sprintf(
                "Generation%d: Best Fitness = %.6e", generation, population[1, end]
            )
        )
        flush(stdout)
        best_indiv = decode_gene2val(population[1, 1:n_gene])
        if population[1, end] < best_fitness
            f = open(
                "./fitparam/$nth_param_set/fit_param$generation.dat", "w"
            )
            for val in best_indiv
                write(f, @sprintf("%.6e\n", val))
            end
            close(f)

            open("./fitparam/$nth_param_set/generation.dat", "w") do f
                write(f, @sprintf("%d", generation))
            end
        end
        best_fitness = population[1, end]

        open("./fitparam/$nth_param_set/best_fitness.dat", "w") do f
            write(f, @sprintf("%.6e", best_fitness))
        end

        if population[1, end] <= allowable_error
            best_indiv = decode_gene2val(population[1, 1:n_gene])
            best_fitness = population[1, end]

            return best_indiv, best_fitness
        end

        open("./fitparam/$nth_param_set/count_num.dat", "w") do f
            write(f, @sprintf("%d", generation))
        end

        generation += 1
    end

    best_indiv = decode_gene2val(population[1, 1:n_gene])
    best_fitness = population[1, end]

    return best_indiv, best_fitness
end

function ga_v2_continue(
        nth_param_set::Int64,
        max_generation::Int64,
        n_population::Int64,
        n_children::Int64,
        n_gene::Int64,
        allowable_error::Float64,
        p0_bounds::Vector{Float64})::Tuple{Array{Float64, 1}, Float64}
    if n_population < n_gene + 2
        error("n_population must be larger than $(n_gene+2)")
    end

    N_iter::Int64 = 1
    N0::Vector{Float64} = zeros(3*n_population)

    count::Int64 = readdlm(
        "./fitparam/$nth_param_set/count_num.dat"
    )[1, 1]
    best_generation::Int64 = readdlm(
        "./fitparam/$nth_param_set/generation.dat"
    )[1, 1]
    best_indiv::Vector{Float64} = readdlm(
        @sprintf(
            "./fitparam/%d/fit_param%d.dat", nth_param_set, best_generation
        )
    )[:, 1]
    best_indiv_gene::Vector{Float64} = encode_val2gene(best_indiv)
    best_fitness::Float64 = objective(best_indiv_gene)

    population = get_initial_population_continue(
        nth_param_set, n_population, n_gene, p0_bounds
    )
    if best_fitness < population[1, end]
        for i in 1:n_gene
            @inbounds population[1, i] = best_indiv_gene[i]
        end
        population[1, end] = best_fitness
    else
        best_indiv = decode_gene2val(population[1, 1:n_gene])
        best_fitness = population[1, end]
        open("./fitparam/$nth_param_set/fit_param$count.dat", "w") do f
            for i=1:n_gene
                write(f, @sprintf("%.6e", best_indiv[i]))
            end
        end
    end

    N0[1] = population[1, end]

    println(
        @sprintf(
            "Generation%d: Best Fitness = %.6e", count + 1, population[1, end]
        )
    )
    flush(stdout)

    if population[1, end] <= allowable_error
        best_indiv = decode_gene2val(population[1, 1:n_gene])
        best_fitness = population[1, end]

        return best_indiv, best_fitness
    end

    generation::Int64 = 2 + count
    while generation <= max_generation
        ip = randperm(n_population)[1:n_gene+2]
        population = converging!(ip, population, n_population, n_gene)
        population = local_search!(ip, population, n_population, n_children, n_gene)
        if N_iter > 1
            for _ in 1:N_iter
                ip = randperm(n_population)[1:n_gene+2]
                population = converging!(ip, population, n_population, n_gene)
            end
        end

        if generation%length(N0) == 0
            N0[end] = population[1, end]
            if N0[1] == N0[end]
                N_iter *= 2
            else
                N_iter = 1
            end
        else
            N0[generation%length(N0)] = population[1, end]
        end

        println(
            @sprintf(
                "Generation%d: Best Fitness = %.6e", generation, population[1, end]
            )
        )
        flush(stdout)
        best_indiv = decode_gene2val(population[1, 1:n_gene])
        if population[1, end] < best_fitness
            f = open(
                @sprintf(
                    "./fitparam/%d/fit_param%d.dat", nth_param_set, generation
                ), "w"
            )
            for val in best_indiv
                write(f, @sprintf("%.6e\n", val))
            end
            close(f)

            open("./fitparam/$nth_param_set/generation.dat", "w") do f
                write(f, @sprintf("%d", generation))
            end
        end
        best_fitness = population[1, end]

        open("./fitparam/$nth_param_set/best_fitness.dat", "w") do f
            write(f, @sprintf("%.6e", best_fitness))
        end

        if population[1, end] <= allowable_error
            best_indiv = decode_gene2val(population[1, 1:n_gene])
            best_fitness = population[1, end]

            return best_indiv, best_fitness
        end

        open("./fitparam/$nth_param_set/count_num.dat", "w") do f
            write(f, @sprintf("%d", generation))
        end

        generation += 1
    end

    best_indiv = decode_gene2val(population[1, 1:n_gene])
    best_fitness = population[1, end]

    return best_indiv, best_fitness
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder
using Pkg
name = "CImPlot"
version = v"0.3.0" # tracking tags of epezent/implot

# Collection of sources required to build CImGui
sources = [
    GitSource("https://github.com/cimgui/cimplot.git",
              "929f61d27b8b4e3c899b2a386679a7f4a82826ce"),

    GitSource("https://github.com/epezent/implot.git",
              "b95b23fec80dfe547da8905b570a9db30becf904"),

    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
mv cmake/CMakeLists.txt ./cimplot/
mv implot cimplot/
cd cimplot/
mkdir build && cd build
cmake .. -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release
make -j${nproc}
make install
mkdir -p ${prefix}/share/licenses/${SRC_NAME}
mv ../LICENSE ${prefix}/share/licenses/${SRC_NAME}/LICENSE_cimplot
mv ../implot/LICENSE ${prefix}/share/licenses/${SRC_NAME}/LICENSE_implot
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libcimplot", :libcimplot),
    FileProduct("share/compile_commands.json", :compile_commands)
]

# Dependencies that must be installed before this package can be built
dependencies = [
                Dependency(PackageSpec(name="CImGui_jll", version=v"1.77.0")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
println("\n\n\nStarting runtests.jl $(join(ARGS, " ")) ...")

using Tests
using KShiftsClustering

getdata(n) = 10*rand(1, n) .+ 0.5

centers = kshifts(getdata(1_000_000), 10)
@test all(round.(Int, sort(vec(centers)))  .==  collect(1:10))

centers = kshifts(getdata(1000), 10)
for i = 1:100
    kshifts!(centers, getdata(1000))
end
@test all(round.(Int, sort(vec(centers)))  .==  collect(1:10))


text:
using Test, Spotify

#=
   Test behaviour can be further customized by providing
   additional arguments to Pkg as shown here:
   Pkg.test("HelloWorld"; test_args = ["foo", "baz"])
=#

# Argument for quiet mode
quiet = length(ARGS) > 0 && ARGS[1] == "q"

errors = false

all_tests = ["test_int_format_strings.jl",
             "generalize_calls.jl",
             "test_playlists.jl",
             "test_users.jl",
             "test_search.jl",
             "test_markets.jl",
             "test_genres.jl",
             "test_player.jl",
             "test_albums.jl",
             "test_artist.jl",
             "test_tracks.jl",
             "test_shows.jl",
             "test_browse.jl",
             "test_episodes.jl"]

Spotify.refresh_spotify_credentials()

println("Running full test suite:")

@time for file in all_tests
    try
        include(file)
        println("\t\033[1m\033[32mPASSED\033[0m: $(file)")
    catch e
        println("\t\033[1m\033[31mFAILED\033[0m: $(file)")
        global errors = true

        # Show detailed tracing when not in quiet mode
        if ~quiet
            showerror(stdout, e, backtrace())
            println()
        end
    end
end

if errors
    @warn "Some tests have failed! Check the results summary above."
end


text:
using Plots
x = 1:10; y = rand(10); # These are the plotting data
plot(x,y, label="my label")


text:
# ucal.jl - Wrapper for ICU (International Components for Unicode) library

# Some content of the documentation strings was derived from the ICU header files ucal.h
# (Those portions copyright (C) 1996-2015, International Business Machines Corporation and others)

export UCAL

"""
Defines the following constants for ICU Calendar functions (not exported):
ERA, YEAR, MONTH, WEEK_OF_YEAR, WEEK_OF_MONTH,
DATE, DAY_OF_YEAR, DAY_OF_WEEK, DAY_OF_WEEK_IN_MONTH,
AM_PM, HOUR, HOUR_OF_DAY, MINUTE, SECOND, MILLISECOND,
ZONE_OFFSET, DST_OFFSET, YEAR_WOY, DOW_LOCAL, EXTENDED_YEAR,
JULIAN_DAY, MILLISECONDS_IN_DAY, IS_LEAP_MONTH
"""
module UCAL
for (i,a) in enumerate([
        :ERA,
        :YEAR,
        :MONTH,
        :WEEK_OF_YEAR,
        :WEEK_OF_MONTH,
        :DATE,
        :DAY_OF_YEAR,
        :DAY_OF_WEEK,
        :DAY_OF_WEEK_IN_MONTH,
        :AM_PM,
        :HOUR,
        :HOUR_OF_DAY,
        :MINUTE,
        :SECOND,
        :MILLISECOND,
        :ZONE_OFFSET,
        :DST_OFFSET,
        :YEAR_WOY,
        :DOW_LOCAL,
        :EXTENDED_YEAR,
        :JULIAN_DAY,
        :MILLISECONDS_IN_DAY,
        :IS_LEAP_MONTH
    ])
    @eval const $a = Int32($i - 1)
end
end # module UCAL

const UDate = Float64

macro libcal(s)     ; _libicu(s, iculibi18n, "ucal_")     ; end

export UCalendar
mutable struct UCalendar
    ptr::Ptr{Cvoid}

    function UCalendar(tz::Ptr{UInt16}, len)
        err = Ref{UErrorCode}(0)
        p = ccall(@libcal(open), Ptr{Cvoid},
                  (Ptr{UChar}, Int32, Ptr{UInt8}, Int32, Ptr{UErrorCode}),
                  tz, len, locale[], 0, err)
        self = new(p)
        finalizer(self, close)
        self
    end
    UCalendar(tz::WordStrings) = @preserve tz UCalendar(pointer(tz), ncodeunits(tz))
    UCalendar(tz::Vector{UInt16}) = @preserve tz UCalendar(pointer(tz), length(tz))
    function UCalendar()
        err = Ref{UErrorCode}(0)
        p = ccall(@libcal(open), Ptr{Cvoid},
                  (Ptr{UChar}, Int32, Ptr{UInt8}, Int32, Ptr{UErrorCode}),
                  C_NULL, 0, locale[], 0, err)
        self = new(p)
        finalizer(self, close)
        self
    end
end

UCalendar(timezone::AbstractString) = UCalendar(cvt_utf16(timezone))

close(c::UCalendar) =
    c.ptr == C_NULL || (ccall(@libcal(close), Cvoid, (Ptr{Cvoid},), c.ptr) ; c.ptr = C_NULL)

getnow() = ccall(@libcal(getNow), UDate, ())

function get_millis(cal::UCalendar)
    err = Ref{UErrorCode}(0)
    ccall(@libcal(getMillis), UDate, (Ptr{Cvoid}, Ptr{UErrorCode}), cal.ptr, err)
end

function set_millis!(cal::UCalendar, millis::UDate)
    err = Ref{UErrorCode}(0)
    ccall(@libcal(setMillis), Cvoid,
          (Ptr{Cvoid}, UDate, Ptr{UErrorCode}),
          cal.ptr, millis, err)
end

function set_date!(cal::UCalendar, y::Integer, m::Integer, d::Integer)
    err = Ref{UErrorCode}(0)
    ccall(@libcal(setDate), Cvoid,
          (Ptr{Cvoid}, Int32, Int32, Int32, Ptr{UErrorCode}),
          cal.ptr, y, m-1, d, err)
end

function set_datetime!(cal::UCalendar, y::Integer, mo::Integer, d::Integer,
                     h::Integer, mi::Integer, s::Integer)
    err = Ref{UErrorCode}(0)
    ccall(@libcal(setDateTime), Cvoid,
          (Ptr{Cvoid}, Int32, Int32, Int32, Int32, Int32, Int32, Ptr{UErrorCode}),
          cal.ptr, y, mo-1, d, h, mi, s, err)
end

function clear!(cal::UCalendar)
    err = Ref{UErrorCode}(0)
    ccall(@libcal(clear), Cvoid,
          (Ptr{Cvoid}, Ptr{UErrorCode}),
          cal.ptr, err)
end

function get(cal::UCalendar, field::Int32)
    err = Ref{UErrorCode}(0)
    ccall(@libcal(get), Int32,
          (Ptr{Cvoid},Int32,Ptr{UErrorCode}),
          cal.ptr, field, err)
end
get(cal::UCalendar, fields::Vector{Int32}) = [get(cal,f) for f in fields]

function add!(cal::UCalendar, field::Int32, amount::Integer)
    err = Ref{UErrorCode}(0)
    ccall(@libcal(add), Int32,
          (Ptr{Cvoid},Int32,Int32,Ptr{UErrorCode}),
          cal.ptr, field, amount, err)
end

set!(cal::UCalendar, field::Int32, val::Integer) =
    ccall(@libcal(set), Cvoid, (Ptr{Cvoid}, Int32, Int32), cal.ptr, field, val)

function get_timezone_displayname(cal::UCalendar)
    bufsz = 64
    buf, pnt = _allocate(UInt16, bufsz)
    err = Ref{UErrorCode}(0)
    len = ccall(@libcal(getTimeZoneDisplayName), Int32,
                (Ptr{Cvoid}, Int32, Ptr{UInt8}, Ptr{UChar}, Int32, Ptr{UErrorCode}),
                cal.ptr, 1, locale[], pnt, bufsz, err)
    Str(UTF16CSE, buf[1:len])
end

function get_default_timezone()
    bufsz = 64
    buf, pnt = _allocate(UInt16, bufsz)
    err = Ref{UErrorCode}(0)
    len = ccall(@libcal(getDefaultTimeZone), Int32,
                (Ptr{UChar}, Int32, Ptr{UErrorCode}),
                pnt, bufsz, err)
    Str(UTF16CSE, buf[1:len])
end


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

@doc raw"""EndpointPort represents a Port used by an EndpointSlice

    IoK8sApiDiscoveryV1beta1EndpointPort(;
        appProtocol=nothing,
        name=nothing,
        port=nothing,
        protocol=nothing,
    )

    - appProtocol::String : The application protocol for this port. This field follows standard Kubernetes label syntax. Un-prefixed names are reserved for IANA standard service names (as per RFC-6335 and http://www.iana.org/assignments/service-names). Non-standard protocols should use prefixed names. Default is empty string.
    - name::String : The name of this port. All ports in an EndpointSlice must have a unique name. If the EndpointSlice is dervied from a Kubernetes service, this corresponds to the Service.ports[].name. Name must either be an empty string or pass DNS_LABEL validation: * must be no more than 63 characters long. * must consist of lower case alphanumeric characters or '-'. * must start and end with an alphanumeric character. Default is empty string.
    - port::Int32 : The port number of the endpoint. If this is not specified, ports are not restricted and must be interpreted in the context of the specific consumer.
    - protocol::String : The IP protocol for this port. Must be UDP, TCP, or SCTP. Default is TCP.
"""
mutable struct IoK8sApiDiscoveryV1beta1EndpointPort <: SwaggerModel
    appProtocol::Any # spec type: Union{ Nothing, String } # spec name: appProtocol
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    port::Any # spec type: Union{ Nothing, Int32 } # spec name: port
    protocol::Any # spec type: Union{ Nothing, String } # spec name: protocol

    function IoK8sApiDiscoveryV1beta1EndpointPort(;appProtocol=nothing, name=nothing, port=nothing, protocol=nothing)
        o = new()
        validate_property(IoK8sApiDiscoveryV1beta1EndpointPort, Symbol("appProtocol"), appProtocol)
        setfield!(o, Symbol("appProtocol"), appProtocol)
        validate_property(IoK8sApiDiscoveryV1beta1EndpointPort, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiDiscoveryV1beta1EndpointPort, Symbol("port"), port)
        setfield!(o, Symbol("port"), port)
        validate_property(IoK8sApiDiscoveryV1beta1EndpointPort, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        o
    end
end # type IoK8sApiDiscoveryV1beta1EndpointPort

const _property_map_IoK8sApiDiscoveryV1beta1EndpointPort = Dict{Symbol,Symbol}(Symbol("appProtocol")=>Symbol("appProtocol"), Symbol("name")=>Symbol("name"), Symbol("port")=>Symbol("port"), Symbol("protocol")=>Symbol("protocol"))
const _property_types_IoK8sApiDiscoveryV1beta1EndpointPort = Dict{Symbol,String}(Symbol("appProtocol")=>"String", Symbol("name")=>"String", Symbol("port")=>"Int32", Symbol("protocol")=>"String")
Base.propertynames(::Type{ IoK8sApiDiscoveryV1beta1EndpointPort }) = collect(keys(_property_map_IoK8sApiDiscoveryV1beta1EndpointPort))
Swagger.property_type(::Type{ IoK8sApiDiscoveryV1beta1EndpointPort }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiDiscoveryV1beta1EndpointPort[name]))}
Swagger.field_name(::Type{ IoK8sApiDiscoveryV1beta1EndpointPort }, property_name::Symbol) =  _property_map_IoK8sApiDiscoveryV1beta1EndpointPort[property_name]

function check_required(o::IoK8sApiDiscoveryV1beta1EndpointPort)
    true
end

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


text:
export FFT, u_of_t, uv_nonrotating_frame

"""Fourier transform the timeseries of a simulation in the rotating frame and calculate the quadratures and freqeuncies in the non-rotating frame."""
function FFT(soln_u, soln_t; window = DSP.Windows.hanning)
    "Input: solution object of DifferentialEquation (positions array and corresponding time)
    Output: Fourier transform and frequencies, where window function window was used"
    w = window(length(soln_t))
    dt = soln_t[2]-soln_t[1]

    soln_tuples = Tuple.(zip(soln_u, soln_t))

    fft_u =length(soln_t)/sum(w)*[fft(w.*[u[j] for (u,t) in soln_tuples])|> fftshift for j in 1:length(soln_u[1])];
    fft_f = fftfreq(length(soln_t), 1/dt) |> fftshift

    # normalize fft_u
    return(fft_u / length(fft_f), 2*pi*fft_f)
end

FFT(soln::OrdinaryDiffEq.ODECompositeSolution; window=DSP.Windows.hanning) = FFT(soln.u, soln.t, window=window)

function FFT_analyze(fft_u::Vector{ComplexF64}, fft_f)
    "finds peaks in the spectrum and returns corresponding frequency, amplitude and phase.
    Frequency and phase are corrected according to Huang Dishan, Mechanical Systems and Signal Processing (1995) 9(2), 113–118
    This correction works for a rectangular window."

    # retaining more sigdigits gives more ''spurious'' peaks
    max_indices, mxval = peakprom(round.( abs.(fft_u), sigdigits=3),minprom =1)
    Del = fft_f[2]-fft_f[1] # frequency spacing
    A1= abs.(fft_u)[max_indices]
    df = zeros(length(max_indices))

    # correction to the amplitude and phase of the peak
    for i in 1:length(max_indices)
        if abs.(fft_u)[max_indices[i]-1]<abs.(fft_u)[max_indices[i]+1]
            A2= abs.(fft_u)[max_indices[i]+1]
            df[i] = -Del /(A1[i]/A2+1)
        else
            A2= abs.(fft_u)[max_indices[i]-1]
            df[i] = Del /(A1[i]/A2+1)
        end
    end
    return 2*pi*(fft_f[max_indices]-df),  A1.*2, angle.(fft_u)[max_indices]+pi*df/Del
end

function u_of_t(omegas_peak,As_peak,phis_peak,t)
    "Calculate us or vs as a function of time from the Fourier components."
    N = length(omegas_peak)
    u = zeros(length(t))
    for m in Int(N/2 - mod(N/2,1))+1:N
        u .+= As_peak[m]*cos.(omegas_peak[m].*t.+phis_peak[m]);
    end
    return u
end

function uv_nonrotating_frame(omega_rot,omega_peak,A_u_peak,phi_u_peak,A_v_peak,phi_v_peak)
    "calculates amplitudes and frequencies of the position in the nonrotating frame from the
    amplitudes and frequencies in the rotating frame."
    omega_nr = [omega_rot-omega_peak, omega_rot+omega_peak]
    u_nr = [-A_u_peak*cos(phi_u_peak) + A_v_peak*sin(phi_v_peak); -A_u_peak*cos(phi_u_peak) - A_v_peak*sin(phi_v_peak)]./2
    v_nr = [A_v_peak*cos(phi_v_peak) + A_u_peak*sin(phi_u_peak) ; A_v_peak*cos(phi_v_peak) - A_u_peak*sin(phi_u_peak)]./2
    return omega_nr, u_nr, v_nr
end


text:
# NOTE: Luxor is using JuliaGraphics/Luxor.jl/pull/48
using Luxor, BitBasis

"""
   btable(n, center=O; len=ndigits(n), width=20, title=false)

Draw a table of binary digits with `n` which can be [`BitStr`](@ref) or an `Integer`.
"""
function btable(n::Integer, center = O; len = ndigits(n; base = 2), width = 20, title = false)
    sethue("black")
    if title
        t = Table(2, len, width, width, center)
    else
        t = Table(1, len, width, width, center)
    end

    bstr = string(n, base = 2, pad = len)
    fontsize(width ÷ 2)

    # start plotting
    col = 1
    if title
        for k in 1:size(t, 2)
            text(string(len - k + 1), t[1, k], halign = :center, valign = :middle)
        end
        col += 1
    end

    for j in 1:size(t, 2)
        box(t[col, j], width, width, :stroke)
        text(string(bstr[j]), t[col, j], halign = :center, valign = :middle)
    end

    return t
end

function btable(
    f::Function,
    n::Integer,
    center = O;
    len = ndigits(n; base = 2),
    width = 20,
    title = false,
)
    gsave()
    # before
    before = btable(n, O - (0, 50); width = 50, title = true)
    sethue("blue")
    arrow(Point(0, 10), Point(0, 50); linewidth = 2)
    # after
    after = btable(f(n), O + (0, 80); width = 50, title = false)
    grestore()
    return (before, after)
end

function btable(f::Function, n::BitStr, center = O; width = 50, title = false)
    gsave()
    # before
    before = btable(n, O - (0, 50); width = 50, title = true)
    sethue("blue")
    arrow(Point(0, 10), Point(0, 50); linewidth = 2)
    # after
    after = btable(f(n), O + (0, 80); width = 50, title = false)
    grestore()
    return (before, after)
end

btable(n::BitStr, center = O; width = 50, title = false) =
    btable(n.val, center; len = length(n), width = width, title = title)

# TODO: push this to upstream: Luxor
function swaparrow(
    left::Point,
    right::Point,
    action = :fill;
    linewidth = 1.0,
    headlength = 10,
    headangle = pi / 8,
    target_angle = pi / 3,
    fraction = 1 / 5,
    debug = false,
)
    gsave()
    v = right - left
    delta = fraction * v.x
    C1 = Point(left.x + delta, left.y - delta * tan(target_angle))
    C2 = Point(right.x - delta, right.y - delta * tan(target_angle))
    arrow(left, C1, C2, right, action; linewidth = linewidth, rightarrow = true, leftarrow = true)
    if debug
        line(left, C1, :stroke)
        line(right, C2, :stroke)
    end
    grestore()
end

@png begin
    t, _ = btable(breflect, bit"1011"; width = 50, title = true)
    lineheight = 20
    sethue("red")
    swaparrow(
        t[1, 1] - (0, lineheight),
        t[1, 4] - (0, lineheight),
        :stroke;
        linewidth = 2,
        fraction = 0.25,
        target_angle = deg2rad(60),
    )
    swaparrow(
        t[1, 2] - (0, lineheight),
        t[1, 3] - (0, lineheight),
        :stroke;
        linewidth = 2,
        target_angle = deg2rad(80),
    )
end 400 400 "docs/src/assets/breflect"


text:
# # Load MEDLINE

#md # [![nbviewer](https://img.shields.io/badge/jupyter_notebook-nbviewer-orange.svg)](http://nbviewer.jupyter.org/github/bcbi/BioMedQuery.jl/tree/master/docs/src/notebooks/5_load_medline.ipynb)

# The MEDLINE loader process in BioMedQuery saves the MEDLINE baseline files to a
# MySQL database and saves the raw (xml.gz) and parsed (csv) files to a ```medline```
# directory that will be created in the provided ```output_dir```.
#
# **WARNING:** There are 900+ medline files each with approximately 30,000 articles.
# This process will take hours to run for the full baseline load.
#
# The baseline files can be found [here](ftp://ftp.ncbi.nlm.nih.gov/pubmed/baseline/).

# ### Set Up
# The database and tables must already be created before loading the medline files.
# This process is set up for parallel processing.  To take advantage of this, workers
# can be added before loading the BioMedQuery package using the ```addprocs``` function.

using BioMedQuery.DBUtils
using BioMedQuery.PubMed
using BioMedQuery.Processes

# BioMedQuery has utility functions to create the database and tables. *Note: creating
# the tables using this function will drop any tables that already exist in the target
# database.*

const conn = DBUtils.init_mysql_database("127.0.0.1","root","","test_db", true);
PubMed.create_tables!(conn);

# ### Load a Test File
# As the full medline load is a large operation, it is recommended that a test run
# be completed first.

@time Processes.load_medline!(conn, pwd(), test=true)

# Review the output of this run in MySQL to make sure that it ran as expected.
# Additionally, the sample raw and parsed file should be in the new ```medline```
# directory in the current directory.

# ### Performing a Full Load
# To run a full load, use the same code as above, but do not pass the test variable.
# It is also possible to break up the load by passing which files to start and stop at -
# simply pass ```start_file=n``` and ```end_file=p```.
#
# After loading, it is recommended you add indexes to the tables, the ```add_mysql_keys!```
# function can be used to add a standard set of indexes.

add_mysql_keys!(conn)


text:
using Klara

function plogtarget(x::Vector{Float64})
  s2 = exp(x[1])
  nx = length(x)-1
  -0.5*((x[1]-0.0)^2/9.0+dot(x[2:end], x[2:end])/s2+nx*log(s2))
end

p = BasicContMuvParameter(:p, logtarget=plogtarget)

model = likelihood_model(p, false)

sampler = SliceSampler(1., 5)

mcrange = BasicMCRange(nsteps=100000, burnin=10000)

v0 = Dict(:p=>[10.; zeros(4)])

job = BasicMCJob(model, sampler, mcrange, v0)

run(job)

chain = output(job)

mean(chain)


text:
using TonieBox
using Test

@testset "TonieBox.jl" begin
    # Write your tests here.
end


text:
using POMDPLinter
using Test
import POMDPs: POMDP, solve, DDNOut
import POMDPs: states, stateindex, transition, discount, reward, observation, initialstate, initialobs, gen
import POMDPs
using Random

mightbemissing(x) = ismissing(x) || x

tcall = Meta.parse("f(arg1::T1, arg2::T2)")
@test POMDPLinter.unpack_typedcall(tcall) == (:f, [:arg1, :arg2], [:T1, :T2])

# tests case where types aren't specified
@POMDP_require tan(s) begin
    @req sin(::typeof(s))
    @req cos(::typeof(s))
end

module MyModule
    import POMDPs: POMDP, Solver, statetype, actiontype, states, actions, observations, transition, gen, DDNOut, solve
    import POMDPs
    using Random
    using POMDPLinter

    export CoolSolver, solve

    mutable struct CoolSolver <: Solver end

    p = nothing # to test hygeine
    @POMDP_require solve(s::CoolSolver, p::POMDP) begin
        PType = typeof(p)
        S = statetype(PType)
        A = actiontype(PType)
        @req states(::PType)
        @req actions(::PType)
        @req transition(::PType, ::S, ::A)
        @subreq util2(p)
        s = first(states(p))
        @subreq util1(s)
        a = first(actions(p))
        t_dist = transition(p, s, a)
        @req rand(::AbstractRNG, ::typeof(t_dist))
        @req gen(::DDNOut{:o}, ::PType, ::S, ::A, ::MersenneTwister)
    end

    function POMDPs.solve(s::CoolSolver, problem::POMDP{S,A,O}) where {S,A,O}
        @warn_requirements solve(s, problem)
        reqs = @get_requirements solve(s,problem)
        @assert p==nothing
        return check_requirements(reqs)
    end

    util1(x) = abs(x)

    util2(p::POMDP) = observations(p)
    @POMDP_require util2(p::POMDP) begin
        P = typeof(p)
        @req observations(::P)
    end
end

using Main.MyModule

mutable struct SimplePOMDP <: POMDP{Float64, Bool, Int} end
POMDPs.actions(SimplePOMDP) = [true, false]

POMDPs.discount(::SimplePOMDP) = 0.9

let a = 0.0, f(x) = x^2
    @test @req(f(a, 4)) == @req(f(::typeof(a), ::typeof(4)))
end

reqs = nothing # to check the hygeine of the macro
println("There should be a warning about no @reqs here:")
# 27 minutes has been spent trying to suppress this warning and automate a test for it. If you work more on it, please update this counter. The following things have been tried
# - @test_logs (:warn, "No") @POMDP_requirements ...
# - @capture_err @POMDP_requirements ... # From Suppressor.jl
# - @capture_out @POMDP_requirements ... # From Suppressor.jl
@POMDP_requirements "Warn none" begin
    1+1
end
@test reqs == nothing
@test_throws LoadError macroexpand(Main, quote @POMDP_requirements "Malformed" begin
        @req iterator(typeof(as))
    end
end)

# solve(CoolSolver(), SimplePOMDP())
@test_throws MethodError solve(CoolSolver(), SimplePOMDP())

POMDPs.states(::SimplePOMDP) = [1.4, 3.2, 5.8]
struct SimpleDistribution
    ss::Vector{Float64}
    b::Vector{Float64}
end
POMDPs.transition(p::SimplePOMDP, s::Float64, ::Bool) = SimpleDistribution(states(p), [0.2, 0.2, 0.6])

@test (solve(CoolSolver(), SimplePOMDP()) & false) == false

POMDPs.observations(p::SimplePOMDP) = [1,2,3]

Random.rand(rng::AbstractRNG, d::SimpleDistribution) = sample(rng, d.ss, WeightVec(d.b))
POMDPs.gen(::DDNOut{:o}, m::SimplePOMDP, s, a, rng) = 1

@test solve(CoolSolver(), SimplePOMDP())

@testset "show_requirements" begin
    @test @show_requirements(solve(CoolSolver(), SimplePOMDP()))==true
end

@testset "requirements_info" begin
    @test_skip @requirements_info(CoolSolver(), SimplePOMDP())==true
end

struct A <: POMDP{Int,Bool,Bool} end
struct B <: POMDP{Int, Bool, Bool} end
struct W <: POMDP{Int, Bool, Int} end
@testset "implement" begin

    # broken tests should probably only be fixed when POMDPs 1.0 is released

    @test_throws MethodError length(states(A()))
    @test_throws MethodError stateindex(A(), 1)

    @test !@implemented transition(::A, ::Int, ::Bool)
    POMDPs.transition(::A, s, a) = [s+a]
    @test @implemented transition(::A, ::Int, ::Bool)

    @test !@implemented discount(::A)
    POMDPs.discount(::A) = 0.95
    @test @implemented discount(::A)

    @test_broken !@implemented POMDPs.reward(::A,::Int,::Bool,::Int)
    @test !@implemented POMDPs.reward(::A,::Int,::Bool)
    POMDPs.reward(::A,::Int,::Bool) = -1.0
    @test @implemented POMDPs.reward(::A,::Int,::Bool,::Int)
    @test @implemented POMDPs.reward(::A,::Int,::Bool)

    @test_broken !@implemented POMDPs.observation(::A,::Int,::Bool,::Int)
    @test_broken !@implemented POMDPs.observation(::A,::Bool,::Int)
    POMDPs.observation(::A,::Bool,::Int) = [true, false]
    @test @implemented POMDPs.observation(::A,::Int,::Bool,::Int)
    @test @implemented POMDPs.observation(::A,::Bool,::Int)

    @test_broken !@implemented POMDPs.initialstate(::W)
    @test !@implemented POMDPs.initialobs(::W, ::Int)

    POMDPs.transition(b::B, s::Int, a::Bool) = Deterministic(s+a)
    @test mightbemissing(implemented(gen, Tuple{DDNOut{:sp}, B, Int, Bool, MersenneTwister}))

    reward(b::B, s::Int, a::Bool, sp::Int) = -1.0
    observation(b::B, s::Int, a::Bool, sp::Int) = Deterministic(sp)
    @test mightbemissing(@implemented(gen(::DDNOut{(:sp,:o,:r)}, ::B, ::Int, ::Bool, ::MersenneTwister)))
    @test mightbemissing(@implemented(gen(::DDNOut{(:sp,:o)}, b::B, s::Int, a::Bool, rng::MersenneTwister)))
    @test mightbemissing(@implemented gen(::DDNOut{(:sp,:o,:r)}, b::B, s::Int, a::Bool, rng::MersenneTwister))

    initialstate_distribution(b::B) = Int[1,2,3]
    @test @implemented initialstate(::B, ::MersenneTwister)

    POMDPs.observation(b::B, s::Int) = Bool[s]
    @test @implemented initialobs(::B, ::Int, ::MersenneTwister)
end