text:
using FourierFlows: stepforward!
using Printf
using NPZ, MAT

include("nondimlineartwolayerqg.jl")
set_q! = NondimLinearTwoLayerQG.set_q!
updatevars! = NondimLinearTwoLayerQG.updatevars!
energies = NondimLinearTwoLayerQG.energies
fluxes = NondimLinearTwoLayerQG.fluxes
get_bickley = NondimLinearTwoLayerQG.get_bickley
get_topo = NondimLinearTwoLayerQG.get_topo
get_σ = NondimLinearTwoLayerQG.get_σ

include("common_params_twolay.jl")

for F₁ in all_F₁s, ht in all_hts

println("")
println("F₁ = ", string(F₁), ", ht = ", string(ht), ", topography = ", ttype)
println("")
F₁, ht = float(F₁), float(ht)

if ttype!="rand"
    h, hx, hy = get_topo(gr, ht, kt, type=ttype)
else
    T = eltype(gr)
    nx, ny = gr.nx, gr.ny
    h, hx, hy = zeros(T, nx, ny), zeros(T, nx, ny), zeros(T, nx, ny)
    mfile = matopen(fname_hrand)
    h = read(mfile, "h")
    close(mfile)
    @. h = h*ht*0.5 # ht is rms in random topography. The factor 1/2 is the rms of sin(x)sin(y).
    hh = rfft(h)
    hxh = @. im * gr.kr * hh
    hyh = @. im * gr.l * hh
    hx = irfft(hxh, nx)
    hy = irfft(hyh, nx)
end

prob = NondimLinearTwoLayerQG.Problem(dev;
                                      nx = N,
                                      ny = N,
                                      Lx = L,
                                      Ly = L,
                                      x0 = 0.0,
                                      F₁ = F₁,
                                       δ = δ,
                                       U = U,
                                     Uyy = Uyy,
                                       h = h,
                                      hx = hx,
                                      hy = hy,
                                       α = α,
                                       β = β,
                                      dt = dt,
                                 stepper = stepper,
                                  effort = effort)

sol, clock, params, vars, grid = prob.sol, prob.clock, prob.params, prob.vars, prob.grid
set_q!(prob, q₀)

if PSI_HOVMOLLER thovm, ψ₁hovm, ψ₂hovm = nothing, nothing, nothing; fm = Int(round(N/2)) end
startwalltime = time()
j = 0
R = 1.0
cyc = 0
Eo = energies(prob) # Reference upper layer KE level for renormalization.
ke1o, ke2o, peo = Eo[1][1], Eo[1][2], Eo[2][1]
t = Array([clock.t])
KE1 = Array([ke1o])
KE2 = Array([ke2o])
PE = Array([peo])

while cyc<cycles
  stepforward!(prob)
  updatevars!(prob)
  E = energies(prob)
  ke1, ke2, pe = E[1][1], E[1][2], E[2][1]

  if j % nlog == 0
    cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy])
    log = @sprintf("step: %05d, t: %2.1f, R: %1.2f, cfl: %1.2f, KE₁/KE₂: %1.2f, walltime: %1.2f min", clock.step, clock.t, R, cfl, ke1/ke2, (time()-startwalltime)/60)
    println(log)

    if cyc==cyclesm && PSI_HOVMOLLER
      if isnothing(thovm)
        ψ₁hovm = vars.ψ[:, fm, 1]
        ψ₂hovm = vars.ψ[:, fm, 2]
        thovm = Array([clock.t])
      else
        ψ₁hovm = cat(ψ₁hovm, vars.ψ[:, fm, 1], dims=2)
        ψ₂hovm = cat(ψ₂hovm, vars.ψ[:, fm, 2], dims=2)
        append!(thovm, clock.t)
      end
    end
  end

  append!(t, clock.t)
  append!(KE1, ke1)
  append!(KE2, ke2)
  append!(PE, pe)

  j += 1
  R = ke1o/ke1 # renormalize (KE_{ref}/KE).
  if R<Rthresh
    set_q!(prob, vars.q*R)
    cyc += 1
    println("")
    println("Completed renormalization cycle", " ", cyc, "/", cycles)
    println("")
  end
end

# Calculate bulk metrics after final renormalization.
E, F = energies(prob), fluxes(prob)
ke1, ke2, pe = E[1][1], E[1][2], E[2][1]
keratio = ke1/ke2
mom1, mom2, thick = F[1][1], F[1][2], F[2][1]

# Get growth rate from exponential fit to upper-layer KE time series.
σ = get_σ(t, KE1, KE2, PE)

if SAVE_OUTPUT
  println("Saving output data.")
  tk = Int(kt/π)
  Nstr = string(N)
  F1str = string(Int(F₁))
  htstr = string(Int(ht))

  if ttype=="rand"
      npzname = "lin_N" * Nstr * "_ht" * htstr * "_F1" * F1str * "_" * split(fname_hrand, "/")[end][1:end-4]
  elseif ttype=="slop"
      npzname = "lin_N" * Nstr * "_al" * string(α) * "_F1" * F1str * "_" * ttype
  else
      npzname = "lin_N" * Nstr * "_ht" * htstr * "_F1" * F1str * "_" * ttype * string(tk)
  end

  npzname = "../simulations/" * npzname * ".npz"
  ψ₁, ψ₂ = vars.ψ[:, :, 1], vars.ψ[:, :, 2]
  npzdata = Dict("p1" => ψ₁', "p2" => ψ₂', "thick" => thick, "mom1" => mom1, "mom2" => mom2, "ke1" => ke1, "ke2" => ke2, "pe" => pe, "N" => N, "L" => L, "dt" => dt, "F1" => F₁, "U2fac" => U₂fac, "al" => α, "bet" => β, "d12" => δ, "ht" => ht, "tk" => tk, "sigma" => σ, "keratio" => keratio)

  if PSI_HOVMOLLER
    npzdata = merge(npzdata, Dict("p1hovm" => ψ₁hovm', "p2hovm" => ψ₂hovm', "thovm" => thovm))
  end

  npzwrite(npzname, npzdata)
end

end # all_F₁s, all_hts

max_stars_repo_name: apaloczy/RoughOceanQGEnergetics

max_issues_repo_name: apaloczy/RoughOceanQGEnergetics

max_forks_repo_name: apaloczy/RoughOceanQGEnergetics


text:
using LocalApproximationPolicyEvaluation
using Test
using POMDPModels
using POMDPs
using GridInterpolations
using LocalFunctionApproximation

# Step 1 - Setup the gridworld problem
g_size = (9,9)
g = SimpleGridWorld(size = g_size, rewards = Dict(GWPos(g_size...) => 1, GWPos(1,1) => 0), tprob = 1., discount=1)
action_probability(g::SimpleGridWorld, s, a) = 0.25
POMDPs.convert_s(::Type{AbstractArray}, s::GWPos, g::SimpleGridWorld) = SVector{4, Float64}(1., s[1], s[2], s[1]*s[2])
POMDPs.gen(mdp::SimpleGridWorld, s, a, rng = Random.GLOBAL_RNG) = (sp =rand(transition(g, s, a )), r=reward(g, s, a))
POMDPs.initialstate(g::SimpleGridWorld) = initialstate_distribution(g)

# Step 2 - Solve the problem semi-exactly using local approximation
grid = RectangleGrid([1:g_size[1] ...], [1:g_size[2]...])
interp = LocalGIFunctionApproximator(grid)
solver = LocalPolicyEvalSolver(interp, action_probability, is_mdp_generative = true, n_generative_samples = 1, max_iterations = 2000, belres = 1e-6)
policy = solve(solver, g)
values = [value(policy, s) for s in states(g)]
@test value(policy, GWPos(9,9)) == 1.0
@test value(policy, GWPos(8,9)) > .25
@test value(policy, GWPos(1,1)) == 0

ext: jl


text:
"""
    nfan(n_out, n_in=1) -> Tuple
    nfan(dims...)
    nfan(dims::Tuple)

For a layer characterized by dimensions `dims`, return a tuple `(fan_in, fan_out)`, where `fan_in`
is the number of input neurons connected to an output one, and `fan_out` is the number of output neurons
connected to an input one.

This function is mainly used by weight initializers, e.g., [`kaiming_normal`](@ref Flux.kaiming_normal).

# Examples

```jldoctest
julia> layer = Dense(10, 20);

julia> Flux.nfan(size(layer.weight))
(10, 20)

julia> layer = Conv((3, 3), 2=>10);

julia> Flux.nfan(size(layer.weight))
(18, 90)
```
"""
nfan() = 1, 1 # fan_in, fan_out
nfan(n) = 1, n # A vector is treated as a n×1 matrix
nfan(n_out, n_in) = n_in, n_out # In case of Dense kernels: arranged as matrices
nfan(dims::Tuple) = nfan(dims...)
nfan(dims...) = prod(dims[1:end-2]) .* (dims[end-1], dims[end]) # In case of convolution kernels

ofeltype(x, y) = convert(float(eltype(x)), y)
epseltype(x) = eps(float(eltype(x)))

"""
    rng_from_array([x])

Create an instance of the RNG most appropriate for `x`.
The current defaults are:
- `x isa CuArray`: `CUDA.default_rng()`, else:
- `x isa AbstractArray`, or no `x` provided:
  - Julia version is < 1.7: `Random.GLOBAL_RNG`
  - Julia version is >= 1.7: `Random.default_rng()`
"""
rng_from_array(::AbstractArray) = rng_from_array()
rng_from_array(::CuArray) = CUDA.default_rng()
if VERSION >= v"1.7"
  rng_from_array() = Random.default_rng()
else
  rng_from_array() = Random.GLOBAL_RNG
end

"""
    glorot_uniform([rng=GLOBAL_RNG], size...; gain = 1) -> Array
    glorot_uniform([rng]; kw...) -> Function

Return an `Array{Float32}` of the given `size` containing random numbers drawn from a uniform
distribution on the interval ``[-x, x]``, where `x = gain * sqrt(6 / (fan_in + fan_out))`.

This method is described in [1] and also known as Xavier initialization.

# Examples
```jldoctest; setup = :(using Random; Random.seed!(0))
julia> Flux.glorot_uniform(3, 4) |> summary
"3×4 Matrix{Float32}"

julia> round.(extrema(Flux.glorot_uniform(10, 100)), digits=3)
(-0.232f0, 0.234f0)

julia> round.(extrema(Flux.glorot_uniform(100, 10)), digits=3)
(-0.233f0, 0.233f0)

julia> round.(extrema(Flux.glorot_uniform(100, 100)), digits=3)
(-0.173f0, 0.173f0)

julia> Dense(3 => 2, tanh; init = Flux.glorot_uniform(MersenneTwister(1)))
Dense(3 => 2, tanh)  # 8 parameters

julia> ans.bias
2-element Vector{Float32}:
 0.0
 0.0
```

# References

[1] Glorot, Xavier, and Yoshua Bengio. "Understanding the difficulty of training deep feedforward neural networks." _Proceedings of the thirteenth international conference on artificial intelligence and statistics_. 2010.
"""
function glorot_uniform(rng::AbstractRNG, dims::Integer...; gain::Real=1)
  scale = Float32(gain) * sqrt(24.0f0 / sum(nfan(dims...)))
  (rand(rng, Float32, dims...) .- 0.5f0) .* scale
end
glorot_uniform(dims::Integer...; kw...) = glorot_uniform(rng_from_array(), dims...; kw...)
glorot_uniform(rng::AbstractRNG=rng_from_array(); init_kwargs...) = (dims...; kwargs...) -> glorot_uniform(rng, dims...; init_kwargs..., kwargs...)

ChainRulesCore.@non_differentiable glorot_uniform(::Any...)

"""
    glorot_normal([rng=GLOBAL_RNG], size...; gain = 1) -> Array
    glorot_normal([rng]; kw...) -> Function

Return an `Array{Float32}` of the given `size` containing random numbers drawn from a normal
distribution with standard deviation `gain * sqrt(2 / (fan_in + fan_out))`,
using [`nfan`](@ref Flux.nfan).

This method is described in [1] and also known as Xavier initialization.

# Examples
```jldoctest; setup = :(using Random; Random.seed!(0))
julia> using Statistics

julia> round(std(Flux.glorot_normal(10, 1000)), digits=3)
0.044f0

julia> round(std(Flux.glorot_normal(1000, 10)), digits=3)
0.044f0

julia> round(std(Flux.glorot_normal(1000, 1000)), digits=3)
0.032f0

julia> Dense(10 => 1000, tanh; init = Flux.glorot_normal(gain=100))
Dense(10 => 1000, tanh)  # 11_000 parameters

julia> round(std(ans.weight), sigdigits=3)
4.45f0
```

# References

[1] Glorot, Xavier, and Yoshua Bengio. "Understanding the difficulty of training deep feedforward neural networks." _Proceedings of the thirteenth international conference on artificial intelligence and statistics_. 2010.
"""
function glorot_normal(rng::AbstractRNG, dims::Integer...; gain::Real=1)
  std = Float32(gain) * sqrt(2.0f0 / sum(nfan(dims...)))
  randn(rng, Float32, dims...) .* std
end
glorot_normal(dims::Integer...; kwargs...) = glorot_normal(rng_from_array(), dims...; kwargs...)
glorot_normal(rng::AbstractRNG=rng_from_array(); init_kwargs...) = (dims...; kwargs...) -> glorot_normal(rng, dims...; init_kwargs..., kwargs...)

ChainRulesCore.@non_differentiable glorot_normal(::Any...)

"""
    kaiming_uniform([rng=GLOBAL_RNG], size...; gain = √2) -> Array
    kaiming_uniform([rng]; kw...) -> Function

Return an `Array{Float32}` of the given `size` containing random numbers drawn from a uniform distribution
on the interval `[-x, x]`, where `x = gain * sqrt(3/fan_in)` using [`nfan`](@ref Flux.nfan).

This method is described in [1] and also known as He initialization.

# Examples
```jldoctest; setup = :(using Random; Random.seed!(0))
julia> round.(extrema(Flux.kaiming_uniform(100, 10)), digits=3)
(-0.774f0, 0.774f0)

julia> round.(extrema(Flux.kaiming_uniform(10, 100)), digits=3)
(-0.245f0, 0.244f0)

julia> round.(extrema(Flux.kaiming_uniform(100, 100)), digits=3)
(-0.245f0, 0.245f0)
```

# References

[1] He, Kaiming, et al. "Delving deep into rectifiers: Surpassing human-level performance on imagenet classification." _Proceedings of the IEEE international conference on computer vision_. 2015.
"""
function kaiming_uniform(rng::AbstractRNG, dims::Integer...; gain::Real = √2)
  bound = Float32(√3 * gain / sqrt(first(nfan(dims...)))) # fan_in
  return (rand(rng, Float32, dims...) .- 0.5f0) .* 2bound
end

kaiming_uniform(dims::Integer...; kwargs...) = kaiming_uniform(rng_from_array(), dims...; kwargs...)
kaiming_uniform(rng::AbstractRNG=rng_from_array(); init_kwargs...) = (dims...; kwargs...) -> kaiming_uniform(rng, dims...; init_kwargs..., kwargs...)

ChainRulesCore.@non_differentiable kaiming_uniform(::Any...)

"""
    kaiming_normal([rng=GLOBAL_RNG], size...; gain = √2) -> Array
    kaiming_normal([rng]; kw...) -> Function

Return an `Array{Float32}` of the given `size` containing random numbers taken from a normal
distribution standard deviation `gain / sqrt(fan_in)`, using [`nfan`](@ref Flux.nfan).

This method is described in [1] and also known as He initialization.

# Examples
```jldoctest; setup = :(using Random; Random.seed!(0))
julia> using Statistics

julia> round(std(Flux.kaiming_normal(10, 1000)), digits=3)
0.045f0

julia> round(std(Flux.kaiming_normal(1000, 10)), digits=3)
0.447f0

julia> round(std(Flux.kaiming_normal(1000, 1000)), digits=3)
0.045f0
```

# References

[1] He, Kaiming, et al. "Delving deep into rectifiers: Surpassing human-level performance on imagenet classification." _Proceedings of the IEEE international conference on computer vision_. 2015.
"""
function kaiming_normal(rng::AbstractRNG, dims::Integer...; gain::Real = √2f0)
  std = Float32(gain / sqrt(first(nfan(dims...)))) # fan_in
  return randn(rng, Float32, dims...) .* std
end

kaiming_normal(dims::Integer...; kwargs...) = kaiming_normal(rng_from_array(), dims...; kwargs...)
kaiming_normal(rng::AbstractRNG; init_kwargs...) = (dims...; kwargs...) -> kaiming_normal(rng, dims...; init_kwargs..., kwargs...)

ChainRulesCore.@non_differentiable kaiming_normal(::Any...)

"""
    truncated_normal([rng=GLOBAL_RNG], size...; mean = 0, std = 1, lo = -2, hi = 2) -> Array
    truncated_normal([rng]; kw...) -> Function

Return an `Array{Float32}` of the given `size` where each element is drawn from a truncated normal distribution.
The numbers are distributed like `filter(x -> lo<=x<=hi, mean .+ std .* randn(100))`.

The values are generated by sampling a Uniform(0, 1) (`rand()`) and then
applying the inverse CDF of the truncated normal distribution.
This method works best when `lo ≤ mean ≤ hi`.

# Examples
```jldoctest
julia> using Statistics

julia> Flux.truncated_normal(3, 4) |> summary
"3×4 Matrix{Float32}"

julia> round.(extrema(Flux.truncated_normal(10^6)); digits=3)
(-2.0f0, 2.0f0)

julia> round(std(Flux.truncated_normal(10^6; lo = -100, hi = 100)))
1.0f0
```
"""
function truncated_normal(rng::AbstractRNG, dims::Integer...; mean = 0, std = 1, lo = -2, hi = 2)
  norm_cdf(x) = 0.5 * (1 + erf(x/√2))
  if (mean < lo - 2 * std) || (mean > hi + 2 * std)
    @warn "Mean is more than 2 std outside the limits in truncated_normal, so the distribution of values may be inaccurate." maxlog=1
  end
  l = norm_cdf((lo - mean) / std)
  u = norm_cdf((hi - mean) / std)
  xs = rand(rng, Float32, dims...)
  broadcast!(xs, xs) do x
    x = x * 2(u - l) + (2l - 1)
    x = erfinv(x)
    x = clamp(x * std * √2 + mean, lo, hi)
  end
  return xs
end

truncated_normal(dims::Integer...; kwargs...) = truncated_normal(rng_from_array(), dims...; kwargs...)
truncated_normal(rng::AbstractRNG=rng_from_array(); init_kwargs...) = (dims...; kwargs...) -> truncated_normal(rng, dims...; init_kwargs..., kwargs...)

ChainRulesCore.@non_differentiable truncated_normal(::Any...)

"""
    orthogonal([rng=GLOBAL_RNG], size...; gain = 1) -> Array
    orthogonal([rng]; kw...) -> Function

Return an `Array{Float32}` of the given `size` which is a (semi) orthogonal matrix, as described in [1].

Cannot construct a vector, i.e. `length(size) == 1` is forbidden.
For `length(size) > 2`, a `prod(size[1:(end - 1)])` by `size[end]` orthogonal matrix
is computed before reshaping it to the original dimensions.

# Examples
```jldoctest; setup = :(using LinearAlgebra)
julia> W = Flux.orthogonal(5, 7);

julia> summary(W)
"5×7 Matrix{Float32}"

julia> W * W' ≈ I(5)
true

julia> W2 = Flux.orthogonal(7, 5);

julia> W2 * W2' ≈ I(7)
false

julia> W2' * W2 ≈ I(5)
true

julia> W3 = Flux.orthogonal(3, 3, 2, 4);

julia> transpose(reshape(W3, :, 4)) * reshape(W3, :, 4) ≈ I(4)
true
```

# References

[1] Saxe, McClelland, Ganguli. "Exact solutions to the nonlinear dynamics of learning in deep linear neural networks", ICLR 2014, https://arxiv.org/abs/1312.6120

"""
function orthogonal(rng::AbstractRNG, rows::Integer, cols::Integer; gain::Real = 1)
  if rows < cols
    return permutedims(orthogonal(rng, cols, rows; gain))
  end
  mat = randn(rng, Float32, rows, cols)
  Q, R = LinearAlgebra.qr(mat)
  mat .= Array(Q) * sign.(LinearAlgebra.Diagonal(R)) .* Float32(gain)
  return mat
end

function orthogonal(rng::AbstractRNG, d1::Integer, ds::Integer...; kwargs...)
  dims = (d1, ds...)
  rows = prod(dims[1:end-1])
  cols = dims[end]
  return reshape(orthogonal(rng, rows, cols; kwargs...), dims)
end

orthogonal(dims::Integer...; kwargs...) = orthogonal(rng_from_array(), dims...; kwargs...)
orthogonal(rng::AbstractRNG=rng_from_array(); init_kwargs...) = (dims::Integer...; kwargs...) -> orthogonal(rng, dims...; init_kwargs..., kwargs...)

ChainRulesCore.@non_differentiable orthogonal(::Any...)

"""
    sparse_init([rng=GLOBAL_RNG], rows, cols; sparsity, std = 0.01) -> Array
    sparse_init([rng]; kw...) -> Function

Return a `Matrix{Float32}` of size `rows, cols` where each column contains a fixed fraction of
zero elements given by `sparsity`. Non-zero elements are normally distributed
with a mean of zero and standard deviation `std`.

This method is described in [1].

# Examples
```jldoctest; setup = :(using Random; Random.seed!(0))
julia> count(iszero, Flux.sparse_init(10, 10, sparsity=1/5))
20

julia> sum(0 .== Flux.sparse_init(10, 11, sparsity=0.9), dims=1)
1×11 Matrix{Int64}:
 9  9  9  9  9  9  9  9  9  9  9

julia> Dense(3 => 10, tanh; init=Flux.sparse_init(sparsity=0.5))
Dense(3 => 10, tanh)  # 40 parameters

julia> count(iszero, ans.weight, dims=1)
1×3 Matrix{Int64}:
 5  5  5
```

# References

[1] Martens, J, "Deep learning via Hessian-free optimization" _Proceedings of the 27th International Conference on International Conference on Machine Learning_. 2010.
"""
function sparse_init(rng::AbstractRNG, dims::Integer...; sparsity, std = 0.01)
  if length(dims) != 2
    throw(ArgumentError("Only 2-dimensional outputs are supported for sparse initialization."))
  end
  rows, cols = dims
  prop_zero = min(1.0, sparsity)
  num_zeros = ceil(Integer, prop_zero * rows)
  sparse_array = randn(rng, Float32, dims...) .* Float32(std)
  sparse_array[1:num_zeros, :] .= 0f0
  return mapslices(shuffle, sparse_array, dims=1)
end

sparse_init(dims::Integer...; kwargs...) = sparse_init(rng_from_array(), dims...; kwargs...)
sparse_init(rng::AbstractRNG=rng_from_array(); init_kwargs...) = (dims...; kwargs...) -> sparse_init(rng, dims...; init_kwargs..., kwargs...)

ChainRulesCore.@non_differentiable sparse_init(::Any...)

"""
    identity_init(size...; gain=1, shift=0) -> Array
    identity_init(; kw...) -> Function

Return an `Array{Float32}` of the given `size` which yields an identity mapping when used as parameters in
most Flux layers. Use `gain` to scale the identity by a constant.

Often useful in the context of transfer learning, i.e when one wants to add more capacity to
a model but start from the same mapping.

Has the following behaviour
*  1D: A `Vector` of `zeros` (useful for an identity bias)
*  2D: An identity matrix (useful for an identity matrix multiplication)
*  More than 2D: A dense block array of center tap spatial filters (useful for an identity convolution)

Some caveats:
* Not all layers will be identity mapping when used with this init. Exceptions
  include recurrent layers and normalization layers.

* Layers must have `input_size == output_size` for identity mapping to be
  possible. When this is not the case, extra dimensions of the array are padded with zeros.

* For convolutional layers, in addition to the above, the kernel sizes must also be odd and
  padding must be applied so that output feature maps have the same size as input feature maps,
  e.g by using [`SamePad`](@ref).

Use keyword `shift` (integer or tuple) to apply circular shift to the output,
equivalent to `Base.circshift(identity_init(size...), shift)`.

For consistency with other initialisers, it accepts `rng::AbstractRNG` as an optional
first argument. But this is ignored, since the result is not random.

# Examples
```jldoctest
julia> Flux.identity_init(3,5)
3×5 Matrix{Float32}:
 1.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0

julia> Dense(5 => 3, relu, init=Flux.identity_init)([1,-2,3,-4,5])
3-element Vector{Float32}:
 1.0
 0.0
 3.0

julia> Flux.identity_init(3,3,2; gain=100)
3×3×2 Array{Float32, 3}:
[:, :, 1] =
   0.0  0.0  0.0
 100.0  0.0  0.0
   0.0  0.0  0.0

[:, :, 2] =
 0.0    0.0  0.0
 0.0  100.0  0.0
 0.0    0.0  0.0

julia> x4 = cat([1 2 3; 4 5 6; 7 8 9]; dims=4);

julia> Conv((2,2), 1 => 1, init=Flux.identity_init(gain=10), pad=SamePad())(x4)
3×3×1×1 Array{Float32, 4}:
[:, :, 1, 1] =
 10.0  20.0  30.0
 40.0  50.0  60.0
 70.0  80.0  90.0
```
"""
identity_init(cols::Integer; gain::Real=1, shift=0) = zeros32(cols) # Assume bias

# Assume matrix multiplication
identity_init(rows::Integer, cols::Integer; gain::Real=1, shift=0) = circshift(Matrix{Float32}(I * gain, rows,cols), shift)

# Assume convolution
function identity_init(dims::Integer...; gain::Real=1, shift=0)
  nin, nout = dims[end-1], dims[end]
  centers = map(d -> cld(d, 2), dims[1:end-2])
  weights = zeros32(dims...)
  for i in 1:min(nin,nout)
    weights[centers..., i, i] = gain
  end
  return circshift(weights, shift)
end

# For consistency, it accepts an RNG, but ignores it:
identity_init(::AbstractRNG, dims::Integer...; kwargs...) = identity_init(dims...; kwargs...)
identity_init(rng::AbstractRNG=rng_from_array(); init_kwargs...) = (args...;kwargs...) -> identity_init(rng, args...; init_kwargs..., kwargs...)

ChainRulesCore.@non_differentiable identity_init(::Any...)

ones32(dims::Integer...) = Base.ones(Float32, dims...)
zeros32(dims::Integer...) = Base.zeros(Float32, dims...)

"""
    ones32(size...) = ones(Float32, size...)
    zeros32(size...) = zeros(Float32, size...)

Return an `Array{Float32}` of the given `size`.
"""
ones32(dims...) = Base.ones(Float32, dims...)

@doc @doc(ones32)
zeros32(dims...) = Base.zeros(Float32, dims...)

"""
    rand32([rng], size...)
    randn32([rng], size...)

Return an `Array{Float32}` of the given `size`, filled like `rand` or `randn`.
When the size is not provided, `rand32(rng::AbstractRNG)` returns a function.
"""
rand32(dims::Integer...) = Base.rand(Float32, dims...)
rand32(rng::AbstractRNG, dims::Integer...) = Base.rand(rng, Float32, dims...)
rand32(rng::AbstractRNG) = (dims...,) -> Base.rand(rng, Float32, dims...)

@doc @doc(rand32)
randn32(dims::Integer...) = Base.randn(Float32, dims...)
randn32(rng::AbstractRNG, dims::Integer...) = Base.randn(rng, Float32, dims...)
randn32(rng::AbstractRNG) = (dims...,) -> Base.randn(rng, Float32, dims...)

"""
    create_bias(weights, bias, size...)

Return a bias parameter for a layer, based on the value given
to the constructor's keyword `bias=bias`.

* `bias == true` creates a trainable array of the given size, of the same type as `weights`, initialised to zero.
* `bias == false` returns `false`, which is understood by AD to be non-differentiable.
* `bias::AbstractArray` uses the array provided, provided it has the correct size.
  It does not at present correct the `eltype` to match that of `weights`.
"""
function create_bias(weights::AbstractArray, bias::Bool, dims::Integer...)
  bias ? fill!(similar(weights, dims...), 0) : false
end
function create_bias(weights::AbstractArray, bias::AbstractArray, dims::Integer...)
  size(bias) == dims || throw(DimensionMismatch("expected bias of size $(dims), got size $(size(bias))"))
  bias
end

# Other

"""
    throttle(f, timeout; leading=true, trailing=false)

Return a function that when invoked, will only be triggered at most once
during `timeout` seconds.

Normally, the throttled function will run as much as it can, without ever
going more than once per `wait` duration; but if you'd like to disable the
execution on the leading edge, pass `leading=false`. To enable execution on
the trailing edge, pass `trailing=true`.
"""
function throttle(f, timeout; leading=true, trailing=false)
  cooldown = true
  later = nothing
  result = nothing

  function throttled(args...; kwargs...)
    yield()

    if cooldown
      if leading
        result = f(args...; kwargs...)
      else
        later = () -> f(args...; kwargs...)
      end

      cooldown = false
      @async try
        while (sleep(timeout); later != nothing)
          later()
          later = nothing
        end
      finally
        cooldown = true
      end
    elseif trailing
      later = () -> (result = f(args...; kwargs...))
    end

    return result
  end
end

"""
    modules(m)

Return an iterator over non-leaf objects
that can be reached by recursing `m` over
the children given by [`functor`](@ref).

Useful for applying a function (e.g. a regularizer)
over specific modules or subsets of the parameters
(e.g. the weights but not the biases).

# Examples

```jldoctest
julia> m1 = Chain(Dense(28^2, 64), BatchNorm(64, relu));

julia> m2 = Chain(m1, Dense(64, 10))
Chain(
  Chain(
    Dense(784 => 64),                   # 50_240 parameters
    BatchNorm(64, relu),                # 128 parameters, plus 128
  ),
  Dense(64 => 10),                      # 650 parameters
)         # Total: 6 trainable arrays, 51_018 parameters,
          # plus 2 non-trainable, 128 parameters, summarysize 200.312 KiB.

julia> Flux.modules(m2)
7-element Vector{Any}:
 Chain(Chain(Dense(784 => 64), BatchNorm(64, relu)), Dense(64 => 10))  # 51_018 parameters, plus 128 non-trainable
 (Chain(Dense(784 => 64), BatchNorm(64, relu)), Dense(64 => 10))
 Chain(Dense(784 => 64), BatchNorm(64, relu))  # 50_368 parameters, plus 128 non-trainable
 (Dense(784 => 64), BatchNorm(64, relu))
 Dense(784 => 64)    # 50_240 parameters
 BatchNorm(64, relu)  # 128 parameters, plus 128 non-trainable
 Dense(64 => 10)     # 650 parameters

julia> L2(m) = sum(sum(abs2, l.weight) for l in Flux.modules(m) if l isa Dense)
L2 (generic function with 1 method)

julia> L2(m2) isa Float32
true
```
"""
modules(m) = [x for x in Functors.fcollect(m) if !isleaflike(x)]

@nograd modules # TODO: is this correct? might fail with explicit parameters.
function ChainRulesCore.rrule(::typeof(modules), m)
  modules(m), dm -> error("Flux.modules is not at present differentiable, sorry")
end

isleaflike(x) = Functors.isleaf(x)
isleaflike(::Tuple{Vararg{<:Number}}) = true
isleaflike(::Tuple{Vararg{<:AbstractArray{<:Number}}}) = true

"""
    patience(predicate, wait)

Return a function that internally counts by one when
`predicate(...) == true`, otherwise the count is reset to zero.
If the count is greater than or equal to `wait`,
the function returns `true`, otherwise it returns `false`.

# Examples
```jldoctest
julia> loss() = rand();

julia> trigger = Flux.patience(() -> loss() < 1, 3);

julia> Flux.@epochs 10 begin
         trigger() && break
       end
[ Info: Epoch 1
[ Info: Epoch 2
[ Info: Epoch 3
```
"""
function patience(predicate, wait)
  let count = 0
    function on_trigger(args...; kwargs...)
      count = predicate(args...; kwargs...) ? count + 1 : 0

      return count >= wait
    end
  end
end

"""
    early_stopping(f, delay; distance = -, init_score = 0, min_dist = 0)

Return a function that internally counts by one when
`distance(best_score, f(...)) <= min_dist`, where
`best_score` is the last seen best value of `f(...)`.
If the count is greater than or equal to `delay`,
the function returns `true`, otherwise it returns `false`.
The count is reset when `distance(best_score, f(...)) > min_dist`.

# Examples
```jldoctest
julia> loss = let l = 0
         () -> l += 1
       end; # pseudo loss function that returns increasing values

julia> es = Flux.early_stopping(loss, 3);

julia> Flux.@epochs 10 begin
         es() && break
       end
[ Info: Epoch 1
[ Info: Epoch 2
[ Info: Epoch 3
```
"""
function early_stopping(f, delay; distance = -, init_score = 0, min_dist = 0)
  trigger = let best_score = init_score
    (args...; kwargs...) -> begin
      score = f(args...; kwargs...)
      Δ = distance(best_score, score)
      best_score = Δ < 0 ? best_score : score

      return Δ < min_dist
    end
  end

  return patience(trigger, delay)
end

"""
    plateau(f, width; distance = -, init_score = 0, min_dist = 1f-6)

Return a function that internally counts by one when
`abs(distance(last_score, f(...))) <= min_dist`, where
`last_score` holds the last value of `f(...)`.
If the count is greater than or equal to `width`,
the function returns `true`, otherwise it returns `false`.
The count is reset when `abs(distance(last_score, f(...))) > min_dist`.

# Examples
```jldoctest
julia> f = let v = 10
         () -> v = v / abs(v) - v
       end; # -9, 8, -7, 6, ...

julia> trigger = Flux.plateau(f, 3; init_score=10, min_dist=18);

julia> Flux.@epochs 10 begin
         trigger() && break
       end
[ Info: Epoch 1
[ Info: Epoch 2
[ Info: Epoch 3
[ Info: Epoch 4
```
"""
function plateau(f, width; distance = -, init_score = 0, min_dist = 1f-6)
  is_plateau = let last_score = init_score
    (args...; kwargs...) -> begin
      score = f(args...; kwargs...)
      Δ = abs(distance(last_score, score))
      last_score = score

      return Δ < min_dist
    end
  end

  return patience(is_plateau, width)
end


text:
export add_meta!, locations, locations!, loci, genotypes, get_genotypes, get_genotype, populations, population, populations!, population!, exclude, remove, omit, exclude!, remove!, omit!, samples

"""
    add_meta!(popdata::PopData, metadata::T; name::String, loci::Bool = true, categorical::Bool = true) where T <: AbstractVector
Add an additional metadata information to a `PopData` object. Mutates `PopData` in place. Metadata
must be in the same order as the samples in `PopData.meta`.

#### Arguments
- `popdata` : The `PopData` object to add information to
- `metadata` : A `Vector` with the metadata you wish to add to the `PopData`, in the same order as the names appear in `PopData.meta`

#### Keyword Arguments
- `name` : String of the name of this new column
- `loci` : Boolean of whether to also add this information to `PopData.loci` (default: `true`)
- `categorical` : Boolean of whether the metadata being added is categorical aka "factors" (default: `true`)
"""
function add_meta!(popdata::PopData, metadata::T; name::String, loci::Bool = true, categorical::Bool = true) where T <: AbstractVector
    length(metadata) != length(popdata.meta.name) && error("Provided metadata vector (n = $length(metadata)) and samples in PopData (n = $length(popdata.meta.name)) have different lengths")
    @info "Adding $Symbol(name) column to .meta dataframe"
    # add to meta
    insertcols!(popdata.meta, Symbol(name) => metadata)

    # add to loci
    if loci
        @info "Adding $Symbol(name) column to .meta and .loci dataframes"
        tmp = DataFrame(:name => popdata.meta.name, Symbol(name) => metadata)
        popdata.loci = outerjoin(popdata.loci, tmp, on = :name)
        if categorical == true
            popdata.loci[name] = PooledArray(popdata.loci[name])
        end
    end
    return
end

"""
    add_meta!(popdata::PopData, samples::Vector{String}, metadata::T; name::String, loci::Bool = true, categorical::Bool = true) where T <: AbstractVector
Add an additional metadata information to a `PopData` object. Mutates `PopData` in place. Takes a vector of
sample names if the metadata is not in the same order as samples appear in `PopData.meta`.

#### Arguments
- `popdata` : The `PopData` object to add information to
- `sample` : A `Vector{String}` of sample names corresponding to the order of the `metadata`
- `metadata` : A `Vector` with the metadata you wish to add to the `PopData`, in the same order as the names appear in `PopData.meta`

#### Keyword Arguments
- `name` : String of the name of this new column
- `loci` : Boolean of whether to also add this information to `PopData.loci` (default: `true`)
- `categorical` : Boolean of whether the metadata being added is categorical aka "factors" (default: `true`)
"""
function add_meta!(popdata::PopData, samples::Vector{String}, metadata::T; name::String, loci::Bool = true) where T <: AbstractVector
    length(samples) != length(popdata.meta.name) && error("Provided sample vector (n = $length(samples)) and samples in PopData (n = $length(popdata.meta.name)) have different lengths")
    length(samples) != length(metadata) && error("Sample names (n = $length(samples)) and metadata vectors (n = $length(metadata)) have different lengths")
    sort(samples) != sort(popdata.meta.name) && error("Sample names are not identical")
    @info "Adding $Symbol(name) column to .meta dataframe"

    tmp = DataFrame(:name => samples, Symbol(name) => metadata)
    popdata.meta = outerjoin(popdata.meta, tmp, on = :name)

    if loci
        @info "Adding $Symbol(name) column to .meta and .loci dataframes"
        popdata.loci = outerjoin(popdata.loci, tmp, on = :name)
        if categorical == true
            popdata.loci[name] = PooledArray(popdata.loci[name])
        end
    end
    return
end

"""
    locations(data::PopData)
View the longitude and latitude data in a `PopData` object. Returns a table
derived from the PopData. Changes made to this table will not alter the source
`PopData` object.

Use `locations!` to add spatial data to a `PopData` object.
"""
function locations(data::PopData)
    @view data.meta[!, [:longitude, :latitude]]
end

"""
    locations!(data::PopData; long::Vector{Float64}, lat::Vector{Float64})
Replaces existing `PopData` location data (longitude `long`, latitude `lat`).
Takes **decimal degrees** as a `Vector` of any `AbstractFloat`.
## Formatting requirements
- Decimal Degrees format: `-11.431`
- **Must** use negative sign `-` instead of cardinal directions
- Location data must be in the order that samples appear in your `PopData`
- Missing data should be coded as `missing` values of type `Missing` (can be accomplished with `replace!()`)
### Example
```
ncats = @nancycats ;
x = rand(237) ; y = rand(237)
locations!(ncats, long = x, lat = y)
```
"""
function locations!(data::PopData, long::Vector{Union{Missing,T}}, lat::Vector{Union{Missing,T}}) where T <: AbstractFloat
    long_len = length(long)
    lat_len = length(lat)
    long_len != lat_len && error("latitude ($lat_len) and longitude ($long_len) arrays not equal in length")
    long_len != length(data.meta.name) && error("lat/long array length ($long_len) and number of samples in PopData ($long_len) are not equal")

    data.meta.longitude .= long
    data.meta.latitude .= lat
    return
end

function locations!(data::PopData, long::Vector{T}, lat::Vector{T}) where T <: AbstractFloat
    # convert to the right type and use locations!()
    lat_adjust = lat |> Vector{Union{Missing, Float32}}
    long_adjust = long |> Vector{Union{Missing, Float32}}
    locations!(data, long = long_adjust, lat = lat_adjust)
end

"""
    locations!(data::PopData; long::Vector{String}, lat::Vector{String})
Replaces existing `PopData` location data (longitude `long`, latitude `lat`). Takes
**decimal minutes** or **degrees minutes seconds** format as a `Vector` of `String`. Recommended to use `CSV.read`
from `CSV.jl` to import your spatial coordinates from a text file.
## Formatting requirements
- Coordinates as a `String` separated by spaces (`"11 43 41"`) or colons (`"11:43:41"`)
- Must use negative sign (`"-11 43.52"`) or single-letter cardinal direction (`"11 43.52W"`)
- Missing data should be coded as the string `"missing"` (can be accomplished with `replace!()`)
- Can mix colons and spaces (although it's bad practice)
### NOTE
If you read in the coordinate data as 4 vectors (longitude degrees, longitude minutes, latitude degrees, latitude minutes),
then the easiest course of action would be to merge them into two vectors of strings
(one for longitude, one for latitude):
```
long_string = string.(lat_deg, " ", lat_min)
lat_string = string.(long_deg, " ", long_min)
```
and use these as inputs into `locations!`

### Example
```
ncats = @nancycats;
x = fill("11 22.33W", 237) ; y = fill("-41 31.52", 237)
locations!(ncats, long = x, lat = y)
```
"""
function locations!(data::PopData, long::Vector{String}, lat::Vector{String})
    long_len = length(long)
    lat_len = length(lat)
    lat_len != long_len && error("latitude ($lat_len) and longitude ($long_len) arrays not equal in length")
    lat_len != length(data.meta.name) && error("lat/long array length ($lat_len) and number of samples in PopData ($long_len) are not equal")
    @info "Converting decimal minutes to decimal degrees"
    # convert coordinates to decimal degrees
    data.meta.longitude .= convert_coord.(long)
    data.meta.latitude .= convert_coord.(lat)
    return
end

function locations!(data::PopData; kwargs...)
    kwargs = Dict(kwargs)
    # check for matching lat and long keywords
    if all([haskey(kwargs, :lat), haskey(kwargs, :long)])
        locations!(data, kwargs[:long], kwargs[:lat])
    else
        error("keyword arguments \"lat\" and \"long\" must be supplied together")
    end
end

"""
    loci(data::PopData)
Returns an array of strings of the loci names in a `PopData` object.
"""
function loci(data::PopData)
    unique(data.loci.locus)
end

"""
    get_genotype(data::PopObj; sample::String, locus::String)
Return the genotype of one sample at one locus in a `PopData` object.
### Example
```
cats = @nancycats;
get_genotype(cats, sample = "N115", locus = "fca8")
```
"""
function get_genotype(data::PopData; sample::String, locus::String)
    @views data.loci[(data.loci.name .== sample) .& (data.loci.locus .== locus), :genotype][1]
end

"""
    get_genotypes(data::PopData, sample::String)
Return a vector of all the genotypes of a sample in a `PopData` object. To return a
single genotype at a locus, see `get_genotype`.
```
cats = @nancycats
get_genotypes(cats, "N115")
```
"""
function get_genotypes(data::PopObj, sample::String)
    data.loci[data.loci.name .== sample, :genotype]
end

"""
    get_genotypes(data::PopObj; name::Union{String, Vector{String}}, locus::Union{String, Vector{String}})
Return a table of the genotype(s) of one or more `samples` for one or more
specific `loci` (both as keywords) in a `PopData` object.
### Examples
```
cats = @nancycats;
get_genotypes(cats, name = "N115" , locus = "fca8")
get_genotypes(cats, name = ["N115", "N7"] , locus = "fca8")
get_genotypes(cats, name = "N115" , locus = ["fca8", "fca37"])
get_genotypes(cats, name = ["N1", "N2"] , locus = ["fca8", "fca37"])
```
"""
function get_genotypes(data::PopData; name::Union{String, Vector{String}}, locus::Union{String, Vector{String}})
    if typeof(name) == String
        name = [name]
    end
    if typeof(locus) == String
        locus = [locus]
    end
    @view data.loci[(data.loci.name .∈ Ref(name)) .& (data.loci.locus .∈ Ref(locus)), :]
end

"""
    genotypes(data::PopData, locus::Union{String, Symbol})
Convenience wrapper to return a vector of all the genotypes of a single locus

### Example
```
genotypes(@gulfsharks, "contig_475")
```
"""
function genotypes(data::PopData, locus::String)
    @view data.loci[data.loci.locus .== locus, :genotype]
end

"""
    populations(data::PopData; counts::Bool = false)
View unique population ID's and/or their counts in `PopData`.

- `counts` returns a dataframe of samples per `population` instead (default = `false`)
"""
@inline function populations(data::PopData; counts::Bool = false)
    if all(ismissing.(data.meta.population)) == true
        @info "no population data present in PopData"
        return
    end

    uniq_pops = unique(data.meta.population)

    if counts == false
        return uniq_pops
    else
        pops = countmap(data.meta.population)
        return DataFrame(:population => uniq_pops, :count => [pops[i] for i in uniq_pops])
    end
end

"""
```
# Replace by matching
populations!(data::PopData, rename::Dict)
populations!(data::PopData, rename::Vector{String})
populations!(data::PopData, samples::Vector{String}, populations::Vector{String})
```
Multiple methods to rename or reassign population names to a `PopData`.

## Rename using a Dictionary
Rename existing population ID's of `PopData` using a `Dict` of
`population_name => replacement`
\n**Example**
```
potatopops = Dict("1" => "Idaho", "2" => "Russet")
populations!(potatoes, potatopops)
```

## Rename using a Vector of Strings
`Vector` of new unique population names in the order that they appear in the PopData.meta
\n**Example**
```
potatopops = ["Idaho", "Russet"]
populations!(potatoes, potatopops)
```
## Reassign using samples and new population assignments
Completely reassign populations for each individual. Takes two vectors of strings
as input: one of the sample names, and the other with their new corresponding
population name. This can be useful to change population names for only some individuals.

\n**Example**
```
populations!(potatoes, ["potato_1", "potato_2"], ["north_russet", "south_russet"])
```

"""
function populations!(data::PopData, rename::Dict)
    msg = ""
    @inbounds for key in keys(rename)
        if key ∉ unique(data.meta.population)
            msg *= "  Population \"$key\" not found in PopData\n"
        else
            replace!(data.meta.population, key => rename[key])
            replace!(data.loci.population.pool, key => rename[key])
        end
    end
    msg != "" && printstyled("Warnings:", color = :yellow) ; print("\n"*msg)
    return
end

function populations!(data::PopData, rename::Vector{String})
    current_popnames = unique(data.meta.population)
    rn_dict = Dict{String, String}()
    [rn_dict[string(i)] = j for (i,j) in zip(current_popnames, rename)]
    populations!(data, rn_dict)
    return
end

function populations!(data::PopData, samples::Vector{String}, populations::Vector{String})
    meta_df = groupby(data.meta, :name)
    loci_df = groupby(data.loci, :name)
    for (sample, new_pop) in zip(samples, populations)
        meta_df[(name = sample,)].population .= new_pop
        loci_df[(name = sample,)].population .= new_pop
    end
    # drop old levels
    data.loci.population = data.loci.population |> Array |> PooledArray
    return
end

const population = populations
const population! = populations!
const popnames! = populations!

##### Exclusion #####
"""
    exclude!(data::PopData, kwargs...)
Edit a `PopData` object in-place by excluding all occurences of the specified information.
The keywords can be used in any combination. Synonymous with `omit!` and `remove!`. All
values are converted to `String` for filtering, so `Symbol` and numbers will also work.

### Keyword Arguments
#### `locus`
A `String` or `Vector{String}` of loci you want to remove from the `PopData`.
The keyword `loci` also works.

#### `population`
A `String` or `Vector{String}` of populations you want to remove from the `PopData`
The keyword `populations` also works.

#### `name`
A `String` or `Vector{String}` of samples you want to remove from the `PopData`
The keywords `names`, `sample`, and `samples` also work.

**Examples**
```
cats = @nancycats;
exclude!(cats, name = "N100", population = 1:5)
exclude!(cats, name = ["N100", "N102", "N211"], locus = ["fca8", "fca23"])
exclude!(cats, name = "N102", locus = :fca8, population = "3")
```
"""
function exclude!(data::PopData; population::Any = nothing, locus::Any = nothing, name::Any = nothing)
    filter_by = Dict{Symbol,Vector{String}}()

    if !isnothing(population)
        filter_by[:population] = typeof(population) <: AbstractArray ? string.(population) : [string(population)]
        err = filter_by[:population][filter_by[:population] .∉ Ref(unique(data.meta.population))]
        if length(err) > 0
            printstyled("Populations not found: ", bold = true)
            print("\"" * err[1] * "\"")
            if length(err) > 1
                [print(", \"$i\"") for i in err[2:end]] ; print("\n")
            end
        end
    end
    if !isnothing(locus)
        filter_by[:locus] = typeof(locus) <: AbstractArray ? string.(locus) : [string(locus)]
        err = filter_by[:locus][filter_by[:locus] .∉ Ref(loci(data))]
        if length(err) > 0
            printstyled("Loci not found: ", bold = true)
            print("\"" * err[1] * "\"")
            if length(err) > 1
                [print(", \"$i\"") for i in err[2:end]] ; print("\n")
            end
        end
    end
    if !isnothing(name)
        filter_by[:name] = typeof(name) <: AbstractArray ? string.(name) : [string(name)]
        err = filter_by[:name][filter_by[:name] .∉ Ref(data.meta.name)]
        if length(err) > 0
            printstyled("Samples not found: ", bold = true)
            print("\"" * err[1] * "\"")
            if length(err) > 1
                [print(", \"$i\"") for i in err[2:end]] ; print("\n")
            end
        end
    end

    filter_keys = Symbol.(keys(filter_by))
    meta_keys = filter_keys[filter_keys .!= :locus]

    if length(filter_keys) == 1
        filter!(filter_keys[1] => x -> x ∉ filter_by[filter_keys[1]] , data.loci)
        if !isempty(meta_keys)
            filter!(meta_keys[1] => x -> x ∉ filter_by[meta_keys[1]] , data.meta)
        end
    elseif length(filter_keys) == 2
        filter!([filter_keys[1], filter_keys[2]] => (x,y) -> x ∉ filter_by[filter_keys[1]] && y ∉ filter_by[filter_keys[2]] , data.loci)
        if !isempty(meta_keys)
            [filter!(i => x -> x ∉ filter_by[i] , data.meta) for i in meta_keys]
        end
    elseif length(filter_keys) == 3
        filter!([filter_keys[1], filter_keys[2], filter_keys[3]] => (x,y,z) -> x ∉ filter_by[filter_keys[1]] && y ∉ filter_by[filter_keys[2]] && z ∉ filter_by[filter_keys[3]] , data.loci)
        if !isempty(meta_keys)
            [filter!(i => x -> x ∉ filter_by[i] , data.meta) for i in meta_keys]
        end
    else
        throw(ArgumentError("Please specify at least one filter parameter of population, locus, or name"))
    end
    # make sure to update all the PooledArray pools
    [data.loci[!, i] = PooledArray(Array(data.loci[!, i])) for i in [:population, :locus, :name]]
    return
end

const omit! = exclude!
const remove! = exclude!

"""
    exclude(data::PopData, kwargs...)
Returns a new `PopData` object excluding all occurrences of the specified keywords.
The keywords can be used in any combination. Synonymous with `omit` and `remove`. All
values are converted to `String` for filtering, so `Symbol` and numbers will also work.

### Keyword Arguments
#### `locus`
A `String` or `Vector{String}` of loci you want to remove from the `PopData`.

#### `population`
A `String` or `Vector{String}` of populations you want to remove from the `PopData`.

#### `name`
A `String` or `Vector{String}` of samples you want to remove from the `PopData`.

**Examples**
```
cats = @nancycats;
exclude(cats, name = "N100", population = 1:5)
exclude(cats, name = ["N100", "N102", "N211"], locus = ["fca8", "fca23"])
exclude(cats, name = "N102", locus = :fca8, population = "3")
```
"""
function exclude(data::PopData; population::Any = nothing, locus::Any = nothing, name::Any = nothing)
    tmp = copy(data)
    exclude!(tmp; population = population, locus = locus, name = name)
    return tmp
end

const omit = exclude
const remove = exclude

"""
    keep(data::PopData, kwargs...)
Edit a `PopData` object in-place by keeping only the occurrences of the specified keyword.
Unlike `exclude!()`. only one keyword can be used at a time. All values are
converted to `String` for filtering, so `Symbol` and numbers will also work.

### Keyword Arguments
#### `locus`
A `String` or `Vector{String}` of loci you want to keep in the `PopData`.

#### `population`
A `String` or `Vector{String}` of populations you want to keep in the `PopData`.

#### `name`
A `String` or `Vector{String}` of samples you want to keep in the `PopData`.

**Examples**
```
cats = @nancycats;
keep(cats, population = 1:5)
keep(cats, name = ["N100", "N102", "N211"])
keep(cats, locus = [:fca8, "fca37"])
"""
function keep!(data::PopData; population::Any = nothing, locus::Any = nothing, name::Any = nothing)
    count(!isnothing, [population, locus, name]) > 1 && throw(ArgumentError("Please specify only one of \`population\`, \`locus\`, or \`name\` keyword arguments"))

    filter_by = Dict{Symbol,Vector{String}}()

    if !isnothing(population)
        filter_by[:population] = typeof(population) <: AbstractArray ? string.(population) : [string(population)]
        err = filter_by[:population][filter_by[:population] .∉ Ref(unique(data.meta.population))]
        if length(err) > 0
            notice = "Criteria not found in PopData\nPopulations: "
            notice *= "\"" * err[1] * "\""
            if length(err) > 1
                [notice *= ", \"$i\"" for i in err[2:end]]
            end
            throw(ArgumentError(notice))
        end
    end
    if !isnothing(locus)
        filter_by[:locus] = typeof(locus) <: AbstractArray ? string.(locus) : [string(locus)]
        err = filter_by[:locus][filter_by[:locus] .∉ Ref(loci(data))]
        if length(err) > 0
            notice = "Criteria not found in PopData\nLoci: "
            notice *= "\"" * err[1] * "\""
            if length(err) > 1
                [notice *= ", \"$i\"" for i in err[2:end]]
            end
            throw(ArgumentError(notice))
        end
    end
    if !isnothing(name)
        filter_by[:name] = typeof(name) <: AbstractArray ? string.(name) : [string(name)]
        err = filter_by[:name][filter_by[:name] .∉ Ref(data.meta.name)]
        if length(err) > 0
            notice = "Criteria not found in PopData\nSamples: "
            notice *= "\"" * err[1] * "\""
            if length(err) > 1
                [notice *= ", \"$i\"" for i in err[2:end]]
            end
            throw(ArgumentError(notice))
        end
    end

    filter_keys = Symbol.(keys(filter_by))
    meta_keys = filter_keys[filter_keys .!= :locus]

    filter!(filter_keys[1] => x -> x ∈ filter_by[filter_keys[1]] , data.loci)
    if !isempty(meta_keys)
        filter!(meta_keys[1] => x -> x ∈ filter_by[meta_keys[1]] , data.meta)
    end
    [data.loci[!, i] = PooledArray(Array(data.loci[!, i])) for i in [:population, :locus, :name]]
    return
end

"""
    keep(data::PopData, kwargs...)
Returns a new `PopData` object keeping only the occurrences of the specified keyword.
Unlike `exclude()`. only one keyword can be used at a time. All values are
converted to `String` for filtering, so `Symbol` and numbers will also work.

### Keyword Arguments
#### `locus`
A `String` or `Vector{String}` of loci you want to keep in the `PopData`.

#### `population`
A `String` or `Vector{String}` of populations you want to keep in the `PopData`.

#### `name`
A `String` or `Vector{String}` of samples you want to keep in the `PopData`.

**Examples**
```
cats = @nancycats;
keep(cats, population = 1:5)
keep(cats, name = ["N100", "N102", "N211"])
keep(cats, locus = [:fca8, "fca37"])
```
"""
function keep(data::PopData; population::Any = nothing, locus::Any = nothing, name::Any = nothing)
    tmp = copy(data)
    keep!(tmp; population = population, locus = locus, name = name)
    return tmp
end

"""
    samples(data::PopData)
View individual/sample names in a `PopData`
"""
function samples(data::PopData)
    @view data.meta[!, :name]
end


text:
using Wflow
using Documenter

pages = [
    "Home" => "index.md",
    "Quick start" => "quick-start.md",
    "Model structure" => "structure.md",
    "Models" => [
        "wflow_sbm" => "model/sbm.md",
        "wflow_hbv" => "model/hbv.md",
    ],
    "Vertical components" => [
        "SBM" => "vertical/sbm.md",
        "HBV" => "vertical/hbv.md",
        "Vertical processes" => "vertical/process.md",
    ],
    "Lateral components" => [
        "Kinematic wave" => "lateral/kinwave.md",
        "Groundwater flow" => "lateral/gwf.md",
    ],
    "Developer documentation" => "devdocs.md",
]

makedocs(;
    modules = [Wflow],
    authors = "Deltares and contributors",
    repo = "https://github.com/Deltares/Wflow.jl/blob/{commit}{path}#L{line}",
    sitename = "Wflow.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://deltares.github.io/Wflow.jl",
        assets = String[],
    ),
    pages = pages
)

deploydocs(; repo = "github.com/Deltares/Wflow.jl")


text:
__precompile__()

"""
Finder is a module that find freelancers near
"""
module Finder

using MyGraph, Person, DataStructures, Geodesy

export  FinderGraph, Person,
        availableworkers, availableworkersbyprof, findworkers!

include("findergraph.jl")

end


text:
@testset "Run Function Tests" begin

  @test isdefined(Julz, :run) == true

  run_examples = Dict(
    "main" => "done.\n",
    "example_task" => "404\n"
  )

  output_file_name = "output.jl"

  for (cur_task, cur_output) in run_examples

    if isfile(output_file_name) ; rm(output_file_name) ; end

    originalSTDOUT = STDOUT

    (outRead, outWrite) = redirect_stdout()

    Julz.run(cur_task)

    close(outWrite)

    data = readavailable(outRead)

    close(outRead)

    redirect_stdout(originalSTDOUT)

    @test String(data) == cur_output

    if cur_task == "main"
      @test isfile(output_file_name)
    end

  end

  test_input = ["\"woof\""]
  test_output = "404 - woof\n"

  originalSTDOUT = STDOUT

  (outRead, outWrite) = redirect_stdout()

  Julz.run("example_task", test_input)

  close(outWrite)

  data = readavailable(outRead)

  close(outRead)

  redirect_stdout(originalSTDOUT)

  @test String(data) == test_output

end


text:
using DoseResponse
using Test

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


text:
macro c_cmd(exprs...) return clang_cmd(__module__, __source__, :c, exprs...) end
macro c_str(exprs...) return clang_str(__module__, __source__, :c, exprs...) end

const C_HEADER = joinpath(@__DIR__, "$(nameof(@__MODULE__)).h")
header(::Type{Context{:c}}) = C_HEADER
Base.String(::Type{Context{:c}}) = "C"

const C_PRIMS   = ("void", "bool", "char", "short", "int", "long", "float", "double", "signed", "unsigned", "_Bool", "_Complex", "const", "restrict", "volatile")
const C_PRIM    = Regex(raw"^\s*((?:"*join(C_PRIMS, '|')*raw")(?:\s+(?:"*join(C_PRIMS, '|')*raw"))*)\s*$")
const C_USER    = r"^\s*(?:(const|restrict|volatile)\s+)?(?:(struct|union|enum)\s+)?(\w+)(?:\s+(const|restrict|volatile))?\s*$"
const C_POINTER = r"^(.+)\s*\*(?:\s*(const|restrict|volatile))?\s*$"
const C_ARRAY   = r"^(.+)\s*\[(\d+)\]\s*$"
const C_TYPES   = raw"\w+(?:(?:[*\s]+)\w+)*(?:[*\s]+)?"
const C_MODULE  = r"\s*(\w+)\.(.*)$"
const C_TYPEVAR = Regex(raw"^\s*("*C_TYPES*raw")\s*$")
const C_FUNCPTR = Regex(raw"^\s*("*C_TYPES*raw")\s*\(\*\)\(([^)]*)\)(const)?(?:\s+__attribute__\(\(([^)]+)\)\))?\s*$")
function getref(ctx::Type{Context{:c}}, str::AbstractString, mod::Union{Expr, Nothing} = nothing)
	function decorate(t, set)
		t = "volatile" in set ? :(Cvolatile{$(t)}) : t
		t = "restrict" in set ? :(Crestrict{$(t)}) : t
		t = "const"    in set ?    :(Cconst{$(t)}) : t
		return t
	end

	function sign(t, set)
		t =   "signed" in set ?   :(signed($(t))) : t
		t = "unsigned" in set ? :(unsigned($(t))) : t
		return t
	end

	function array(t, szs = [])
		m = match(C_ARRAY, t)
		if !isnothing(m)
			push!(szs, parse(Int, m.captures[2]))
			return array(m.captures[1], szs)
		end

		isempty(szs) && return nothing

		inner = getref(ctx, t, mod)
		isnothing(inner) && return nothing

		for sz in szs
			inner = decorate(:(Carray{$(inner), $(sz)}), ())
		end
		return inner
	end

	arr = array(str)
	isnothing(arr) || return arr

	m = match(C_POINTER, str)
	if !isnothing(m)
		inner = getref(ctx, m.captures[1], mod)
		isnothing(inner) || return decorate(:(Cptr{$(inner)}), (m.captures[2],))
	end

	m = match(C_FUNCPTR, str)
	if !isnothing(m)
		rettype  = getref(ctx, m.captures[1], mod)
		argtypes = isempty(m.captures[2]) ? () : map(arg -> getref(ctx, arg, mod), split(m.captures[2], ','))
		decor    = (m.captures[3],)
		conv     = isnothing(m.captures[4]) ? () : (QuoteNode(Symbol(m.captures[4])),)

		isnothing(rettype) || any(isnothing, argtypes) || return decorate(:(Cptr{Cfunction{$(rettype), Tuple{$(argtypes...)}, $(conv...)}}), decor)
	end

	m = match(C_PRIM, str)
	if !isnothing(m)
		prims = split(m.captures[1])

		if "_Complex" in prims
			if "double" in prims
				if "long" in prims
					return decorate(:(Ccomplex{Cdouble}), prims)
				else
					return decorate(:(Ccomplex{Clongdouble}), prims)
				end
			elseif "float" in prims
				return decorate(:(Ccomplex{Cfloat}), prims)
			end
		elseif "double" in prims
			if "long" in prims
				return decorate(:(Cdouble), prims)
			else
				return decorate(:(Clongdouble), prims)
			end
		elseif "float" in prims
			return decorate(:(Cfloat), prims)
		elseif "void" in prims
			return decorate(:(Cvoid), prims)
		elseif "bool" in prims || "_Bool" in prims
			return decorate(:(Cbool), prims)
		elseif "char" in prims
			return decorate(sign(:(Cchar), prims), prims)
		elseif "short" in prims
			return decorate(sign(:(Cshort), prims), prims)
		elseif "int" in prims
			return decorate(sign(:(Cint), prims), prims)
		elseif "long" in prims
			if count(isequal("long"), prims) > 1
				return decorate(sign(:(Clonglong), prims), prims)
			else
				return decorate(sign(:(Clong), prims), prims)
			end
		elseif "signed" in prims
			return decorate(:(signed(Cint)), prims)
		elseif "unsigned" in prims
			return decorate(:(unsigned(Cint)), prims)
		end
	end

	m = match(C_USER, str)
	if !isnothing(m)
		decor = (m.captures[1], m.captures[4])
		inner = getname(ctx, isnothing(m.captures[2]) ? "$(m.captures[3])" : "$(m.captures[2]) $(m.captures[3])")
		isnothing(inner) || return decorate(isnothing(mod) ? esc(inner) : :($(mod).$(inner)), decor)
	end

	strip(str) == "..." && return :(Cvariadic)

	m = match(C_MODULE, str)
	if !isnothing(m)
		str = m.captures[2]
		m   = Symbol(m.captures[1])
		return getref(ctx, str, isnothing(mod) ? esc(m) : :($(mod).$(m)))
	end

	return nothing
end

function gettype(ctx::Type{Context{:c}}, type::CXType; kwargs...)
	if type.kind == CXType_Void
		result = :(Cvoid)
	elseif type.kind == CXType_Bool
		result = :(Cbool)
	elseif type.kind == CXType_Char_S
		result = :(Cchar)
	elseif type.kind == CXType_Char_U
		result = :(Cuchar)
	elseif type.kind == CXType_SChar
		result = :(Cchar)
	elseif type.kind == CXType_UChar
		result = :(Cuchar)
	elseif type.kind == CXType_Short
		result = :(Cshort)
	elseif type.kind == CXType_Int
		result = :(Cint)
	elseif type.kind == CXType_Long
		result = :(Clong)
	elseif type.kind == CXType_LongLong
		result = :(Clonglong)
	elseif type.kind == CXType_UShort
		result = :(Cushort)
	elseif type.kind == CXType_UInt
		result = :(Cuint)
	elseif type.kind == CXType_ULong
		result = :(Culong)
	elseif type.kind == CXType_ULongLong
		result = :(Culonglong)
	elseif type.kind == CXType_Float
		result = :(Cfloat)
	elseif type.kind == CXType_Double
		result = :(Cdouble)
	elseif type.kind == CXType_LongDouble
		result = :(Clongdouble)
	elseif type.kind in (
		CXType_Typedef,
		CXType_Record,
		CXType_Elaborated,
		CXType_Enum,
	)
		decl = clang_getTypeDeclaration(type)
		t = clang_getCursorType(decl)
		result = gettype(ctx, string(t); kwargs...)
	elseif type.kind == CXType_Complex
		cplxtype = clang_getElementType(type)
		cplxtype = gettype(ctx, cplxtype; kwargs...)
		result = :(Ccomplex{$(cplxtype)})
	elseif type.kind == CXType_Pointer
		ptrtype = clang_getPointeeType(type)
		ptrtype = gettype(ctx, ptrtype; kwargs...)
		result = :(Cptr{$(ptrtype)})
	elseif type.kind == CXType_ConstantArray
		num = clang_getNumElements(type)
		arrtype = clang_getElementType(type)
		arrtype = gettype(ctx, arrtype; kwargs...)
		result = :(Carray{$(arrtype), $(num)})
	elseif type.kind == CXType_IncompleteArray
		arrtype = clang_getElementType(type)
		arrtype = gettype(ctx, arrtype; kwargs...)
		result = :(Carray{$(arrtype), 0})
	elseif type.kind in (
		CXType_FunctionProto,
		CXType_FunctionNoProto,
	)
		rettype = clang_getResultType(type)
		rettype = gettype(ctx, rettype; kwargs...)

		num = clang_getNumArgTypes(type)
		argtypes = map(1:num) do ind
			argtype = clang_getArgType(type, ind-1)
			return gettype(ctx, argtype; kwargs...)
		end
		argtypes = Bool(clang_isFunctionTypeVariadic(type)) ? [argtypes..., Cvariadic] : argtypes

		conv = clang_getFunctionTypeCallingConv(type)
		if conv == CXCallingConv_C
			conv = :(cdecl)
		else
			error("Unhandled calling convention $(conv)")
		end

		result = :(Cfunction{$(rettype), Tuple{$(argtypes...)}, $(QuoteNode(conv))})
	else
		error("Unhandle type $(type.kind)")
	end

	if Bool(clang_isVolatileQualifiedType(type))
		# TODO: not handled in get/set property etc. yet
		# result = :(Cvolatile{$(result)})
	end

	if Bool(clang_isRestrictQualifiedType(type))
		# TODO: not handled in get/set property etc. yet
		# result = :(Crestrict{$(result)})
	end

	# NOTE: always want Cconst around everything else to make life easier in get/set
	if Bool(clang_isConstQualifiedType(type))
		result = :(Cconst{$(result)})
	end

	return result
end

function getexprs(ctx::Context{:c}, cursor::CXCursor)
	exprs = []

	getblock(ctx).flags.skip && return exprs

	if cursor.kind == CXCursor_TranslationUnit
		append!(exprs, getexprs_tu(ctx, cursor))
	elseif cursor.kind == CXCursor_TypedefDecl
		push!(exprs, getexprs_typedef(ctx, cursor))
	elseif cursor.kind in (
		CXCursor_StructDecl,
		CXCursor_UnionDecl,
		CXCursor_EnumDecl,
	)
		append!(exprs, getexprs_opaque(ctx, cursor))
	elseif cursor.kind in (
		CXCursor_VarDecl,
		CXCursor_FunctionDecl,
	)
		append!(exprs, getexprs_binding(ctx, cursor))
	elseif cursor.kind == CXCursor_ParmDecl
		for child in children(cursor)
			if child.kind in (
				CXCursor_StructDecl,
				CXCursor_UnionDecl,
				CXCursor_EnumDecl,
				CXCursor_TypedefDecl,
			)
				append!(exprs, getexprs(ctx, child))
			end
		end
	elseif cursor.kind == CXCursor_InclusionDirective
		append!(exprs, getexprs_include(ctx, cursor))
	elseif !(cursor.kind in (
		CXCursor_FieldDecl,
		CXCursor_EnumConstantDecl,
		CXCursor_ParmDecl,
		CXCursor_MacroDefinition,  # TODO: support this
		CXCursor_MacroInstantiation,  # TODO: support this
	))
		@warn "Not implemented $(cursor.kind):  $(string(cursor))"
	end

	return exprs
end

function getexprs_typedef(ctx::Context{:c}, cursor::CXCursor)
	type  = clang_getTypedefDeclUnderlyingType(cursor)
	if type.kind == CXType_Elaborated
		# get from elaborated to actual record definition
		decl  = clang_getTypeDeclaration(type)
		type  = clang_getCursorType(decl)
	end

	tname = gettype(ctx, type)
	name  = gettype(ctx, string(cursor))
	jlsym = getjl(ctx, string(cursor))
	docs  = getdocs(ctx, cursor)

	# nested anonymous opaque type had this typedef's name and provided these exprs
	return name == tname ? quote end : getexprs(ctx, ((name, jlsym, docs),), quote
		const $(name) = $(tname)
	end)
end

function getexprs_opaque(ctx::Context{:c}, cursor::CXCursor)
	exprs = []

	type = clang_getCursorType(cursor)
	decl = clang_getTypeDeclaration(type)

	size = clang_Type_getSizeOf(type)
	isdecl = Bool(clang_isCursorDefinition(cursor))
	isanon = Bool(clang_Cursor_isAnonymous(cursor))

	kind =
		cursor.kind == CXCursor_EnumDecl ? Cenum :
		cursor.kind == CXCursor_StructDecl ? Cstruct : Cunion
	(name, absname) = isanon ?
		(gettype(ctx, type), kind) :
		(gettype(ctx, type, isbitstype = true), gettype(ctx, type))
	(jlsym, jlabssym, docs) = isanon ?
		(getjl(ctx, string(type)), nothing, nothing) :
		(nothing, getjl(ctx, string(type)), getdocs(ctx, cursor))

	cursorLoc = getlocation(cursor)
	declLoc   = getlocation(decl)
	islater = !isdecl && first(cursorLoc).file == first(declLoc).file
	if !isanon && (isdecl || islater)
		push!(exprs, getexprs(ctx, ((absname, jlabssym, islater ? nothing : docs),), quote
			abstract type $(absname) <: $(kind) end
		end))
	end

	if isdecl
		common = quote
			C.bitstype(::Type{$(name)}) = $(name)
			$(!isanon ? :(C.bitstype(::Type{$(absname)}) = $(name)) : :())
		end

		if cursor.kind == CXCursor_EnumDecl
			push!(exprs, getexprs_enum(ctx, cursor, name, absname, jlsym, jlabssym, docs, common))
		else
			for child in children(cursor)
				if child.kind in (
					CXCursor_StructDecl,
					CXCursor_UnionDecl,
					CXCursor_EnumDecl,
					CXCursor_TypedefDecl,
				)
					append!(exprs, getexprs(ctx, child))
				end
			end
			push!(exprs, getexprs_aggregate(ctx, cursor, name, absname, jlsym, jlabssym, docs, common))
		end
	end

	return exprs
end

function getexprs_enum(ctx::Context{:c}, cursor::CXCursor, name, absname, jlsym, jlabssym, docs, common)
	type = clang_getCursorType(cursor)
	size = clang_Type_getSizeOf(type)
	inttype = clang_getEnumDeclIntegerType(cursor)
	itype = gettype(ctx, inttype)

	values = map(filter(child -> child.kind == CXCursor_EnumConstantDecl, children(cursor))) do val
		j = Symbol(string(val))
		v = inttype.kind in (
			CXType_UChar,
			CXType_Char_U,
			CXType_UShort,
			CXType_UInt,
			CXType_ULong,
			CXType_ULongLong,
		) ? clang_getEnumConstantDeclUnsignedValue(val) : clang_getEnumConstantDeclValue(val)
		d = getdocs(ctx, val)

		return (j, v, d)
	end

	return getexprs(ctx, (
		(!isnothing(jlabssym) ? ((absname, jlabssym, nothing),) : ())...,
		(name, jlsym, docs),
		map(((j, v, d),) -> (esc(getname(ctx, j)), esc(j), d), values)...,
	), quote
		primitive type $(name) <: $(absname) $(size*8) end
		(::$(Type){$(name)})($(esc(:val))::Integer = zero($(itype))) = Core.Intrinsics.bitcast($(name), convert($(itype), $(esc(:val))))

		# values = Tuple{
		# 	Tuple{:name, value},  # name and literal value
		# }
		C.values(::Type{$(name)}) = Tuple{$(map(((j, v, d),) -> :(Tuple{$(QuoteNode(getname(ctx, j))), $(v)}), values)...)}
		Base.Enums.basetype(::Type{$(name)}) = $(itype)
	end, common, map(((j, v, d),) -> quote
		const $(esc(getname(ctx, j))) = $(name)($(v))
	end, values)...)
end

function getexprs_aggregate(ctx::Context{:c}, cursor::CXCursor, name, absname, jlsym, jlabssym, docs, common)
	type = clang_getCursorType(cursor)
	size = clang_Type_getSizeOf(type)

	# children are ordered with type defined before the field using it, so we iterate in reverse to add the field first and then not add the type's fields in that case
	fields = []
	lastField = nothing
	function getfields(kids)
		for c in reverse(eachindex(kids))
			child = kids[c]
			if child.kind == CXCursor_FieldDecl
				t = clang_getCursorType(child)
				d = clang_getTypeDeclaration(t)

				if t.kind in (
					CXType_Void, CXType_Bool,
					CXType_Char_S, CXType_Char_U,
					CXType_SChar, CXType_UChar,
					CXType_Short, CXType_UShort,
					CXType_Int, CXType_UInt,
					CXType_Long, CXType_ULong,
					CXType_LongLong, CXType_ULongLong,
					CXType_Float, CXType_Double, CXType_LongDouble,
					CXType_Complex,
				)
					lastField = nothing
				elseif t.kind in (
					CXType_Typedef,
					CXType_Elaborated,
				)
					lastField = clang_getTypeDeclaration(clang_getCursorType(clang_getTypeDeclaration(t)))
				elseif t.kind in (
					CXType_ConstantArray,
					CXType_IncompleteArray,
				)
					lastField = clang_getTypeDeclaration(clang_getElementType(t))
				elseif t.kind == CXType_Pointer
					lastField = clang_getTypeDeclaration(clang_getPointeeType(t))
				else
					error("Unhandled field nested type $(t.kind)")
				end

				n = string(child)
				isempty(n) && continue

				b = -1
				o = clang_Type_getOffsetOf(type, n)
				if Bool(clang_Cursor_isBitField(child))
					b = clang_getFieldDeclBitWidth(child)
					b > 0 || continue
				end

				t = gettype(ctx, t)
				push!(fields, :(Tuple{$(QuoteNode(Symbol(n))), Tuple{$(t), bitstype($(t)), integraltype($(t)), $(o÷8), $(o & 0b111), $(b)}}))
			elseif child.kind in (
				CXCursor_StructDecl,
				CXCursor_UnionDecl,
			)
				if Bool(clang_Cursor_isAnonymous(child))
					if isnothing(lastField) || !Bool(clang_equalCursors(lastField, child))
						getfields(children(child))
					end
				end
			end
		end
	end
	getfields(children(cursor))
	reverse!(fields)

	return getexprs(ctx, (
		(!isnothing(jlabssym) ? ((absname, jlabssym, nothing),) : ())...,
		(name, jlsym, docs),
	), quote
		struct $(name) <: $(absname)
			mem::NTuple{$(size), UInt8}

			$(name)(::UndefInitializer, $(esc(:mem))::NTuple{$(size), UInt8}) = new($(esc(:mem)))
		end

		# fields = Tuple{
		# 	# a field's name, type, type with const/volatile/restrict removed, integral equivalent, bytes+bits offset, and bitfield size (-1 for non-bitfield)
		# 	Tuple{:name, Tuple{FieldType, UnqualifiedFieldType, IntegerType, bytes, bits, size}},
		# }
		C.fields(::Type{$(name)}) = Tuple{$(fields...)}
	end, common)
end

function getexprs_binding(ctx::Context{:c}, cursor::CXCursor)
	exprs = []

	for child in children(cursor)
		if child.kind in (
			CXCursor_StructDecl,
			CXCursor_UnionDecl,
			CXCursor_EnumDecl,
			CXCursor_TypedefDecl,
			CXCursor_ParmDecl,
		)
			append!(exprs, getexprs(ctx, child))
		end
	end

	exposed = clang_Cursor_getStorageClass(cursor)
	visible = clang_getCursorVisibility(cursor)
	linkage = clang_getCursorLinkage(cursor)
	if exposed in (
		CX_SC_None,
		CX_SC_Extern,
	) && visible in (
		CXVisibility_Default,
		CXVisibility_Protected,
	) && linkage in (
		CXLinkage_External,
	)
		name = string(cursor)
		lib = getlib(ctx, name)
		jlsym = getjl(ctx, name)
		sym = gettype(ctx, name)
		docs = getdocs(ctx, cursor)

		type = clang_getCursorType(cursor)
		if cursor.kind == CXCursor_VarDecl
			cb = :(Cbinding{$(gettype(ctx, type)), $(QuoteNode(Symbol(lib))), $(QuoteNode(Symbol(name)))})
			push!(exprs, getexprs(ctx, ((sym, jlsym, docs),), quote
				const $(sym) = $(cb)()
			end))
		elseif Bool(clang_Cursor_isFunctionInlined(cursor))
			getblock(ctx).flasg.quiet || @warn "Skipping inline function `$(name)`"
		else
			rettype = clang_getResultType(type)

			num = clang_Cursor_getNumArguments(cursor)
			args = map(1:num) do ind
				arg = clang_Cursor_getArgument(cursor, ind-1)

				argname = string(arg)
				argname = isempty(argname) ? getjl(ctx, "arg$(ind)") : gettype(ctx, argname)
				argtype = clang_getCursorType(arg)

				return argname => gettype(ctx, argtype)
			end
			args = Bool(clang_Cursor_isVariadic(cursor)) ? [args..., :($(getjl(ctx, :vararg))...) => :(Cvariadic)] : args

			conv = clang_getFunctionTypeCallingConv(type)
			if conv == CXCallingConv_C
				conv = :(cdecl)
			else
				error("Unhandled calling convention $(conv)")
			end

			func = getjl(ctx, :func)
			cb = :(Cbinding{Cfunction{$(gettype(ctx, rettype)), Tuple{$(map(last, args)...)}, $(QuoteNode(conv))}, $(lib), $(QuoteNode(Symbol(name)))})
			push!(exprs, getexprs(ctx, ((sym, jlsym, docs),), quote
				const $(sym) = $(cb)()
				($(func)::typeof($(sym)))($(map(first, args)...),) = funccall($(func), $(map(first, args)...),)
			end))
		end
	end

	return exprs
end


text:
# This file contains functions to get the list of all available games in a season

"""
    PublicWyscoutGames
"""
Base.@kwdef struct PublicWyscoutGames
    data::DataFrame
end

"""
function games()
Return a dataframe with all available games in a season.
"""
function games(events_data::PublicWyscoutLoader, competition_id::Int, season_id::Int)::PublicWyscoutGames

    data = filter("competition_id" => ==(competition_id),events_data.match_index)
    data = filter("season_id"=> ==(season_id), data)
    data = unique(data[:,[:competition_id, :season_id, :db_matches]])

    match_data = get_matchs(events_data, data[1,:db_matches])
    rename!(match_data, :seasonId => :season_id,
                        :gameweek => :game_day,
                        :competitionId => :competition_id,
                        :wyId => :game_id,
                        :dateutc => :game_date)

    match_data = leftjoin(match_data,data, on = [:competition_id, :season_id])
    insertcols!(match_data, :home_team_id => 0, :away_team_id => 0)
    for i in eachindex(match_data.season_id)
        teams_id = names(DataFrame(match_data[i, :teamsData]))
        match_data[i, :home_team_id] = parse(Int,teams_id[1])
        match_data[i, :away_team_id] = parse(Int,teams_id[2])
    end

    match_data = match_data[:,[:season_id, :competition_id, :game_day, :game_id, :game_date,:home_team_id, :away_team_id]]
    match_data = PublicWyscoutGames(data  = match_data)

    return match_data
end

"""
    function get_matchs()
"""
function get_matchs(loader::PublicWyscoutLoader,matches::String)

    zarchive = ZipFile.Reader(joinpath(loader.path,"matches.zip"))
    dictio = Dict(zarchive.files[i].name => i for i in eachindex(zarchive.files))
    file_num = dictio[matches]
    data = DataFrame(jsontable(JSON3.read(read(zarchive.files[file_num]))))

    return data
end


text:
int main()
{
  printInt(+32 % -10);
  printInt(-32 % -10);
  printInt(+32 % +10);
  printInt(-32 % +10);
  printString("");
  printInt(+38 % -10);
  printInt(-38 % -10);
  printInt(+38 % +10);
  printInt(-38 % +10);
  printString("");
  printDouble(+32.0 % -10.0);
  printDouble(-32.0 % -10.0);
  printDouble(+32.0 % +10.0);
  printDouble(-32.0 % +10.0);
  printString("");
  printDouble(+38.0 % -10.0);
  printDouble(-38.0 % -10.0);
  printDouble(+38.0 % +10.0);
  printDouble(-38.0 % +10.0);
  return 0;
}

/* vim:set ts=2 sts=2 sw=2 et ft=c: */

max_stars_repo_name: jwilk/jtc

max_issues_repo_name: jwilk/jtc

max_forks_repo_name: jwilk/jtc


text:
include(joinpath("..", "src", "axelrodnet.jl"))

Axelrodnet.run_experiment(
    experiment_name="test_run",
    agentcount=100,
    n_iter=1_000_000,
    nettopology="grid",
    networkprops=Dict("grid_height" => 10),
    stubbornfrac=0.0,
    rndseed=0,
    repcount=10,
    export_every_n=100_000,
    export_experiment=false,
    exportmode="aggregated",
    keep_rawnet=true
)

max_stars_repo_name: JohannesNakayama/axelrodnet

max_issues_repo_name: JohannesNakayama/axelrodnet

max_forks_repo_name: JohannesNakayama/axelrodnet


text:
using PH2514
using Test

@testset "PH2514.jl" begin
end


text:
function SmallManeuverAnalyticalModel(arg0::SpacecraftState, arg1::Frame, arg2::Vector3D, arg3::jdouble)
    return SmallManeuverAnalyticalModel((SpacecraftState, Frame, Vector3D, jdouble), arg0, arg1, arg2, arg3)
end

function SmallManeuverAnalyticalModel(arg0::SpacecraftState, arg1::Vector3D, arg2::jdouble)
    return SmallManeuverAnalyticalModel((SpacecraftState, Vector3D, jdouble), arg0, arg1, arg2)
end

function apply(obj::SmallManeuverAnalyticalModel, arg0::Orbit)
    return jcall(obj, "apply", Orbit, (Orbit,), arg0)
end

function apply(obj::SmallManeuverAnalyticalModel, arg0::SpacecraftState)
    return jcall(obj, "apply", SpacecraftState, (SpacecraftState,), arg0)
end

function get_date(obj::SmallManeuverAnalyticalModel)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_inertial_dv(obj::SmallManeuverAnalyticalModel)
    return jcall(obj, "getInertialDV", Vector3D, ())
end

function get_inertial_frame(obj::SmallManeuverAnalyticalModel)
    return jcall(obj, "getInertialFrame", Frame, ())
end

function get_jacobian(obj::SmallManeuverAnalyticalModel, arg0::Orbit, arg1::PositionAngle, arg2::Vector{Vector{jdouble}})
    return jcall(obj, "getJacobian", void, (Orbit, PositionAngle, Vector{Vector{jdouble}}), arg0, arg1, arg2)
end

function update_mass(obj::SmallManeuverAnalyticalModel, arg0::jdouble)
    return jcall(obj, "updateMass", jdouble, (jdouble,), arg0)
end