text:
struct KernelFunctionOperation{LX, LY, LZ, P, G, T, K, D} <: AbstractOperation{LX, LY, LZ, G, T}
    kernel_function :: K
    computed_dependencies :: D
    parameters :: P
    grid :: G

    function KernelFunctionOperation{LX, LY, LZ}(kernel_function::K, computed_dependencies::D,
                                                 parameters::P, grid::G) where {LX, LY, LZ, K, G, D, P}
        T = eltype(grid)
        return new{LX, LY, LZ, P, G, T, K, D}(kernel_function, computed_dependencies, parameters, grid)
    end

end

"""
    KernelFunctionOperation{LX, LY, LZ}(kernel_function, grid;
                                        computed_dependencies=(), parameters=nothing)

Constructs a `KernelFunctionOperation` at location `(LX, LY, LZ)` on `grid` an with
an optional iterable of `computed_dependencies` and arbitrary `parameters`.

With `isnothing(parameters)` (the default), `kernel_function` is called with

```julia
kernel_function(i, j, k, grid, computed_dependencies...)
```

Otherwise `kernel_function` is called with

```julia
kernel_function(i, j, k, grid, computed_dependencies..., parameters)
```

Examples
========

Construct a kernel function operation that returns random numbers:

```julia
random_kernel_function(i, j, k, grid) = rand() # use CUDA.rand on the GPU

kernel_op = KernelFunctionOperation{Center, Center, Center}(random_kernel_function, grid)
```

Construct a kernel function operation using the vertical vorticity operator
valid on curvilinear and cubed sphere grids:

```julia
using Oceananigans.Operators: ζ₃ᶠᶠᵃ # called with signature ζ₃ᶠᶠᵃ(i, j, k, grid, u, v)

grid = model.grid
u, v, w = model.velocities

ζ_op = KernelFunctionOperation{Face, Face, Center}(ζ₃ᶠᶠᵃ, grid, computed_dependencies=(u, v))
```
"""
function KernelFunctionOperation{LX, LY, LZ}(kernel_function, grid;
                                            computed_dependencies = (),
                                            parameters = nothing) where {LX, LY, LZ}

    return KernelFunctionOperation{LX, LY, LZ}(kernel_function, computed_dependencies, parameters, grid)
end

@inline Base.getindex(κ::KernelFunctionOperation, i, j, k) = κ.kernel_function(i, j, k, κ.grid, κ.computed_dependencies..., κ.parameters)
@inline Base.getindex(κ::KernelFunctionOperation{LX, LY, LZ, <:Nothing}, i, j, k) where {LX, LY, LZ} = κ.kernel_function(i, j, k, κ.grid, κ.computed_dependencies...)

# Compute dependencies
compute_at!(κ::KernelFunctionOperation, time) = Tuple(compute_at!(d, time) for d in κ.computed_dependencies)

"Adapt `KernelFunctionOperation` to work on the GPU via CUDAnative and CUDAdrv."
Adapt.adapt_structure(to, κ::KernelFunctionOperation{LX, LY, LZ}) where {LX, LY, LZ} =
    KernelFunctionOperation{LX, LY, LZ}(Adapt.adapt(to, κ.kernel_function),
                                        Adapt.adapt(to, κ.computed_dependencies),
                                        Adapt.adapt(to, κ.parameters),
                                        Adapt.adapt(to, κ.grid))


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

name = "Libgpg_error"
version = v"1.36"

# Collection of sources required to build Libgpg-Error
sources = [
    ArchiveSource("https://gnupg.org/ftp/gcrypt/libgpg-error/libgpg-error-$(version.major).$(version.minor).tar.bz2",
                  "babd98437208c163175c29453f8681094bcaf92968a15cafb1a276076b33c97c"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libgpg-error-*/

for p in ${WORKSPACE}/srcdir/patches/*; do
    atomic_patch -p1 "${p}"
done

# "Generate" new arch platform definitions because they're actually really simple. :P
cp -iv src/syscfg/lock-obj-pub.i686-unknown-linux-gnu.h src/syscfg/lock-obj-pub.i686-unknown-linux-musl.h
cp -iv src/syscfg/lock-obj-pub.aarch64-unknown-linux-gnu.h src/syscfg/lock-obj-pub.aarch64-unknown-linux-musl.h
cp -iv src/syscfg/lock-obj-pub.arm-unknown-linux-gnueabi.h src/syscfg/lock-obj-pub.arm-unknown-linux-musleabihf.h
cp -iv src/syscfg/lock-obj-pub.x86_64-unknown-linux-gnu.h src/syscfg/lock-obj-pub.freebsd11.1.h

# Use libgpg-specific mapping for triplets
TARGET="${target}"
if [[ "${target}" == x86_64-apple-darwin* ]]; then
    TARGET=x86_64-apple-darwin
fi

# We need `msgformat`
apk add gettext

./configure --prefix=${prefix} --host=${TARGET} --build=${MACHTYPE}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.  We are manually disabling
# many platforms that do not seem to work.
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct(["libgpg-error", "libgpg-error6"], :libgpg_error),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

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

ext: jl


text:
using Oceananigans.TurbulenceClosures: ExplicitTimeDiscretization, VerticallyImplicitTimeDiscretization, z_viscosity
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBoundary

function relative_error(u_num, u, time)
    u_ans = Field(location(u_num), architecture(u_num), u_num.grid, nothing)
    set!(u_ans, (x, y, z) -> u(x, y, z, time))
    return mean((interior(u_num) .- interior(u_ans)).^2 ) / mean(interior(u_ans).^2)
end

function test_diffusion_simple(fieldname, timestepper, time_discretization)

    model = NonhydrostaticModel(timestepper = timestepper,
                                       grid = RegularRectilinearGrid(size=(1, 1, 16), extent=(1, 1, 1)),
                                    closure = IsotropicDiffusivity(ν=1, κ=1, time_discretization=time_discretization),
                                   coriolis = nothing,
                                    tracers = :c,
                                   buoyancy = nothing)

    field = get_model_field(fieldname, model)

    value = π
    interior(field) .= value

    for n in 1:10
        ab2_or_rk3_time_step!(model, 1, n)
    end

    field_data = interior(field)

    return !any(@. !isapprox(value, field_data))
end

function test_isotropic_diffusion_budget(fieldname, model)
    set!(model; u=0, v=0, w=0, c=0)
    set!(model; Dict(fieldname => (x, y, z) -> rand())...)

    field = get_model_field(fieldname, model)

    ν = z_viscosity(model.closure, nothing) # for generalizing to isotropic AnisotropicDiffusivity

    return test_diffusion_budget(fieldname, field, model, ν, model.grid.Δz)
end

function test_biharmonic_diffusion_budget(fieldname, model)
    set!(model; u=0, v=0, w=0, c=0)
    set!(model; Dict(fieldname => (x, y, z) -> rand())...)

    field = get_model_field(fieldname, model)

    return test_diffusion_budget(fieldname, field, model, model.closure.νz, model.grid.Δz, 4)
end

function test_diffusion_budget(fieldname, field, model, κ, Δ, order=2)
    init_mean = mean(interior(field))

    for n in 1:10
        # Very small time-steps required to bring error under machine precision
        ab2_or_rk3_time_step!(model, 1e-4 * Δ^order / κ, n)
    end

    final_mean = mean(interior(field))

    @info @sprintf("    Initial <%s>: %.16f, final <%s>: %.16f, final - initial: %.4e",
                   fieldname, init_mean, fieldname, final_mean, final_mean - init_mean)

    return isapprox(init_mean, final_mean)
end

function test_diffusion_cosine(fieldname, timestepper, grid, time_discretization)
    κ, m = 1, 2 # diffusivity and cosine wavenumber

    model = NonhydrostaticModel(timestepper = timestepper,
                                       grid = grid,
                                    closure = IsotropicDiffusivity(ν=κ, κ=κ, time_discretization=time_discretization),
                                   buoyancy = nothing)

    field = get_model_field(fieldname, model)

    zC = znodes(Center, grid, reshape=true)
    interior(field) .= cos.(m * zC)

    diffusing_cosine(κ, m, z, t) = exp(-κ * m^2 * t) * cos(m * z)

    # Step forward with small time-step relative to diff. time-scale
    Δt = 1e-6 * grid.Lz^2 / κ
    for n in 1:10
        ab2_or_rk3_time_step!(model, Δt, n)
    end

     numerical = interior(field)
    analytical = diffusing_cosine.(κ, m, zC, model.clock.time)

    return !any(@. !isapprox(numerical, analytical, atol=1e-6, rtol=1e-6))
end

function passive_tracer_advection_test(timestepper; N=128, κ=1e-12, Nt=100, background_velocity_field=false)
    L, U, V = 1.0, 0.5, 0.8
    δ, x₀, y₀ = L/15, L/2, L/2

    Δt = 0.05 * L/N / sqrt(U^2 + V^2)

    T(x, y, z, t) = exp( -((x - U*t - x₀)^2 + (y - V*t - y₀)^2) / (2*δ^2) )
    u₀(x, y, z) = U
    v₀(x, y, z) = V
    T₀(x, y, z) = T(x, y, z, 0)
    background_fields = Dict()

    if background_velocity_field
        background_fields[:u] = (x, y, z, t) -> U
        background_fields[:v] = (x, y, z, t) -> V
        u₀ = 0
        v₀ = 0
    end

    background_fields = NamedTuple{Tuple(keys(background_fields))}(values(background_fields))

    grid = RegularRectilinearGrid(size=(N, N, 2), extent=(L, L, L))
    closure = IsotropicDiffusivity(ν=κ, κ=κ)
    model = NonhydrostaticModel(timestepper=timestepper, grid=grid, closure=closure,
                                background_fields=background_fields)

    set!(model, u=u₀, v=v₀, T=T₀)

    for n in 1:Nt
        ab2_or_rk3_time_step!(model, Δt, n)
    end

    # Error tolerance is a bit arbitrary
    return relative_error(model.tracers.T, T, model.clock.time) < 1e-4
end

"""
Taylor-Green vortex test
See: https://en.wikipedia.org/wiki/Taylor%E2%80%93Green_vortex#Taylor%E2%80%93Green_vortex_solution
     and p. 310 of "Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Application"
     by Hesthaven & Warburton.
"""
function taylor_green_vortex_test(arch, timestepper, time_discretization; FT=Float64, N=64, Nt=10)
    Nx, Ny, Nz = N, N, 2
    Lx, Ly, Lz = 1, 1, 1
    ν = 1

    # Choose a very small time step as we are diffusion-limited in this test: Δt ≤ Δx² / 2ν
    Δx = Lx / Nx
    Δt = (1/10π) * Δx^2 / ν

    # Taylor-Green vortex analytic solution.
    @inline u(x, y, z, t) = -sin(2π*y) * exp(-4π^2 * ν * t)
    @inline v(x, y, z, t) =  sin(2π*x) * exp(-4π^2 * ν * t)

    model = NonhydrostaticModel(
        architecture = arch,
         timestepper = timestepper,
                grid = RegularRectilinearGrid(FT, size=(Nx, Ny, Nz), extent=(Lx, Ly, Lz)),
             closure = IsotropicDiffusivity(FT, ν=1, time_discretization=time_discretization),
             tracers = nothing,
            buoyancy = nothing)

    u₀(x, y, z) = u(x, y, z, 0)
    v₀(x, y, z) = v(x, y, z, 0)
    set!(model, u=u₀, v=v₀)

    for n in 1:Nt
        ab2_or_rk3_time_step!(model, Δt, n)
    end

    xF, yC, zC = nodes(model.velocities.u, reshape=true)
    xC, yF, zC = nodes(model.velocities.v, reshape=true)

    t = model.clock.time
    i = model.clock.iteration

    # Calculate relative error between model and analytic solutions for u and v.
    u_rel_err = abs.((interior(model.velocities.u) .- u.(xF, yC, zC, t)) ./ u.(xF, yC, zC, t))
    u_rel_err_avg = mean(u_rel_err)
    u_rel_err_max = maximum(u_rel_err)

    v_rel_err = abs.((interior(model.velocities.v) .- v.(xC, yF, zC, t)) ./ v.(xC, yF, zC, t))
    v_rel_err_avg = mean(v_rel_err)
    v_rel_err_max = maximum(v_rel_err)

    @info "Taylor-Green vortex test [$arch, $FT, Nx=Ny=$N, Nt=$Nt]: " *
          @sprintf("Δu: (avg=%6.3g, max=%6.3g), Δv: (avg=%6.3g, max=%6.3g)",
                   u_rel_err_avg, u_rel_err_max, v_rel_err_avg, v_rel_err_max)

    return u_rel_err_max < 5e-6 && v_rel_err_max < 5e-6
end

function stratified_fluid_remains_at_rest_with_tilted_gravity_buoyancy_tracer(arch, FT; N=32, L=2000, θ=60, N²=1e-5)
    topo = (Periodic, Bounded, Bounded)
    grid = RegularRectilinearGrid(FT, topology=topo, size=(1, N, N), extent=(L, L, L))

    g̃ = (0, sind(θ), cosd(θ))
    buoyancy = Buoyancy(model=BuoyancyTracer(), vertical_unit_vector=g̃)

    y_bc = GradientBoundaryCondition(N² * g̃[2])
    z_bc = GradientBoundaryCondition(N² * g̃[3])
    b_bcs = FieldBoundaryConditions(bottom=z_bc, top=z_bc, south=y_bc, north=y_bc)

    model = NonhydrostaticModel(
               architecture = arch,
                       grid = grid,
                   buoyancy = buoyancy,
                    tracers = :b,
                    closure = nothing,
        boundary_conditions = (b=b_bcs,)
    )

    b₀(x, y, z) = N² * (x*g̃[1] + y*g̃[2] + z*g̃[3])
    set!(model, b=b₀)

    simulation = Simulation(model, Δt=10minutes, stop_time=1hour)
    run!(simulation)

    ∂y_b = ComputedField(∂y(model.tracers.b))
    ∂z_b = ComputedField(∂z(model.tracers.b))

    compute!(∂y_b)
    compute!(∂z_b)

    mean_∂y_b = mean(∂y_b)
    mean_∂z_b = mean(∂z_b)

    Δ_y = N² * g̃[2] - mean_∂y_b
    Δ_z = N² * g̃[3] - mean_∂z_b

    @info "N² * g̃[2] = $(N² * g̃[2]), mean(∂y_b) = $(mean_∂y_b), Δ = $Δ_y at t = $(prettytime(model.clock.time)) with θ=$(θ)°"
    @info "N² * g̃[3] = $(N² * g̃[3]), mean(∂z_b) = $(mean_∂z_b), Δ = $Δ_z at t = $(prettytime(model.clock.time)) with θ=$(θ)°"

    @test N² * g̃[2] ≈ mean(∂y_b)
    @test N² * g̃[3] ≈ mean(∂z_b)

    CUDA.@allowscalar begin
        @test all(N² * g̃[2] .≈ interior(∂y_b))
        @test all(N² * g̃[3] .≈ interior(∂z_b))
    end

    return nothing
end

function stratified_fluid_remains_at_rest_with_tilted_gravity_temperature_tracer(arch, FT; N=32, L=2000, θ=60, N²=1e-5)
    topo = (Periodic, Bounded, Bounded)
    grid = RegularRectilinearGrid(FT, topology=topo, size=(1, N, N), extent=(L, L, L))

    g̃ = (0, sind(θ), cosd(θ))
    buoyancy = Buoyancy(model=SeawaterBuoyancy(), vertical_unit_vector=g̃)

    α  = buoyancy.model.equation_of_state.α
    g₀ = buoyancy.model.gravitational_acceleration
    ∂T∂z = N² / (g₀ * α)

    y_bc = GradientBoundaryCondition(∂T∂z * g̃[2])
    z_bc = GradientBoundaryCondition(∂T∂z * g̃[3])
    T_bcs = FieldBoundaryConditions(bottom=z_bc, top=z_bc, south=y_bc, north=y_bc)

    model = NonhydrostaticModel(
               architecture = arch,
                       grid = grid,
                   buoyancy = buoyancy,
                    tracers = (:T, :S),
                    closure = nothing,
        boundary_conditions = (T=T_bcs,)
    )

    T₀(x, y, z) = ∂T∂z * (x*g̃[1] + y*g̃[2] + z*g̃[3])
    set!(model, T=T₀)

    simulation = Simulation(model, Δt=10minute, stop_time=1hour)
    run!(simulation)

    ∂y_T = ComputedField(∂y(model.tracers.T))
    ∂z_T = ComputedField(∂z(model.tracers.T))

    compute!(∂y_T)
    compute!(∂z_T)

    mean_∂y_T = mean(∂y_T)
    mean_∂z_T = mean(∂z_T)

    Δ_y = ∂T∂z * g̃[2] - mean_∂y_T
    Δ_z = ∂T∂z * g̃[3] - mean_∂z_T

    @info "∂T∂z * g̃[2] = $(∂T∂z * g̃[2]), mean(∂y_T) = $(mean_∂y_T), Δ = $Δ_y at t = $(prettytime(model.clock.time)) with θ=$(θ)°"
    @info "∂T∂z * g̃[3] = $(∂T∂z * g̃[3]), mean(∂z_T) = $(mean_∂z_T), Δ = $Δ_z at t = $(prettytime(model.clock.time)) with θ=$(θ)°"

    @test ∂T∂z * g̃[2] ≈ mean(∂y_T)
    @test ∂T∂z * g̃[3] ≈ mean(∂z_T)

    CUDA.@allowscalar begin
        @test all(∂T∂z * g̃[2] .≈ interior(∂y_T))
        @test all(∂T∂z * g̃[3] .≈ interior(∂z_T))
    end

    return nothing
end

function inertial_oscillations_work_with_rotation_in_different_axis(arch, FT)
    grid = RegularRectilinearGrid(FT, size=(), topology=(Flat, Flat, Flat))

    f₀ = 1
    ū = 1
    Δt = 1e-3
    T_inertial = 2π/f₀
    stop_time = T_inertial / 2

    zcoriolis = FPlane(f=f₀)
    xcoriolis = ConstantCartesianCoriolis(f=f₀, rotation_axis=(1,0,0))

    model_x =  NonhydrostaticModel(grid=grid, buoyancy=nothing, tracers=nothing, closure=nothing,
                                   timestepper = :RungeKutta3,
                                   coriolis=xcoriolis,
                                  )
    set!(model_x, v=ū)
    simulation_x = Simulation(model_x, Δt=Δt, stop_time=stop_time)
    run!(simulation_x)

    model_z =  NonhydrostaticModel(grid=grid, buoyancy=nothing, tracers=nothing, closure=nothing,
                                   timestepper = :RungeKutta3,
                                   coriolis=zcoriolis,
                                   )
    set!(model_z, u=ū)
    simulation_z = Simulation(model_z, Δt=Δt, stop_time=stop_time)
    run!(simulation_z)

    u_x = model_x.velocities.u.data[1,1,1]
    v_x = model_x.velocities.v.data[1,1,1]
    w_x = model_x.velocities.w.data[1,1,1]

    u_z = model_z.velocities.u.data[1,1,1]
    v_z = model_z.velocities.v.data[1,1,1]
    w_z = model_z.velocities.w.data[1,1,1]

    @test w_z == 0
    @test u_x == 0

    @test √(v_x^2 + w_x^2) ≈ 1
    @test √(u_z^2 + v_z^2) ≈ 1

    @test u_z ≈ v_x
    @test v_z ≈ w_x

    return nothing
end

timesteppers = (:QuasiAdamsBashforth2, :RungeKutta3)

@testset "Dynamics" begin
    @info "Testing dynamics..."

    @testset "Simple diffusion" begin
        @info "  Testing simple diffusion..."
        for fieldname in (:u, :v, :c), timestepper in timesteppers
            for time_discretization in (ExplicitTimeDiscretization(), VerticallyImplicitTimeDiscretization())
                @test test_diffusion_simple(fieldname, timestepper, time_discretization)
            end
        end
    end

    @testset "Budgets in isotropic diffusion" begin
        @info "  Testing model budgets with isotropic diffusion..."
        for timestepper in timesteppers
            for topology in ((Periodic, Periodic, Periodic),
                             (Periodic, Periodic, Bounded),
                             (Periodic, Bounded, Bounded),
                             (Bounded, Bounded, Bounded))

                if topology !== (Periodic, Periodic, Periodic) # can't use implicit time-stepping in vertically-periodic domains right now
                    time_discretizations = (ExplicitTimeDiscretization(), VerticallyImplicitTimeDiscretization())
                else
                    time_discretizations = (ExplicitTimeDiscretization(),)
                end

                for time_discretization in time_discretizations

                    for closure in (IsotropicDiffusivity(ν=1, κ=1, time_discretization=time_discretization),
                                    AnisotropicDiffusivity(νh=1, νz=1, κh=1, κz=1, time_discretization=time_discretization))

                        fieldnames = [:c]

                        topology[1] === Periodic && push!(fieldnames, :u)
                        topology[2] === Periodic && push!(fieldnames, :v)
                        topology[3] === Periodic && push!(fieldnames, :w)

                        grid = RegularRectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1), topology=topology)

                        model = NonhydrostaticModel(timestepper = timestepper,
                                                           grid = grid,
                                                        closure = closure,
                                                        tracers = :c,
                                                       coriolis = nothing,
                                                       buoyancy = nothing)

                        td = typeof(time_discretization).name.wrapper
                        closurename = typeof(closure).name.wrapper

                        for fieldname in fieldnames
                            @info "    [$timestepper, $td, $closurename] Testing $fieldname budget in a $topology domain with isotropic diffusion..."
                            @test test_isotropic_diffusion_budget(fieldname, model)
                        end
                    end
                end
            end
        end
    end

    @testset "Budgets in biharmonic diffusion" begin
        @info "  Testing model budgets with biharmonic diffusion..."
        for timestepper in timesteppers
            for topology in ((Periodic, Periodic, Periodic),
                             (Periodic, Periodic, Bounded),
                             (Periodic, Bounded, Bounded),
                             (Bounded, Bounded, Bounded))

                fieldnames = [:c]

                topology[1] === Periodic && push!(fieldnames, :u)
                topology[2] === Periodic && push!(fieldnames, :v)
                topology[3] === Periodic && push!(fieldnames, :w)

                grid = RegularRectilinearGrid(size=(2, 2, 2), extent=(1, 1, 1), topology=topology)

                model = NonhydrostaticModel(timestepper = timestepper,
                                                   grid = grid,
                                                closure = AnisotropicBiharmonicDiffusivity(νh=1, νz=1, κh=1, κz=1),
                                               coriolis = nothing,
                                                tracers = :c,
                                               buoyancy = nothing)

                for fieldname in fieldnames
                    @info "    [$timestepper] Testing $fieldname budget in a $topology domain with biharmonic diffusion..."
                    @test test_biharmonic_diffusion_budget(fieldname, model)
                end
            end
        end
    end

    @testset "Diffusion cosine" begin
        for timestepper in (:QuasiAdamsBashforth2,) #timesteppers
            for fieldname in (:u, :v, :T, :S)
                for time_discretization in (ExplicitTimeDiscretization(), VerticallyImplicitTimeDiscretization())
                    Nz, Lz = 128, π/2
                    grid = RegularRectilinearGrid(size=(1, 1, Nz), x=(0, 1), y=(0, 1), z=(0, Lz))

                    @info "  Testing diffusion cosine [$fieldname, $timestepper, $time_discretization]..."
                    @test test_diffusion_cosine(fieldname, timestepper, grid, time_discretization)

                    @info "  Testing diffusion cosine on ImmersedBoundaryGrid [$fieldname, $timestepper, $time_discretization]..."
                    solid(x, y, z) = false
                    immersed_grid = ImmersedBoundaryGrid(grid, GridFittedBoundary(solid))
                    @test test_diffusion_cosine(fieldname, timestepper, immersed_grid, time_discretization)
                end
            end
        end
    end

    @testset "Passive tracer advection" begin
        for timestepper in (:QuasiAdamsBashforth2,) #timesteppers
            @info "  Testing passive tracer advection [$timestepper]..."
            @test passive_tracer_advection_test(timestepper)
        end
    end

    @testset "Internal wave" begin
        include("test_internal_wave_dynamics.jl")

        Nx = Nz = 128
        Lx = Lz = 2π

        # Regular grid with no flat dimension
        y_periodic_regular_grid = RegularRectilinearGrid(topology=(Periodic, Periodic, Bounded),
                                                         size=(Nx, 1, Nz), x=(0, Lx), y=(0, Lx), z=(-Lz, 0))

        # Regular grid with a flat y-dimension
        y_flat_regular_grid = RegularRectilinearGrid(topology=(Periodic, Flat, Bounded),
                                                     size=(Nx, Nz), x=(0, Lx), z=(-Lz, 0))

        # Vertically stretched grid with regular spacing and no flat dimension
        z_faces = collect(znodes(Face, y_periodic_regular_grid))
        y_periodic_regularly_spaced_vertically_stretched_grid = VerticallyStretchedRectilinearGrid(topology=(Periodic, Periodic, Bounded),
                                                                                size=(Nx, 1, Nz), x=(0, Lx), y=(0, Lx), z_faces=z_faces)

        # Vertically stretched grid with regular spacing and no flat dimension
        y_flat_regularly_spaced_vertically_stretched_grid = VerticallyStretchedRectilinearGrid(topology=(Periodic, Flat, Bounded),
                                                                                size=(Nx, Nz), x=(0, Lx), z_faces=z_faces)

        solution, kwargs, background_fields, Δt, σ = internal_wave_solution(L=Lx, background_stratification=false)

        test_grids = (y_periodic_regular_grid,
                      y_flat_regular_grid,
                      y_periodic_regularly_spaced_vertically_stretched_grid,
                      y_flat_regularly_spaced_vertically_stretched_grid)

        @testset "Internal wave with HydrostaticFreeSurfaceModel" begin
            for grid in test_grids
                grid_name = typeof(grid).name.wrapper
                topo = topology(grid)

                # Choose gravitational acceleration so that σ_surface = sqrt(g * Lx) = 10σ
                g = (10σ)^2 / Lx

                model = HydrostaticFreeSurfaceModel(; free_surface=ImplicitFreeSurface(gravitational_acceleration=g), grid=grid, kwargs...)

                @info "  Testing internal wave [HydrostaticFreeSurfaceModel, $grid_name, $topo]..."
                internal_wave_dynamics_test(model, solution, Δt)
            end
        end

        @testset "Internal wave with NonhydrostaticModel" begin
            for grid in test_grids
                grid_name = typeof(grid).name.wrapper
                topo = topology(grid)

                model = NonhydrostaticModel(; grid=grid, kwargs...)

                @info "  Testing internal wave [NonhydrostaticModel, $grid_name, $topo]..."
                internal_wave_dynamics_test(model, solution, Δt)
            end
        end
    end

    @testset "Taylor-Green vortex" begin
        for timestepper in (:QuasiAdamsBashforth2,) #timesteppers
            for time_discretization in (ExplicitTimeDiscretization(), VerticallyImplicitTimeDiscretization())
                td = typeof(time_discretization).name.wrapper
                @info "  Testing Taylor-Green vortex [$timestepper, $td]..."
                @test taylor_green_vortex_test(CPU(), timestepper, time_discretization)
            end
        end
    end

    @testset "Background fields" begin
        for timestepper in (:QuasiAdamsBashforth2,) #timesteppers
            @info "  Testing dynamics with background fields [$timestepper]..."
            @test_skip passive_tracer_advection_test(timestepper, background_velocity_field=true)

            Nx = Nz = 128
            Lx = Lz = 2π

            # Regular grid with no flat dimension
            y_periodic_regular_grid = RegularRectilinearGrid(topology=(Periodic, Periodic, Bounded),
                                                             size=(Nx, 1, Nz), x=(0, Lx), y=(0, Lx), z=(-Lz, 0))

            solution, kwargs, background_fields, Δt, σ = internal_wave_solution(L=Lx, background_stratification=true)

            model = NonhydrostaticModel(; grid=y_periodic_regular_grid, background_fields=background_fields, kwargs...)
            internal_wave_dynamics_test(model, solution, Δt)
        end
    end

    @testset "Tilted gravity" begin
        for arch in archs
            @info "  Testing tilted gravity [$(typeof(arch))]..."
            for θ in (0, 1, -30, 60, 90, -180)
                stratified_fluid_remains_at_rest_with_tilted_gravity_buoyancy_tracer(arch, Float64, θ=θ)
                stratified_fluid_remains_at_rest_with_tilted_gravity_temperature_tracer(arch, Float64, θ=θ)
            end
        end
    end

    @testset "Background rotation about arbitrary axis" begin
        for arch in archs
            @info "  Testing background rotation about arbitrary axis [$(typeof(arch))]..."
            inertial_oscillations_work_with_rotation_in_different_axis(arch, Float64)
        end
    end

end


text:
# https://html.spec.whatwg.org/multipage/server-sent-events.html

module ExampleSSE

using Bukdu

struct SSEController <: ApplicationController
    conn::Conn
end

function index(c::SSEController)
    render(HTML, """
<!DOCTYPE html>

  <meta charset="UTF-8" />
  <title>Server Sent Events</title>

julia> h = first(Bukdu.sse_streams())
HTTP.Streams.Stream{HTTP.Messages.Request,HTTP.ConnectionPool.Transaction{Sockets.TCPSocket}}(HTTP.Messages.Request:
\"\"\"
GET /sse HTTP/1.1
Host: 127.0.0.1:8080
Connection: keep-alive
Accept: text/event-stream
Cache-Control: no-cache
User-Agent: Mozilla/5.0 (Macintosh; Intel Mac OS X 10_13_3) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/65.0.3325.181 Safari/537.36
Referer: http://127.0.0.1:8080/
Accept-Encoding: gzip, deflate, br
Accept-Language: ko-KR,ko;q=0.9,en-US;q=0.8,en;q=0.7

\"\"\", T1  ⏸    1↑🔒    1↓🔒   127.0.0.1:8080:8080 ≣16 inactive 10.4s, true, false, 0)

julia> write(h, string("data: ", repr("Hello"), "\\r\\n\\r\\n"))
23

""")
end

function sse(c::SSEController)
    render(EventStream)
end

end # module ExampleSSE

if PROGRAM_FILE == basename(@__FILE__)

using Bukdu
import .ExampleSSE: SSEController, index, sse

plug(Plug.Logger, access_log=(path=normpath(@__DIR__, "access.log"),), formatter=Plug.LoggerFormatter.datetime_message)

routes() do
    get("/", SSEController, index)
    get("/sse", SSEController, sse)
    plug(Plug.ServerSentEvents) #
end

Bukdu.start(8080)

# Router.call(get, "/") #
# CLI.routes()

Base.JLOptions().isinteractive==0 && wait()

# Bukdu.stop()

end # if

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
export copy_to_shifted!, copy_from_shifted!

function copy_to_shifted!{T}(backend::GPUBackend, dst::CuTensorBlob{T},
    src::CuTensorBlob{T}, shift::NTuple{2,Int})

  dims = size(src)
  dims2 = size(dst)

  shift_dim, shift_amount = shift

  dim_x  = prod(dims[1:shift_dim-1])
  @assert dim_x == prod(dims2[1:shift_dim-1])
  dim_y  = dims[shift_dim]
  dim_y2 = dims2[shift_dim]
  dim_z  = prod(dims[shift_dim+1:end])
  @assert dim_z == prod(dims2[shift_dim+1:end])

  x_block = round(Int, ceil(float(dim_x)/CUDA.THREADS_PER_BLOCK_X))
  y_block = round(Int, ceil(float(dim_y)/CUDA.THREADS_PER_BLOCK_Y))
  z_block = round(Int, ceil(float(dim_z)/CUDA.THREADS_PER_BLOCK_Z))

  if T == Float32
    kernel = backend.mocha.copy_to_shifted_float
  elseif T == Float64
    kernel = backend.mocha.copy_to_shifted_double
  else
    error("Unsupported data type $T for shifted copy")
  end

  CUDA.launch(kernel, (x_block,y_block,z_block),
      (CUDA.THREADS_PER_BLOCK_X,CUDA.THREADS_PER_BLOCK_Y,CUDA.THREADS_PER_BLOCK_Z),
      (dst.ptr.p,src.ptr.p,dim_x,dim_y,dim_z,dim_y2,shift_amount))
end

function copy_from_shifted!{T}(backend::GPUBackend, dst::CuTensorBlob{T},
    src::CuTensorBlob{T}, shift::NTuple{2,Int})

  dims = size(dst)
  dims2 = size(src)

  shift_dim, shift_amount = shift

  dim_x  = prod(dims[1:shift_dim-1])
  @assert dim_x == prod(dims2[1:shift_dim-1])
  dim_y  = dims[shift_dim]
  dim_y2 = dims2[shift_dim]
  dim_z  = prod(dims[shift_dim+1:end])
  @assert dim_z == prod(dims2[shift_dim+1:end])

  x_block = round(Int, ceil(float(dim_x)/CUDA.THREADS_PER_BLOCK_X))
  y_block = round(Int, ceil(float(dim_y)/CUDA.THREADS_PER_BLOCK_Y))
  z_block = round(Int, ceil(float(dim_z)/CUDA.THREADS_PER_BLOCK_Z))

  if T == Float32
    kernel = backend.mocha.copy_from_shifted_float
  elseif T == Float64
    kernel = backend.mocha.copy_from_shifted_double
  else
    error("Unsupported data type $T for shifted copy")
  end

  CUDA.launch(kernel, (x_block,y_block,z_block),
      (CUDA.THREADS_PER_BLOCK_X,CUDA.THREADS_PER_BLOCK_Y,CUDA.THREADS_PER_BLOCK_Z),
      (dst.ptr.p,src.ptr.p,dim_x,dim_y,dim_z,dim_y2,shift_amount))
end


text:
# functions relating to more general input to the generalized bunching conjecture
# instead of pure state inputs

"""

    schur_matrix(H)

computes the Schur matrix as defined in Eq. 1 of Linear Algebra and its Applications 490 (2016) 196–201
"""
function schur_matrix(H)

    @test is_positive_semidefinite(H)

    n = size(H,1)
    schur_mat = Matrix{eltype(H)}(undef,(factorial(n), factorial(n)))

    @showprogress for (i,sigma) in enumerate(permutations(collect(1:n)))
        for (j,rho) in enumerate(permutations(collect(1:n)))
            schur_mat[i,j] = prod([H[sigma[k], rho[k]] for k in 1:n])
        end
    end

    schur_mat

end

function find_permutation_index(this_perm, permutation_array)
    findfirst(x->x == this_perm, permutation_array)
end

function multiply_permutations(a,b)
    a[b[:]]
end

"""

    J_array(theta, n)

returns the `J` as defined in in Eq.10 of [Universality of Generalized Bunching and Efficient Assessment of Boson Sampling](https://arxiv.org/abs/1509.01561), with the permutations coming in the order given by `permutations(collect(1:n))`
"""
function J_array(theta, n)

    J = zeros(eltype(theta), factorial(n))

    permutation_array = collect(permutations(collect(1:n)))

    @showprogress for (i,sigma) in enumerate(permutations(collect(1:n)))
        for (j,tau) in enumerate(permutations(collect(1:n)))
            J[i] += conj(theta[j]) * theta[find_permutation_index(multiply_permutations(tau, sigma), permutation_array)]
        end
    end

    J

end

"""
    density_matrix_from_J(J,n)

density matrix associated to a J function as defined in [Universality of Generalized Bunching and Efficient Assessment of Boson Sampling](https://arxiv.org/abs/1509.01561), computed through Eq. 46.
"""
function density_matrix_from_J(J,n)

    density_matrix = zeros(eltype(J), factorial(n), factorial(n))
    permutation_array = collect(permutations(collect(1:n)))

    @showprogress for (i,sigma) in enumerate(permutations(collect(1:n)))
        for (j,tau) in enumerate(permutations(collect(1:n)))
            density_matrix[j,find_permutation_index(multiply_permutations(sigma, tau), permutation_array)] += J[i]
        end
    end

    density_matrix/factorial(n)

end


text:
using Zygote
using BenchmarkTools

fn_scalar(a) = 3a ^ 10 - 5
fn_vector(a) = sum(3 .* a .^ 10 .- 5)
fn_vector_vector(a) = sum(a .* a)
fn_matrix(a) = sum(3 .* a .^ 10 .- 5)
fn_matrix_matrix(a) = sum(a * a')

val_scalar = 100.
val_vector = rand(100)
val_matrix = rand(10, 10)

# Precompile
fn_scalar(val_scalar)
fn_vector(val_vector)
fn_vector_vector(val_vector)
fn_matrix(val_matrix)
fn_matrix_matrix(val_matrix)

gradient(fn_scalar, val_scalar)
gradient(fn_vector, val_vector)
gradient(fn_vector_vector, val_vector)
gradient(fn_matrix, val_matrix)
gradient(fn_matrix_matrix, val_matrix)

# Benchmark
res_1 = @benchmark gradient(fn_scalar, val_scalar)
res_2 = @benchmark gradient(fn_vector, val_vector)
res_3 = @benchmark gradient(fn_vector_vector, val_vector)
res_4 = @benchmark gradient(fn_matrix, val_matrix)
res_5 = @benchmark gradient(fn_matrix_matrix, val_matrix)

open("zygote.txt", "w") do file
    write(file, string("Scalar: ", res_1, "\n"))
    write(file, string("Vector: ", res_2, "\n"))
    write(file, string("Vector Vector: ", res_3, "\n"))
    write(file, string("Matrix: ", res_4, "\n"))
    write(file, string("Matrix Matrix: ", res_5, "\n"))
end

max_stars_repo_name: Matthew-Amos/ad-benchmarks

max_issues_repo_name: Matthew-Amos/ad-benchmarks

max_forks_repo_name: Matthew-Amos/ad-benchmarks


text:
"""
    psaverage(psysc) -> sys::DescriptorStateSpace

Compute for the continuous-time periodic system `psysc = (A(t),B(t),C(t),D(t))`
the corresponding time averaged LTI system `sys = (Am,Bm,Cm,Dm)` over one period.
"""
function psaverage(psysc::PeriodicStateSpace{PM}) where {T,PM <: AbstractPeriodicArray{:c,T}}
    return dss(pmaverage(psysc.A), pmaverage(psysc.B), pmaverage(psysc.C), pmaverage(psysc.D))
end
"""
    psteval(psysc,tval) -> sys::DescriptorStateSpace

Compute for the continuous-time periodic system `psysc = (A(t),B(t),C(t),D(t))` and a time value `tval`,
the LTI system `sys = (A(tval),B(tval),C(tval),D(tval))`.
"""
function psteval(psysc::PeriodicStateSpace{PM}, tval::Real = 0) where {T,PM <: AbstractPeriodicArray{:c,T}}
    return dss(tvmeval(psysc.A, tval)[1], tvmeval(psysc.B, tval)[1], tvmeval(psysc.C, tval)[1], tvmeval(psysc.D, tval)[1])
end
"""
     psmrc2d(sys, Ts; ki, ko) -> psys::PeriodicStateSpace{PeriodicMatrix}

Perform the multirate discretization of a linear time-invariant system.

For a continuous-time state-space system `sys = (A,B,C,D)`, a basic sampling time `Ts`,
and the integer vectors `ki` and `ko` containing, respectively, the numbers of
input and output sampling subperiods, the corresponding discretized
periodic system `psys = (Ap,Bp,Cp,Dp)` of period `T = n*Ts` is determined,
where `n` is the least common multiple
of the integers contained in `ki` and `ko`. For a continuous-time system `sys`
a zero-order hold based discretization method is used, such that the `i`-th input
is constant during intervals of length `ki[i]*Ts`.
An output hold device is used to provide constant intersample outputs, such that the `i`-th output
is constant during intervals of length `ko[i]*Ts`.
For a discrete-time system with a defined sample time `Δ`,
an input and output resampling is performed using `Δ` as basic
sample time and the value of `Ts` is ignored.
If the system sample time is not defined, then the value of
`Ts` is used as the basic sample time.
"""
function psmrc2d(sys::DescriptorStateSpace{T}, Ts::Real; ki::Vector{Int} = ones(Int,size(sys,2)), ko::Vector{Int}= ones(Int,size(sys,1))) where T
    Ts > 0 || error("the sampling time Ts must be positive")
    (all(ki .>= 1) && all(ko .>= 1)) || error("the numbers of subperiods must be at least one")
    p, m = size(sys)
    length(ki) == m || error("the number of input components of ki must be $m")
    length(ko) == p || error("the number of ouput components of ko must be $p")
    A, E, B, C, D = dssdata(sys)
    E == I || error("only standard state-space systems supported")
    Δ = sys.Ts
    T1 = T <: BlasFloat ? T : promote_type(Float64,T)
    n = sys.nx
    if Δ != 0
       i1 = 1:n; i2 = n+1:n+m
       Δ > 0 && (Ts = Δ)
       G = exp([ rmul!(A,Ts) rmul!(B,Ts); zeros(T1,m,n+m)])
       A =  G[i1,i1]
       B =  G[i1,i2]
    end

    K = lcm(lcm(ki),lcm(ko))
    Ap = similar(Vector{Matrix},K)
    Bp = similar(Vector{Matrix},K)
    Cp = similar(Vector{Matrix},K)
    Dp = similar(Vector{Matrix},K)
    for i = 1:K
        si = mod.(i-1,ki) .== 0
        Si = Diagonal(si)
        Sineg = Diagonal(.!si)
        Ti = Diagonal(mod.(i-1,ko) .== 0)
        Ap[i] = [A  B*Sineg; zeros(T1,m,n) Sineg]
        Bp[i] = [B*Si; Si]
        Cp[i] = [Ti*C Ti*D*Sineg]
        Dp[i] = Ti*D
    end
    PMT = PeriodicMatrix
    period = K*Ts
    return ps(PMT(Ap,period),PMT(Bp,period),PMT(Cp,period),PMT(Dp,period))
end
"""
     psc2d(psysc, Ts; solver, reltol, abstol, dt) -> psys::PeriodicStateSpace{PeriodicMatrix}

Compute for the continuous-time periodic system `psysc = (A(t),B(t),C(t),D(t))` of period `T` and
for a sampling time `Ts`, the corresponding discretized
periodic system `psys = (Ad,Bd,Cd,Dd)` using a zero-order hold based discretization method.

The discretization is performed by determining the monodromy matrix as a product of
`K = T/Ts` state transition matrices of the extended state-space matrix [A(t) B(t); 0 0]
by integrating numerically the corresponding homogeneous linear ODE.
The ODE solver to be employed can be
specified using the keyword argument `solver`, together with
the required relative accuracy `reltol` (default: `reltol = 1.e-3`),
absolute accuracy `abstol` (default: `abstol = 1.e-7`) and/or
the fixed step length `dt` (default: `dt = Ts/10`) (see [`tvstm`](@ref)).
For large values of `K`, parallel computation of factors can be alternatively performed
by starting Julia with several execution threads.
The number of execution threads is controlled either by using the `-t/--threads` command line argument
or by using the `JULIA_NUM_THREADS` environment variable.
"""
function psc2d(psysc::PeriodicStateSpace{PM}, Ts::Real; solver::String  = "", reltol = 1e-3, abstol = 1e-7, dt = Ts/10) where {T,PM <: AbstractPeriodicArray{:c,T}}

    Ts > 0 || error("the sampling time Ts must be positive")
    period = psysc.period
    r = rationalize(period/Ts)
    denominator(r) == 1 || error("incommensurate period and sample time")
    K = numerator(r)

    T1 = T <: BlasFloat ? T : (T <: Num ? Float64 : promote_type(Float64,T))
    ONE = one(T1)
    n, m = size(psysc.B)

    # quick exit in constant case
    islti(psysc) && (return ps(c2d(psaverage(psysc),Ts)[1],psysc.period))

    ka, na = isconstant(psysc.A) ? (1,K) : (K,1)
    kb, nb = isconstant(psysc.A) && isconstant(psysc.B) ? (1,K) : (K,1)
    kc, nc = isconstant(psysc.C) ? (1,K) : (K,1)
    kd, nd = isconstant(psysc.D) ? (1,K) : (K,1)
    Ap = similar(Vector{Matrix},ka)
    Bp = similar(Vector{Matrix},kb)
    Cp = similar(Vector{Matrix},kc)
    Dp = similar(Vector{Matrix},kd)

    psys1 = typeof(psysc.A) <: PeriodicFunctionMatrix ? psysc :
               convert(PeriodicStateSpace{PeriodicFunctionMatrix},psysc)
    i1 = 1:n; i2 = n+1:n+m

    if isconstant(psysc.A) && isconstant(psysc.B)
        G = exp([ rmul!(psys1.A.f(0),Ts) rmul!(psys1.B.f(0),Ts); zeros(T1,m,n+m)])
        Ap[1] = G[i1,i1]
        Bp[1] = G[i1,i2]
        ts = zero(T1)
        for i = 1:kc
            Cp[i] = psys1.C.f(ts)
            ts += Ts
        end
        ts = zero(T1)
        for i = 1:kd
            Dp[i] = psys1.D.f(ts)
            ts += Ts
        end
    else
        Gfun = PeriodicFunctionMatrix(t -> [psys1.A.f(t) psys1.B.f(t); zeros(T1,m,n+m)], period)
        G = monodromy(Gfun, K; solver, reltol, abstol, dt)
        [Ap[i] = G.M[i1,i1,i] for i in 1:ka]
        [Bp[i] = G.M[i1,i2,i] for i in 1:kb]
        ts = zero(T1)
        for i = 1:kc
            Cp[i] = psys1.C.f(ts)
            ts += Ts
        end
        ts = zero(T1)
        for i = 1:kd
            Dp[i] = psys1.D.f(ts)
            ts += Ts
        end
    end
    PMT = PeriodicMatrix
    return ps(PMT(Ap,period; nperiod = na),PMT(Bp,period; nperiod = nb),PMT(Cp,period; nperiod = nc),PMT(Dp,period; nperiod = nd))
    # end PSC2D
end


text:
# Generate low-level api from METADATA

using MsgPack
import MsgPack: Extension, ExtensionType

# Types defined by the api
struct NvimApiObject{N} <: NvimObject
    client::NvimClient
    # TODO: use a fixarray or UInt64
    hnd::Vector{UInt8}
end

import Base.==
==(a::NvimApiObject{N}, b::NvimApiObject{N}) where {N} = a.hnd == b.hnd

NvimApiObject(c, e::Extension) = NvimApiObject{UInt8(e.type)}(c, e.data)

# Not really module-interface clean, I know...
# TODO(smolck): Was just s with Any type, is this right?
MsgPack.pack(s::IO, o::NvimApiObject{N}) where {N} = MsgPack.pack(s, Extension(N, o.hnd))

symbolize(val::Dict) = Dict{Symbol,Any}([(symbolize(k), symbolize(v)) for (k, v) in val])
symbolize(val::Vector{UInt8}) = Symbol(String(val))
symbolize(val::String) = Symbol(val)
symbolize(val::Vector) = [symbolize(v) for v in val]
symbolize(val) = val

function _get_metadata()
    data = read(`nvim --api-info`, String)
    return symbolize(unpack(data))
end

const _metadata = _get_metadata()
const _types = _metadata[:types]
const _functions = _metadata[:functions]

# will break if the api starts using overloading
const api_methods = [f[:name] => f for f in _functions]

const Bytes = Union{String,Vector{UInt8}}
const typemap = Dict{Symbol,Type}(
    :Integer => Integer,
    :Boolean => Bool,
    :String => Bytes,
    :Array => Vector{Any},
    :Dictionary => Dict,
    :Object => Any,
    :void => Nothing,
)

for (name, info) in _types
    id = info[:id]
    @eval begin
        const $(name) = NvimApiObject{$id}
        # typealias $(name) NvimApiObject{$id}
        typemap[$(Meta.quot(name))] = $name
    end
end

function parsetype(t::Symbol)
    components = split(rstrip(string(t), ')'), '(')
    top = Symbol(components[1])
    params = length(components) > 1 ? split(components[2], ", ") : []
    (top, params)
end

function gettype(t::Symbol)
    top, params = parsetype(t)
    if haskey(typemap, top)
        typemap[top]
    elseif top == :ArrayOf
        tcontains = gettype(Symbol(params[1]))
        if length(params) == 1
            Vector{tcontains}
        else
            Tuple{fill(tcontains, parse(Int, params[2]))...}
        end
    else
        Any
    end
end

# fake covariant Vectors
checkarg(::Type{Vector{T}}, val::Vector{T}) where {T} = val
checkarg(::Type{Vector{T}}, val::Vector) where {T} = convert(Vector{T}, val)

# I'm pretty sure julia somehow can express subdiagnal dispatch,
# but implement the cases manually for now
checkarg(::Type{Vector{Bytes}}, val::Vector{T}) where {T <: Bytes} = val
checkarg(::Type{Vector{Integer}}, val::Vector{T}) where {T <: Integer} = val

# this must be specialcased though
checkarg(::Type{Vector{Integer}}, val::Vector{UInt8}) = Int[val...]

# tuples are covariant, so no problem
checkarg(::Type{Vector{T}}, val::(Tuple{Vararg{T}})) where {T} = [val...]

retconvert(typ::Union{Type{Any},Type{Bytes}}, c, val::Union{String,Vector{UInt8}}) =
    String(val)

# needed for disambiguation
retconvert(typ::Type{Any}, c, val::Vector{UInt8}) = String(val)
retconvert(typ::Union{Type{Any},Type{Bool}}, c, val::Bool) = val

# this assumes the current unpack implementation in MsgPack,
# where all int types get promoted to 64 bit
retconvert(typ::Union{Type{Any},Type{Integer}}, c, val::Int64) = val
retconvert(typ::Union{Type{Any},Type{Nothing}}, c, val::Nothing) = nothing
retconvert(typ::Union{Type{Any},Type{Dict}}, c, val::Dict) =
    Dict{String,Any}([(String(k), retconvert(Any, c, v)) for (k, v) in val])

# we assume Msgpack only generates untyped arrays
retconvert(typ::Type{Vector{T}}, c, val::Vector) where {T} =
    [retconvert(T, c, v) for v in val]
retconvert(typ::Type{Any}, c, val::Vector) = Any[retconvert(Any, c, v) for v in val]

# this is not really strict enough
# retconvert{T<:Tuple}(typ::Type{T}, c, val::Vector) = tuple(val...)::T
# handle the only case, 2-tuples, manually for now
retconvert(::Type{Tuple{T,U}}, c, val::Vector) where {T,U} =
    (retconvert(T, c, val[1]), retconvert(U, c, val[2]))

retconvert(typ::Type{NvimApiObject{N}}, c, val::Extension) where {N} = NvimApiObject(c, val)::NvimApiObject{N}

retconvert(typ::Type{T}, c, val::T) where {T} = val

# a stagedfunction will probably be simpler and better
function build_function(f)
    name = f[:name]
    params = f[:parameters]
    tret = f[:return_type]

    parts = split(string(name), "_"; limit=2)
    receiver = parts[1]
    shortname = Symbol(parts[2])
    if shortname == :eval
        shortname = :vim_eval
    end

    body = Expr[]
    args = Symbol[Symbol(string("a_", p[2])) for p in params]
    j_args = Expr[]

    if length(args) != 0
        for (i, p) in enumerate(params)
            # this is probably too restrictive sometimes,
            # use convert for some types (sequences)?
            t = gettype(p[1])
            # if type is an array
            # we will allow any Vector argument
            # and dynamically check if not a subtype
            erased = t
            if t <: Vector
                erased = Union{Vector,Tuple{Vararg{eltype(t)}}}
                push!(body, :($(args[i]) = checkarg($t, $(args[i]))))
            end
            push!(j_args, :($(args[i])::($erased)))
        end
    end

    if length(j_args) == 0 || receiver == "vim"
        pushfirst!(j_args, :(c::NvimClient))
    else
        pushfirst!(body, :(c = ($(args[1])).client))
    end

    # when array constructor non-concatenating, we could drop the Any
    arglist = :(Any[])
    append!(arglist.args, args)
    push!(body, :(res = send_request(c, $(Meta.quot(name)), $arglist)))

    # "representative" type, not neccessarily the one
    # retconvert will assert
    typ = gettype(tret)
    push!(body, :(retconvert($typ, c, res)))

    j_call = Expr(:call, shortname, j_args...)
    fun = Expr(:function, j_call, Expr(:block, body...))
    fun
end

# TODO: maybe this should be in a submodule
# then one could do `importall Neovim.API`
for f in _functions
    eval(build_function(f))
end


text:
using Test
using HTTP
using HTTP.IOExtras, HTTP.Sockets, HTTP.WebSockets
using Sockets

@testset "WebSockets" begin
    p = 8085 # rand(8000:8999)
    socket_type = ["wss", "ws"]

    function listen_localhost()
        @async HTTP.listen(Sockets.localhost, p) do http
            if WebSockets.is_upgrade(http.message)
                WebSockets.upgrade(http) do ws
                    while !eof(ws)
                        data = readavailable(ws)
                        write(ws, data)
                    end
                end
            end
        end
    end

    @testset "External Host - $s" for s in socket_type
        WebSockets.open("$s://echo.websocket.org") do ws
            write(ws, "Foo")
            @test !eof(ws)
            @test String(readavailable(ws)) == "Foo"

            write(ws, "Foo"," Bar")
            @test !eof(ws)
            @test String(readavailable(ws)) == "Foo Bar"

            # send fragmented message manually with ping in between frames
            WebSockets.wswrite(ws, ws.frame_type, "Hello ")
            WebSockets.wswrite(ws, WebSockets.WS_FINAL | WebSockets.WS_PING, "things")
            WebSockets.wswrite(ws, WebSockets.WS_FINAL, "again!")
            @test String(readavailable(ws)) == "Hello again!"

            write(ws, "Hello")
            write(ws, " There")
            write(ws, " World", "!")
            closewrite(ws)

            buf = IOBuffer()
            write(buf, ws)
            @test String(take!(buf)) == "Hello There World!"
        end
    end

    @testset "Localhost" begin
        listen_localhost()

        WebSockets.open("ws://127.0.0.1:$(p)") do ws
            write(ws, "Foo")
            @test String(readavailable(ws)) == "Foo"

            write(ws, "Bar")
            @test String(readavailable(ws)) == "Bar"

            write(ws, "This", " is", " a", " fragmented", " message.")
            @test String(readavailable(ws)) == "This is a fragmented message."

            # send fragmented message manually with ping in between frames
            WebSockets.wswrite(ws, ws.frame_type, "Ping ")
            WebSockets.wswrite(ws, WebSockets.WS_FINAL | WebSockets.WS_PING, "stuff")
            WebSockets.wswrite(ws, WebSockets.WS_FINAL, "pong!")
            @test String(readavailable(ws)) == "Ping pong!"
        end
    end

    @testset "Extended feature support for listen" begin
        port=UInt16(8086)
        tcpserver = listen(port)
        target = "/query?k1=v1&k2=v2"

        servertask =  @async WebSockets.listen("127.0.0.1", port; server=tcpserver) do ws
            @test ws.request isa HTTP.Request
            write(ws, ws.request.target)
            while !eof(ws)
                write(ws, readavailable(ws))
            end
            close(ws)
        end

        WebSockets.open("ws://127.0.0.1:$(port)$(target)") do ws
            @test String(readavailable(ws)) == target
            @test write(ws, "Bye!") == 4
            @test String(readavailable(ws)) == "Bye!"
            close(ws)
        end

        close(tcpserver)
        @test timedwait(()->servertask.state === :failed, 5.0) === :ok
        @test_throws Exception wait(servertask)
    end
end


text:
##################################################################
##################### Objective with candidate storage
##################################################################

function objective_min_cost_storage(pm::_PM.AbstractPowerModel)
    investment = sum(
        calc_convdc_ne_cost(pm, n)
        + calc_ne_branch_cost(pm, n)
        + calc_branchdc_ne_cost(pm, n)
        + calc_ne_storage_cost(pm, n)
        for n in nw_ids(pm; hour=1)
    )
    operation = sum(
        calc_gen_cost(pm, n)
        for n in nw_ids(pm)
    )
    JuMP.@objective(pm.model, Min, investment + operation)
end

#################################################################
##################### Objective with candidate storage and flexible demand
##################################################################

function objective_min_cost_flex(pm::_PM.AbstractPowerModel)
    investment = sum(
        calc_convdc_ne_cost(pm, n)
        + calc_ne_branch_cost(pm, n)
        + calc_branchdc_ne_cost(pm, n)
        + calc_ne_storage_cost(pm, n)
        + calc_load_investment_cost(pm, n)
        for n in nw_ids(pm; hour=1)
    )
    operation = sum(
        calc_gen_cost(pm, n)
        + calc_load_operational_cost(pm, n)
        for n in nw_ids(pm)
    )
    JuMP.@objective(pm.model, Min, investment + operation)
end

function objective_min_cost_flex(t_pm::_PM.AbstractPowerModel, d_pm::_PM.AbstractPowerModel)
    t_investment = sum(
        calc_convdc_ne_cost(t_pm, n)
        + calc_ne_branch_cost(t_pm, n)
        + calc_branchdc_ne_cost(t_pm, n)
        + calc_ne_storage_cost(t_pm, n)
        + calc_load_investment_cost(t_pm, n)
        for n in nw_ids(t_pm; hour=1)
    )
    t_operation = sum(
        calc_gen_cost(t_pm, n)
        + calc_load_operational_cost(t_pm, n)
        for n in nw_ids(t_pm)
    )
    d_investment = sum( # Note: distribution networks do not have DC components (modeling decision)
        calc_ne_branch_cost(d_pm, n)
        + calc_ne_storage_cost(d_pm, n)
        + calc_load_investment_cost(d_pm, n)
        for n in nw_ids(d_pm; hour=1)
    )
    d_operation = sum(
        calc_gen_cost(d_pm, n)
        + calc_load_operational_cost(d_pm, n)
        for n in nw_ids(d_pm)
    )
    JuMP.@objective(t_pm.model, Min, t_investment + t_operation + d_investment + d_operation) # Note: t_pm.model == d_pm.model
end

##########################################################################
##################### Stochastic objective with storage & flex candidates
##########################################################################

function objective_stoch_flex(pm::_PM.AbstractPowerModel)
    investment = sum(
        calc_convdc_ne_cost(pm, n)
        + calc_ne_branch_cost(pm, n)
        + calc_branchdc_ne_cost(pm, n)
        + calc_ne_storage_cost(pm, n)
        + calc_load_investment_cost(pm, n)
        for n in nw_ids(pm; hour=1, scenario=1)
    )
    operation = sum(scenario["probability"] *
        sum(
            calc_gen_cost(pm, n)
            + calc_load_operational_cost(pm, n)
            for n in nw_ids(pm; scenario=s)
        )
        for (s, scenario) in dim_prop(pm, :scenario)
    )
    JuMP.@objective(pm.model, Min, investment + operation)
end

##########################################################################
##################### Auxiliary functions
##########################################################################

function calc_gen_cost(pm::_PM.AbstractPowerModel, n::Int)

    function calc_single_gen_cost(i, g_cost)
        len = length(g_cost)
        cost = 0.0
        if len >= 1
            cost = g_cost[len] # Constant term
            if len >= 2
                cost += g_cost[len-1] * _PM.var(pm,n,:pg,i) # Adds linear term
            end
        end
        return cost
    end

    gen = _PM.ref(pm, n, :gen)
    cost = sum(calc_single_gen_cost(i,g["cost"]) for (i,g) in gen)
    if get(pm.setting, "add_co2_cost", false)
        cost += sum(g["emission_factor"] * _PM.var(pm,n,:pg,i) * pm.ref[:it][_PM.pm_it_sym][:co2_emission_cost] for (i,g) in gen)
    end
    ndgen = _PM.ref(pm, n, :ndgen)
    if !isempty(ndgen)
        cost += sum(g["cost_curt"] * _PM.var(pm,n,:pgcurt,i) for (i,g) in ndgen)
    end
    return cost
end

function calc_convdc_ne_cost(pm::_PM.AbstractPowerModel, n::Int)
    cost = 0.0
    if haskey(_PM.ref(pm, n), :convdc_ne)
        convdc_ne = _PM.ref(pm, n, :convdc_ne)
        if !isempty(convdc_ne)
            cost = sum(conv["cost"]*_PM.var(pm, n, :conv_ne_investment, i) for (i,conv) in convdc_ne)
            if get(pm.setting, "add_co2_cost", false)
                cost += sum(conv["co2_cost"]*_PM.var(pm, n, :conv_ne_investment, i) for (i,conv) in convdc_ne)
            end
        end
    end
    return cost
end

function calc_ne_branch_cost(pm::_PM.AbstractPowerModel, n::Int)
    cost = 0.0
    if haskey(_PM.ref(pm, n), :ne_branch)
        ne_branch = _PM.ref(pm, n, :ne_branch)
        if !isempty(ne_branch)
            cost = sum(branch["construction_cost"]*_PM.var(pm, n, :branch_ne_investment, i) for (i,branch) in ne_branch)
            if get(pm.setting, "add_co2_cost", false)
                cost += sum(branch["co2_cost"]*_PM.var(pm, n, :branch_ne_investment, i) for (i,branch) in ne_branch)
            end
        end
    end
    return cost
end

function calc_branchdc_ne_cost(pm::_PM.AbstractPowerModel, n::Int)
    cost = 0.0
    if haskey(_PM.ref(pm, n), :branchdc_ne)
        branchdc_ne = _PM.ref(pm, n, :branchdc_ne)
        if !isempty(branchdc_ne)
            cost = sum(branch["cost"]*_PM.var(pm, n, :branchdc_ne_investment, i) for (i,branch) in branchdc_ne)
            if get(pm.setting, "add_co2_cost", false)
                cost += sum(branch["co2_cost"]*_PM.var(pm, n, :branchdc_ne_investment, i) for (i,branch) in branchdc_ne)
            end
        end
    end
    return cost
end

function calc_ne_storage_cost(pm::_PM.AbstractPowerModel, n::Int)
    cost = 0.0
    if haskey(_PM.ref(pm, n), :ne_storage)
        ne_storage = _PM.ref(pm, n, :ne_storage)
        if !isempty(ne_storage)
            cost = sum((storage["eq_cost"] + storage["inst_cost"])*_PM.var(pm, n, :z_strg_ne_investment, i) for (i,storage) in ne_storage)
            if get(pm.setting, "add_co2_cost", false)
                cost += sum(storage["co2_cost"]*_PM.var(pm, n, :z_strg_ne_investment, i) for (i,storage) in ne_storage)
            end
        end
    end
    return cost
end

function calc_load_cost(pm::_PM.AbstractPowerModel, n::Int)
    return calc_load_operational_cost(pm, n) + calc_load_investment_cost(pm, n)
end

function calc_load_operational_cost(pm::_PM.AbstractPowerModel, n::Int)
    cost = 0.0
    flex_load = _PM.ref(pm, n, :flex_load)
    if !isempty(flex_load)
        cost += sum(
            l["cost_shift"]*_PM.var(pm, n, :pshift_up, i) # Do not add `:pshift_down`: would result in double counting.
            + l["cost_red"]*_PM.var(pm, n, :pred, i)
            + l["cost_curt"]*_PM.var(pm, n, :pcurt, i)
            for (i,l) in flex_load
        )
    end
    fixed_load = _PM.ref(pm, n, :fixed_load)
    if !isempty(fixed_load)
        cost += sum(
            l["cost_curt"]*_PM.var(pm, n, :pcurt, i)
            for (i,l) in fixed_load
        )
    end
    return cost
end

function calc_load_investment_cost(pm::_PM.AbstractPowerModel, n::Int)
    cost = 0.0
    flex_load = _PM.ref(pm, n, :flex_load)
    if !isempty(flex_load)
        cost = sum(l["cost_inv"]*_PM.var(pm, n, :z_flex_investment, i) for (i,l) in flex_load)
        if get(pm.setting, "add_co2_cost", false)
            cost += sum(l["co2_cost"]*_PM.var(pm, n, :z_flex_investment, i) for (i,l) in flex_load)
        end
    end
    return cost
end


text:
using PyPlot
using LinearAlgebra
g1 = 1/2-√3/6
g2 = 1/2+√3/6

function calc_basis(ξ, η, x1=1.0, y1=0.0, x2=0.0, y2=1.0, x3=0.0, y3=0.0)
    λ = [ξ;η;1-ξ-η]
    K = 0.5 * abs(det([
        x2-x1 x3-x1
        y2-y1 y3-y1
    ]))
    ∇λ = [
        y2-y3 x3-x2
        y3-y1 x1-x3
        y1-y2 x2-x1
    ]/2K
    ∇λ = [∇λ[:,2] -∇λ[:,1]]
    return λ, ∇λ
end

function φ1(ξ, η)
    λ, ∇λ = calc_basis(ξ, η)
    return λ[1]*∇λ[2,:]
end

function φ2(ξ, η)
    λ, ∇λ = calc_basis(ξ, η)
    return λ[2]*∇λ[1,:]
end

function φ3(ξ, η)
    λ, ∇λ = calc_basis(ξ, η)
    return λ[2]*∇λ[3,:]
end

function φ4(ξ, η)
    λ, ∇λ = calc_basis(ξ, η)
    return λ[3]*∇λ[2,:]
end

function φ5(ξ, η)
    λ, ∇λ = calc_basis(ξ, η)
    return λ[3]*∇λ[1,:]
end

function φ6(ξ, η)
    λ, ∇λ = calc_basis(ξ, η)
    return λ[1]*∇λ[3,:]
end

Basis = [φ1, φ2, φ3, φ4, φ5, φ6]

a = LinRange(0, 1, 10)
ξ = zeros(length(a)^2)
η = zeros(length(a)^2)
k = 0
for i = 1:length(a)
    for j = 1:length(a)
        global k = k + 1
        ξ[k] = a[i]
        η[k] = a[j]
    end
end
Idx = ξ + η .<= 1
ξ = ξ[Idx]
η = η[Idx]

m = Dict(1=>1, 2=>4, 3=>2, 4=>5, 5=>3, 6=>6)
figure(figsize=(15, 10))
for k = 1:6
    subplot(230+m[k])
    U = zeros(length(ξ))
    V = zeros(length(ξ))
    for i = 1:length(ξ)
        U[i], V[i] = Basis[k](ξ[i], η[i])
    end
    plot([0;1;0;0], [0;0;1;0])
    quiver(ξ, η, U, V)
    title("\$\\phi_$k\$")
    axis("equal")
    xlim(-0.25,1.25)
    ylim(-0.25,1.25)
    if m[k] in [4,5,6]
        xlabel("\$\\xi\$")
    end
    if m[k] in [1,4]
        ylabel("\$\\eta\$")
    end
end
savefig("BDM.png")


text:
using Primes

@testset "Rational reconstruction" begin

    ms = BigInt[107, 199, 509, 2^31 - 1, nextprime(BigInt(2)^100), nextprime(BigInt(2)^200)]
    as = [QQ(0), QQ(1), QQ(1, 2), QQ(3, 2), QQ(4, 5), QQ(1, 7), QQ(6)]

    num, den, buf, buf1, buf2, buf3, u1, u2, u3, v1, v2, v3 = [BigInt(0) for _ in 1:13]
    for m in ms
        bnd = Groebner.rational_reconstruction_bound(m)
        for a in as
            ac = numerator(a) * invmod(denominator(a), m)
            ar1 = Groebner.rational_reconstruction!(
                        num, den, bnd, buf, buf1, buf2, buf3,
                        u1, u2, u3, v1, v2, v3, ac, m)
            ar2 = Groebner.rational_reconstruction(ac, m)
            @test Base.unsafe_rational(num, den) == ar2 == a
        end
    end

    # random tests
    ms = Primes.nextprimes(BigInt(2^31-1), 100)
    ases = 100
    for m in ms
        bnd = Groebner.rational_reconstruction_bound(m)
        for i in 1:ases
            ac = rand(0:m-1)
            ar1 = Groebner.rational_reconstruction(ac, m)
            ar2 = Groebner.rational_reconstruction!(
                        num, den, bnd, buf, buf1, buf2, buf3,
                        u1, u2, u3, v1, v2, v3, ac, m)
            @test Base.unsafe_rational(num, den) == ar1
        end
    end

    # Big tests
    ms = Primes.nextprimes(BigInt(2)^100, 10)
    ases = 100
    for m in ms
        bnd = Groebner.rational_reconstruction_bound(m)
        for i in 1:ases
            ac = rand(0:m-1)
            ar1 = Groebner.rational_reconstruction(ac, m)
            ar2 = Groebner.rational_reconstruction!(
                        num, den, bnd, buf, buf1, buf2, buf3,
                        u1, u2, u3, v1, v2, v3, ac, m)
            @test Base.unsafe_rational(num, den) == ar1
        end
    end

end


text:
abstract type Integration end

abstract type Model{I <: Integration, T <: Time} end

"""
	struct TemplateModel{I, T} <: Model{I, T}
		n::Int # state dimension
		m::Int # control dimension
		d::Int # disturbance dimension

		# additional model parameters ...
	end

	function f(model::TemplateModel, x, u, w)
		# continuous-time dynamics
		nothing
	end

	function k(model::TemplateModel, x)
		# kinematics
		nothing
	end

	model = TemplateModel{Integration, Time}(0, 0, 0)

	state_output(model, x) = x
	control_output(model, u) = u
"""

include_model(str::String) = include(joinpath(pwd(), "models", str * ".jl"))

"""
	propagate dynamics with implicit integrator
	- Levenberg-Marquardt
	- Newton
"""

function propagate_dynamics(model, x, u, w, h, t;
		solver = :levenberg_marquardt,
		tol_r = 1.0e-8, tol_d = 1.0e-6)

	res(z) = fd(model, z, x, u, w, h, t)

	return eval(solver)(res, copy(x), tol_r = tol_r, tol_d = tol_d)
end

function propagate_dynamics_jacobian(model, x, u, w, h, t;
		solver = :levenberg_marquardt,
		tol_r = 1.0e-8, tol_d = 1.0e-6)

	y = propagate_dynamics(model, x, u, w, h, t,
			solver = solver,
			tol_r = tol_r, tol_d = tol_d)

    dy(z) = fd(model, z, x, u, w, h, t)
	dx(z) = fd(model, y, z, u, w, h, t)
	du(z) = fd(model, y, x, z, w, h, t)

	Dy = ForwardDiff.jacobian(dy, y)
	A = -1.0 * Dy \ ForwardDiff.jacobian(dx, x)
	B = -1.0 * Dy \ ForwardDiff.jacobian(du, u)

	return y, A, B
end

"""
	get gravity compensating torques
"""
function initial_torque(model, q1, h;
	solver = :newton,
	tol_r = 1.0e-8, tol_d = 1.0e-6)

	x = [q1; q1]
    res(z) = fd(model, x, x, z, zeros(model.d), h, 0)

	return eval(solver)(res, 1.0e-5 * rand(model.m),
		tol_r = tol_r, tol_d = tol_d)
end

# Model conversions

"""
	free final time model
"""
function free_time_model(model::Model{I, FixedTime}) where I <: Integration
	model_ft = typeof(model).name.wrapper{I, FreeTime}([f == :m ? getfield(model,f) + 1 : getfield(model,f)
	 	for f in fieldnames(typeof(model))]...)
	return model_ft
end

function free_time_model(model::Model{I, FreeTime}) where I <: Integration
	@warn "Model is already free time"
	return model
end

"""
	no slip model
"""
function no_slip_model(model)
	# modify parameters
	m = model.nu + model.nc + model.nb + model.ns
	idx_ψ = (1:0)
	idx_η = (1:0)
	idx_s = model.nu + model.nc + model.nb .+ (1:model.ns)

	# assemble update parameters
	params = []
	for f in fieldnames(typeof(model))
		if f == :m
			push!(params, m)
		elseif f == :idx_ψ
			push!(params, idx_ψ)
		elseif f == :idx_η
			push!(params, idx_η)
		elseif f == :idx_s
			push!(params, idx_s)
		else
			push!(params, getfield(model,f))
		end
	end

	return typeof(model)(params...)
end

state_output(model, x) = x


text:
# BEGIN 0.7 deprecations

# PR #22062
function set_remote_url(repo::LibGit2.GitRepo, url::AbstractString; remote::AbstractString="origin")
    Base.depwarn(string(
        "`LibGit2.set_remote_url(repo, url; remote=remote)` is deprecated, use ",
        "`LibGit2.set_remote_url(repo, remote, url)` instead."), :set_remote_url)
    set_remote_url(repo, remote, url)
end

function set_remote_url(path::AbstractString, url::AbstractString; remote::AbstractString="origin")
    Base.depwarn(string(
        "`LibGit2.set_remote_url(path, url; remote=remote)` is deprecated, use ",
        "`LibGit2.set_remote_url(path, remote, url)` instead."), :set_remote_url)
    set_remote_url(path, remote, url)
end

function prompt(msg::AbstractString; default::AbstractString="", password::Bool=false)
    Base.depwarn(string(
        "`LibGit2.prompt(msg::AbstractString; default::AbstractString=\"\", password::Bool=false)` is deprecated, use ",
        "`result = Base.prompt(msg, default=default, password=password); result === nothing ? \"\" : result` instead."), :prompt)
    coalesce(Base.prompt(msg, default=default, password=password), "")
end

# PR #23640
# when this deprecation is deleted, remove all calls to it, and replace all keywords of:
# `payload::Union{CredentialPayload, AbstractCredential, CachedCredentials, Nothing}`
#  with `payload::CredentialPayload` from base/libgit2/libgit2.jl
function deprecate_nullable_creds(f, sig, payload)
    if isa(payload, Union{AbstractCredential, CachedCredentials, Nothing})
        # Note: Be careful not to show the contents of the credentials as it could reveal a
        # password.
        if payload === nothing
            msg = "`LibGit2.$f($sig; payload=nothing)` is deprecated, use "
            msg *= "`LibGit2.$f($sig; payload=LibGit2.CredentialPayload())` instead."
            p = CredentialPayload()
        else
            cred = payload
            C = typeof(cred)
            msg = "`LibGit2.$f($sig; payload=$C(...))` is deprecated, use "
            msg *= "`LibGit2.$f($sig; payload=LibGit2.CredentialPayload($C(...)))` instead."
            p = CredentialPayload(cred)
        end
        Base.depwarn(msg, f)
    else
        p = payload::CredentialPayload
    end
    return p
end

@deprecate get_creds!(cache::CachedCredentials, credid, default) get!(cache, credid, default)

@eval Base @deprecate merge!(repo::$(GitRepo), args...; kwargs...) $(LibGit2.merge!)(repo, args...; kwargs...)
@eval Base @deprecate push!(w::$(GitRevWalker), arg) $(LibGit2.push!)(w, arg)

# PR #24594
@deprecate AbstractCredentials AbstractCredential false
@deprecate UserPasswordCredentials UserPasswordCredential false
@deprecate SSHCredentials SSHCredential false

max_stars_repo_name: BenChung/Julia-BTI

max_issues_repo_name: BenChung/Julia-BTI

max_forks_repo_name: BenChung/Julia-BTI


text: struct div a end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
# source code representation cannot be recovered, but code_lowered yes
# see for example Symbolics.jl or Cassette.jl

# unary minus, plus

# some unary functions:
# sqrt, cbrt, exp, log, log2, log10

# some trig functions:
# sin    cos    tan    cot    sec    csc
# sinh   cosh   tanh   coth   sech   csch
# asin   acos   atan   acot   asec   acsc
# asinh  acosh  atanh  acoth  asech  acsch
# sinc   cosc

# Implement:
# unary minus, plus
# sqrt, cbrt, exp, log, log2, log10
# sin    cos    tan

# simplifying unary minus

module symbolics

function default_variable_value()
    return 0.0
end

abstract type AbstractVariable end
abstract type AbstractSymbolicExpression end
abstract type AbstractConstant end

include("is_operators.jl")
include("constant.jl")
include("variable.jl")
include("symbolic_expression.jl")
include("basic_operations.jl")
include("to_string.jl")
include("derivative.jl")

function is_equal(x, y)
    return false
end

function is_equal(x::Constant, y::Constant)
    return x.value == y.value
end

function is_equal(x::Variable, y::Variable)
    return x.name == y.name
end

function is_equal(x::SymbolicExpression, y::SymbolicExpression)
    # how to deal with permutations on expressions?
    # e.g.
    # a + b ?= b + a, should be, but
    # (a + b) + c ?= a + (b + c); this should be also true, but much harder to operate on
    # a + b + c + d -> possibly could be represented as +(a, b, c, d)
    # but then it is required to rethink the binary operations
    return false
end

# this is not what I want
# I dont want to compare variables by pointer -> already provided by Julia
# but I'd like to compare whether they are equaivalent
#
# function Base.:(==)(x::Constant, y::Constant)
#     # by value
#     return x.value == y.value
# end

# function Base.:(==)(x::Variable, y::Variable)
#     # by pointer
#     return x == y
# end

# function Base.:(==)(x::SymbolicExpression, y::SymbolicExpression)
#     # by pointer
#     return x == y
# end

function __simplify_plus(ex)
    # op +
    # ex + 0 -> ex
    # 0 + ex -> ex
    # c1 + c2 -> c3(c1.value + c2.value)
    # x + x -> 2 * x
    # ex + ex -> 2 * ex

    is_zero = is_zero_constant.(ex.args)
    if all(is_zero)
        return Constant(0.0)
    end

    is_not_zero = .!is_zero
    if any(is_not_zero)
        arg_vec = ex.args[is_not_zero]
        if length(arg_vec) == 1
            simpli_expr = arg_vec[1]
        else
            b_consts = is_constant.(arg_vec)
            idx      = findall(b_consts)
            if !isempty(idx)
                c        = reduce(+, arg_vec[idx])
                arg_vec  = [arg_vec[.!b_consts]; c]
            end
            simpli_expr = SymbolicExpression(:+, arg_vec...)
        end
        return simpli_expr
    end

    return ex
end

function __simplify_minus(ex)
    # op -
    # ex - 0 -> ex
    # 0 - ex -> -ex
    # c1 - c2 -> c3(c1.value - c2.value)
    # x - x -> 0
    # ex - ex -> 0

    L = length(ex.args)
    if L == 1
        if is_constant(ex.args[1])
            return Constant( -ex.args[1].value )
        else
            return ex
        end
    end

    if is_zero_constant(ex.args[1])
        return -ex.args[2]
    end

    if is_zero_constant(ex.args[2])
        return ex.args[1]
    end

    if is_equal(ex.args[1], ex.args[2])
        return Constant(0.0)
    end

    lhs = ex.args[1]
    rhs = ex.args[2]
    if is_constant(lhs) && is_constant(rhs)
        return lhs - rhs
    end

    return ex
end

function __simplify_product(ex)
    # op *
    # 0 * ex -> 0
    # ex * 0 -> 0
    # 1 * ex -> ex
    # ex * 1 -> ex

    is_zero = is_zero_constant.(ex.args)
    if any(is_zero)
        return Constant(0.0)
    end

    is_one     = is_one_constant.(ex.args)
    is_not_one = .!is_one
    if any(is_not_one)
        arg_vec = ex.args[is_not_one]
        if length(arg_vec) == 1
            simpli_expr = arg_vec[1]
        else

            b_consts = is_constant.(arg_vec)
            idx      = findall(b_consts)
            if !isempty(idx)
                c        = reduce(*, arg_vec[idx])
                arg_vec  = [arg_vec[.!b_consts]; c]
            end
            simpli_expr = SymbolicExpression(:*, arg_vec...)
        end
        return simpli_expr
    end

    return ex
end

function __simplify_divide(ex)
    # op /
    # ex / ex -> 1
    # 0 / rhs -> 0
    # rhs / 0 -> inf
    #
    # if (subexpr1 * subexpr2) / (subexpr1 * subexpr3)
    # -> subexpr2 / subexpr3

    if is_zero_constant(ex.args[1])
        return Constant(0.0)
    end

    if is_zero_constant(ex.args[2])
        return Constant(Inf)
    end

    if is_equal(ex.args[1], ex.args[2])
        return Constant(1.0)
    end

    lhs = ex.args[1]
    rhs = ex.args[2]
    if is_constant(lhs) && is_constant(rhs)
        return lhs / rhs
    end

    return ex
end

function __simplify_power(ex)
    # op ^
    # x^i / x^j -> x^(i-j)
    # x^0 -> 1
    # x^1 -> x
    # 1^x

    if is_zero_constant(ex.args[2])
        return simpli_expr = Constant(1.0)
    end

    if is_one_constant(ex.args[1])
        return simpli_expr = Constant(1.0)
    end

    if is_one_constant(ex.args[2])
        return simpli_expr = ex.args[1]
    end

    lhs = ex.args[1]
    rhs = ex.args[2]
    if is_constant(lhs) && is_constant(rhs)
        return lhs^rhs
    end

    return ex
end

function simplify(ex::Constant)
    return ex
end

function simplify(ex::Variable)
    return ex
end

function simplify(ex::SymbolicExpression)

    # rewriting rules
    # zero element
    # one element
    #
    # harder ones
    # commutative laws
    # associative laws
    #
    # a + b -> b + a
    # (a + b) + c -> a + (b + c)

    return __rec_simplify(ex)
end

function __rec_simplify(ex::Constant)
    return ex
end

function __rec_simplify(ex::Variable)
    return ex
end

function __rec_simplify(ex::SymbolicExpression)

    L           = length(ex.args)
    sym         = ex.head
    simpli_args = Vector{Any}(undef, L)

    for ii = 1:L
        simpli_args[ii] = __rec_simplify(ex.args[ii])
    end

    simpli_expr = SymbolicExpression(sym, simpli_args...)

    if is_plus_operator(sym)
        simpli_expr = __simplify_plus(simpli_expr)
    elseif is_minus_operator(sym)
        simpli_expr = __simplify_minus(simpli_expr)
    elseif is_product_operator(sym)
        simpli_expr = __simplify_product(simpli_expr)
    elseif is_divide_operator(sym)
        simpli_expr = __simplify_divide(simpli_expr)
    elseif is_power_operator(sym)
        simpli_expr = __simplify_power(simpli_expr)
    # else
    #     simpli_expr = ex
    end

    return simpli_expr
end

function simplify(ex_vec::Vector)
    L = length(ex_vec)
    simpli_vec = Vector{Any}(undef, L)
    for ii = 1:L
        simpli_vec[ii] = simplify(ex_vec[ii])
    end
    return simpli_vec
end

function simplify(ex_mat::Matrix)
    L = length(ex_mat)
    (n_row, n_col) = size(ex_mat)
    simpli_mat = Matrix{Any}(undef, n_row, n_col)
    for ii = 1:L
        simpli_mat[ii] = simplify(ex_mat[ii])
    end
    return simpli_mat
end

end # module

module stest
using ..symbolics

# symbolics.@variable x 10
# symbolics.@variable y 0.2
# symbolics.@variable z 3
# e1 = x + y
# e2 = x / y
# e3 = z ^ (x + y)
# e4 = -x
# e5 = -x - y
# e6 = sin(x + y - z * x)
# e3_val = symbolics.evaluate(e3)
# str1 = symbolics.to_string( stest.e1 )
# str2 = symbolics.to_string( stest.e2 )
# str3 = symbolics.to_string( stest.e3 )
# str4 = symbolics.to_string( stest.e4 )
# str5 = symbolics.to_string( stest.e5 )
# str6 = symbolics.to_string( stest.e6 )
# e7 = y * 10
# str7 = symbolics.to_string( stest.e7 )
# e8 = 10 * y
# str8 = symbolics.to_string( stest.e8 )
# e9 = 10 + (x - y)
# str9 = symbolics.to_string( stest.e9 )
# # !!!
# # this doesnt work correctly now
# # e10 = y * 10 * (x * 2 - 3*z)
# e10 = x * y * x * (x * 2 - 3*z)
# str10 = symbolics.to_string( stest.e10 )

function paren_test()

    symbolics.@variable x 10
    symbolics.@variable y 0.2
    symbolics.@variable z 3

    e1 = x + y
    e2 = x / y

    e3 = z ^ (x + y)

    e4 = -x
    e5 = -x - y
    e6 = sin(x + y - z * x)

    # e3_val = symbolics.evaluate(e3)

    # str1 = symbolics.to_string( stest.e1 )
    # str2 = symbolics.to_string( stest.e2 )
    # str3 = symbolics.to_string( stest.e3 )
    # str4 = symbolics.to_string( stest.e4 )
    # str5 = symbolics.to_string( stest.e5 )
    # str6 = symbolics.to_string( stest.e6 )

    # e7 = y * 10
    # str7 = symbolics.to_string( stest.e7 )
    # e8 = 10 * y
    # str8 = symbolics.to_string( stest.e8 )

    # e9 = 10 + (x - y)
    # str9 = symbolics.to_string( stest.e9 )

    # !!!
    # this doesnt work correctly now
    # e7 = y * 10 * (x * 2 - 3*z)
    e7 = x * y * x * (x * 2 - 3*z)
    str7 = symbolics.to_string( stest.e7 )

    # str = symbolics.to_string( stest.e6 )

    # x_val = symbolics.evaluate(x)
    # y_val = symbolics.evaluate(x)

    # e1_val = symbolics.evaluate(e1)

end

function test_derivative_simpli()
    symbolics.@variable x 10
    symbolics.@variable y 0.2
    symbolics.@variable z 3
    symbolics.@variable w 5

    u = -x

    ex1 = x + y
    ex2 = x * y + z
    ex3 = x + y + z
    ex4 = x * y * z
    ex5 = sin(x)

    # str = symbolics.to_string(ex2)

    der_ex3       = symbolics.derivative(stest.ex3, stest.x)
    der_str3      = symbolics.to_string(der_ex3)
    @time s_der_ex3     = symbolics.simplify(der_ex3)
    str_s_der_ex3 = symbolics.to_string(s_der_ex3)

    der_ex4       = symbolics.derivative(stest.ex4, stest.x)
    der_str4      = symbolics.to_string(der_ex4)
    @time s_der_ex4     = symbolics.simplify(der_ex4)
    str_s_der_ex4 = symbolics.to_string(s_der_ex4)
end

function harmonic_osc(der_q, time, q)
    der_q[1] = q[2]
    der_q[2] = -5.0 * q[1] - q[2]

    return der_q
end

# symbolics.@variable x 10
# symbolics.@variable v 0.2

# q = Vector{symbolics.Variable}(undef, 2)
# der_q = Vector{Any}(undef, 2)

# q[1] = x
# q[2] = v

# der_q = harmonic_osc(der_q, 0.0, q)

# jac = symbolics.jacobian(der_q, q)
# sjac = symbolics.simplify(jac)
# sjac = symbolics.simplify(sjac)

# symbolics.@variable a 30.0

# expr1 = a + 1.0

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
module ErrorsWrapper

using JavaCall

const Class = @jimport java.lang.Class
const Enum = @jimport java.lang.Enum
const JString = @jimport java.lang.String
const Locale = @jimport java.util.Locale
const Localizable = @jimport org.hipparchus.exception.Localizable
const LocalizedException = @jimport org.orekit.errors.LocalizedException
const Object = @jimport java.lang.Object
const Optional = @jimport java.util.Optional
const OrekitMessages = @jimport org.orekit.errors.OrekitMessages

export LocalizedException
export OrekitMessages
export compare_to
export describe_constable
export equals
export get_class
export get_declaring_class
export get_localized_string
export get_message
export get_parts
export get_source_string
export get_specifier
export hash_code
export name
export notify
export notify_all
export ordinal
export to_string
export value_of
export values
export wait

include("localized_exception.jl")
include("orekit_messages.jl")

end


text:
type NullCrossover <: Crossover
  stage::AbstractString
  rep::Representation

  NullCrossover(stage::AbstractString, rep::Representation) =
    new(stage, rep)
end

"""
Used to provide a definition for a null crossover operation.
"""
type NullCrossoverDefinition <: CrossoverDefinition
  stage::AbstractString
end

"""
Composes a null crossover operator from a provided definition, for a given
species.
"""
function compose!(c::NullCrossoverDefinition, sp::Species)
  c.stage = c.stage == "" ? genotype(sp).label : c.stage
  r = sp.stages[c.stage].representation
  NullCrossover(c.stage, r)
end

"""
Null crossover operators accept a pair of chromosomes and return the same pair
of chromosomes untouched. Useful for debugging purposes.

**Positional Arguments:**

* `stage::AbstractString`, the name of the stage that this operator should use.
  Defaults to the genotype if no stage is specified.
"""
null() = null("")
null(stage::AbstractString) = NullCrossoverDefinition(stage)

"""
Performs a null crossover operation on two provided chromosomes, returning the
same chromosomes unchanged.
"""
operate!{T}(::NullCrossover, outputs::Vector{IndividualStage{T}}, x::T, y::T) =
  outputs


text:
using LegolasFlux
using Test
using Flux, LegolasFlux
using LegolasFlux: Weights, FlatArray, ModelRow
using Arrow
using Random
using StableRNGs

function make_my_model()
    return Chain(Dense(1, 10), Dense(10, 10), Dense(10, 1))
end

function test_weights()
    shapes = [(10, 1), (10,), (10, 10), (10,), (1, 10), (1,)]
    return [reshape(Float32.(1:prod(s)), s) for s in shapes]
end

# This simple model should work with both Flux's `params/loadparams!` and
# our `weights/load_weights!`. The only difference is in layers with `!isempty(other_weights(layer))`.
@testset "using ($get_weights, $load_weights)" for (get_weights, load_weights) in [(fetch_weights, load_weights!, params, Flux.loadparams!)]
    my_model = make_my_model()
    load_weights(my_model, test_weights())

    model_row = ModelRow(; weights=collect(get_weights(my_model)))
    write_model_row("my_model.model.arrow", model_row)

    fresh_model = make_my_model()

    model_row = read_model_row("my_model.model.arrow")
    weights = collect(model_row.weights)
    load_weights(fresh_model, weights)

    @test collect(params(fresh_model)) == weights == test_weights()

    @test all(x -> eltype(x) == Float32, weights)

    rm("my_model.model.arrow")
end

struct MyArrayModel
    dense_array::Array
end
Flux.@functor MyArrayModel

@testset "Non-numeric arrays ignored" begin
    m = MyArrayModel([Dense(1, 10), Dense(10, 10), Dense(10, 1)])
    weights = fetch_weights(m)
    @test length(weights) == 6

    model_row = ModelRow(; weights=collect(weights))
    write_model_row("my_model.model.arrow", model_row)

    new_model_row = read_model_row("my_model.model.arrow")
    new_weights = collect(new_model_row.weights)
    @test new_weights == weights
end

@testset "Errors" begin
    my_model = make_my_model()
    w = test_weights()
    w[end] = []
    @test_throws ArgumentError load_weights!(my_model, w)

    w = test_weights()
    push!(w, [])
    @test_throws ArgumentError load_weights!(my_model, w)
end

@testset "`Weights`" begin
    rng = StableRNG(245)
    v = [rand(rng, Int8, 5), rand(rng, Float32, 5, 5)]
    @test Weights(v) isa Weights{Float32}
    @test Weights(FlatArray{Float32}.(v)) isa Weights{Float32}
    @test Weights(FlatArray{Float64}.(v)) isa Weights{Float64}

    w = Weights(v)
    tbl = [(; weights = w)]
    @test Arrow.Table(Arrow.tobuffer(tbl)).weights[1] == w
end

@testset "`flux_workarounds`" begin
    @testset "layer $layer" for layer in [BatchNorm, InstanceNorm, (c) -> GroupNorm(c, 1), c -> identity]
        mk_model = () -> (Random.seed!(1); Chain(Dense(1, 10), Dense(10, 10), layer(1), Dense(10, 1)))
        model = mk_model()
        trainmode!(model)
        x = reshape([1f0], 1, 1, 1)
        for i in 1:10
            x = model(x)
        end
        testmode!(model)
        w = fetch_weights(model)
        p = collect(params(model))
        output = model(x)

        r1 = mk_model()
        load_weights!(r1, w)
        testmode!(r1)

        @test output ≈ r1(x)

        if layer == BatchNorm
            r2 = mk_model()
            Flux.loadparams!(r2, p)
            testmode!(r2)

            # If this test *fails*, meaning `output ≈ r2(x)`,
            # then perhaps we should revisit `load_weights!`
            # and could consider switching to `Flux.loadparams`.
            # See https://github.com/beacon-biosignals/LegolasFlux.jl/pull/4
            # for more.
            @test_broken output ≈ r2(x)
        end
    end
end

@testset "Example" begin
    include("../examples/digits.jl")
end


text:
"""
    struct GOCLoad <: AbstractNodeLabel

Structure descendant of AbstractNodeLabel

# Fields
- `busname::String`
- `id::String`: necessary if several loads at a bus
- `load::Complex128`: Sd = Pd + im Qd

"""
struct GOCLoad <: AbstractNodeLabel
    busid::Int64
    id::String #label
    load::Complex #Sd = Pd+im*Qd
end

"""
    create_vars!(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}, scenario::String) where T <: GOCLoad

Create load variable at `bus` for `elemid` in `bus_vars` if `elem_formulation == : NewVar`. \n
Return nothing
"""
function create_vars!(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}, scenario::String) where T <: GOCLoad
    if elem_formulation == :NewVar
      bus_vars[elemid] = Variable(variable_name("Sload", bus, elemid, scenario), Complex)
    end
    return
end

"""
    [cstrname, polynom, lb, ub] = Snodal(element::T, busid::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}) where T <: GOCLoad

Return the polynomial contribution in power of load `elemid` at `busid`. Will be used to construct power balance constraints in polynomial problem.\n
If `elem_formulation == :NbMinVar`, return value of Sload `[cstrname, Sload, 0, 0]`\n
If `elem_formulation == :NewVar`, return variable Sload `[cstrname, Polynomial(Sload), 0, 0]`\n
Return `["", Polynomial(), 0, 0]` otherwise.
"""
function Snodal(element::T, busid::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}) where T <: GOCLoad
    cstrname = "LOAD"
    if elem_formulation == :NbMinVar
      Sload = element.load
      return [cstrname, Sload, 0, 0]
    elseif elem_formulation == :NewVar
      return [cstrname, Polynomial(bus_vars[elemid]), 0, 0]
    else
        return ["", Polynomial(), 0, 0]
    end
end

"""
    constraint(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}, scenario::String, OPFpbs::OPFProblems) where T <: GOCLoad

Return a constraint to define load variable for `elemid` at `bus` if `elem_formulation == : NewVar`: `"LOAD"=> Sload==value` \n
Return empty dictionary of constraints otherwise.

"""
function constraint(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::SortedDict{String, Variable}, scenario::String, OPFpbs::OPFProblems) where T <: GOCLoad
    if elem_formulation == :NewVar
      cstrname = get_LoadDef_cstrname()
      Sload = element.load
      return SortedDict{String, Constraint}(cstrname => Polynomial(bus_vars[elemid]) == Sload)
    end
    return SortedDict{String, Constraint}()
end


text:
export readlsm

const _default_delims = [' ','\t']

"""
    readlsm([type,] io, delim=[' ', '\t'], comment='#')
    readlsm([type,] filepath, delim=[' ', '\t'], comment='#')

Reads a Labelled Square Matrix from an IO stream or file, returning a pair of
`(matrix, labels)`.

One may specify the element type of the resulting matrix; this defaults to
Float64.
"""
function readlsm{T}(::Type{T}, io::IO, delim=_default_delims, comment='#')
    matrix = Matrix{T}()
    labels = String[]
    nitems = 0
    row = 0
    for line in eachline(io)
        line = rstrip(line)
        if startswith(line, comment) || isempty(line)
            continue  # skip comments and empty lines
        end
        if row > nitems
            # Non-square error handled outside loop below
            break
        end
        cells = split(line, delim, keep=false)
        if isempty(labels)
            # header
            append!(labels, cells)
            nitems = length(labels)
            matrix = zeros(T, (nitems, nitems))
            continue
        end
        row += 1
        if cells[1] != labels[row]
            error("Cell label mismatch. (at row $row: $(cells[1]) != ",
                  "$(labels[row]))")
        end
        matrix[row, :] = map((s) -> parse(T, s), cells[2:end])
    end
    if row != nitems
        error("Non-square matrix: read $row rows, with $nitems columns.")
    end
    return (matrix, labels)
end

function readlsm(io::IO, delim=_default_delims)
    return readlsm(Float64, io, delim)
end

function readlsm{T}(::Type{T}, path::AbstractString, delim=_default_delims)
    open(path) do file
        return readlsm(T, file, delim)
    end
end

function readlsm(path::AbstractString, delim=_default_delims)
    return readlsm(Float64, path, delim)
end


text:
# Consider moving this into MLLabelUtils.labelenc

_ambiguous() = throw(ArgumentError("Can't infer the label encoding because argument types/values are ambigous. Please specify the desired LabelEnc manually."))

"""
    _labelenc(targets, outputs)

Tries to automatically infer the compare mode for the given
structure and values in `targets` and `outputs`.
"""
_labelenc(target, output) = _ambiguous()

_labelenc(targets::AbstractArray{Bool}, outputs::AbstractArray{Bool}) =
    LabelEnc.FuzzyBinary()

# Generate the Bool combinations to avoid ambiguity warnings
_labelenc(target::Bool, output::Bool) = LabelEnc.FuzzyBinary()
for _T2 in (:Bool, :Real, :Any)
    @eval begin
        _labelenc(targets::AbstractArray{T2}, outputs::AbstractArray{Bool}) where {T2<:$_T2} = LabelEnc.FuzzyBinary()
        _labelenc(targets::AbstractArray{Bool}, outputs::AbstractArray{T2}) where {T2<:$_T2} = LabelEnc.FuzzyBinary()
    end
    if _T2 != :Bool
        @eval begin
            _labelenc(target::$_T2, output::Bool) =
                LabelEnc.FuzzyBinary()
            _labelenc(target::Bool, output::$_T2) =
                LabelEnc.FuzzyBinary()
        end
    end
end

# Not enough information available to decide for fuzzy binary,
# because we don't want to infer an arbitrary label as positive.
# In this case we decide to let labelenc decide.
function _labelenc(targets::AbstractArray, outputs::AbstractArray)
    labelenc(vcat(vec(targets), vec(outputs)))
end


text:
function getclass(obj::JavaObject)
    jcall(obj, "getClass", JClass, ())
end

function conventional_name(name::AbstractString)
    if startswith(name, "[")
        return conventional_name(name[2:end]) * "[]"
    elseif name == "Z"
        return "boolean"
    elseif name == "B"
        return "byte"
    elseif name == "S"
        return "short"
    elseif name == "C"
        return "char"
    elseif name == "I"
        return "int"
    elseif name == "J"
        return "long"
    elseif name == "F"
        return "float"
    elseif name == "D"
        return "double"
    elseif name == "V"
        return "void"
    elseif startswith(name, "L")
        return name[2:end-1]
    else
        return name
    end
end

"""
```
getname(cls::JClass)
```
Returns the fully qualified name of the java class

### Args
* cls: The JClass object

### Returns
The fully qualified name of the java class
"""
function getname(cls::JClass)
    rawname = jcall(cls, "getName", JString, ())
    conventional_name(rawname)
end

"""
```
getname(method::JMethod)
```
Returns the fully qualified name of the java method

### Args
* cls: The JClass method

### Returns
The fully qualified name of the method
"""
function getname(method::JMethod)
    jcall(method, "getName", JString, ())
end

"""
"""
getname(field::JField) = jcall(field, "getName", JString, ())

"""
```
listmethods(obj::JavaObject)
```
Lists the methods that are available on the java object passed

### Args
* obj: The java object

### Returns
List of methods
"""
function listmethods(obj::JavaObject)
    cls = getclass(obj)
    jcall(cls, "getMethods", Vector{JMethod}, ())
end

function listmethods(::Type{JavaObject{C}}) where C
    cls = classforname(string(C))
    jcall(cls, "getMethods", Vector{JMethod}, ())
end

function listmethods(cls::JClass)
    jcall(cls, "getMethods", Vector{JMethod}, ())
end

"""
```
listmethods(obj::JavaObject, name::AbstractString)
```
Lists the methods that are available on the java object passed. The methods are filtered based on the name passed

### Args
* obj: The java object
* name: The filter (e.g. method name)

### Returns
List of methods available on the java object and matching the name passed
"""
function listmethods(obj::Union{JavaObject{C}, Type{JavaObject{C}}}, name::AbstractString) where C
    allmethods = listmethods(obj)
    filter(m -> getname(m) == name, allmethods)
end

"""
```
getreturntype(method::JMethod)
```
Returns the return type of the java method

### Args
* method: The java method object

### Returns
Returns the type of the return object as a JClass
"""
function getreturntype(method::JMethod)
    jcall(method, "getReturnType", JClass, ())
end

"""
gettype(field::JField)

Get type of field
"""
gettype(field::JField) = jcall(field, "getType", JClass, ())

"""
```
listfields(obj::JavaObject)
```
List the fields that are available on the java object passed.
"""
listfields(cls::AbstractString) = listfields(classforname(cls))
listfields(cls::Type{JavaObject{C}}) where C = listfields(classforname(string(C)))
listfields(cls::JClass) = jcall(cls, "getFields", Vector{JField}, ())
listfields(obj::JavaObject) = listfields(getclass(obj))

function listfields(cls::Union{JavaObject{C}, Type{JavaObject{C}}}, name::AbstractString) where C
    allfields = listfields(cls)
    filter(f -> getname(f) == name, allfields)
end

"""
```
getparametertypes(method::JMethod)
```
Returns the parameter types of the java method

### Args
* method: The java method object

### Returns
Vector the parametertypes
"""
function getparametertypes(method::Union{JMethod,JConstructor})
    jcall(method, "getParameterTypes", Vector{JClass}, ())
end

function Base.show(io::IO, method::JMethod)
    name = getname(method)
    rettype = getname(getreturntype(method))
    argtypes = [getname(c) for c in getparametertypes(method)]
    argtypestr = join(argtypes, ", ")
    print(io, "$rettype $name($argtypestr)")
end

function Base.show(io::IO, field::JField)
    name = getname(field)
    fieldtype = getname(gettype(field))
    print(io, "$fieldtype $name")
end

"""
```
classforname(name::String)
```
Create an instance of `Class<name>` (same as `Class.forName(name)` in Java)

### Args
* name: The name of a class to instantiate

### Returns
JavaObject Instance of `Class<name>`
"""
function classforname(name::String)
    thread = jcall(JThread, "currentThread", JThread, ())
    loader = jcall(thread, "getContextClassLoader", JClassLoader, ())
    return jcall(JClass, "forName", JClass, (JString, jboolean, JClassLoader),
                 name, true, loader)
end


text:
#== import modules here ==#
using PyCall
psi4 = pyimport("psi4")

using HDF5

#== globals ==#
const THRESHOLD = 0.00001

#== psi4 setup ==#
psi4.set_memory("500 MB")
psi4.core.set_output_file("psi4_compare.dat")

geom = psi4.geometry("
  O      1.00339        0.03093        0.09503
  O      2.32339        0.03093        0.09503
  H      0.68006        -0.87208        -0.04962
  H      2.64672        0.93395        0.23969
  symmetry=c1
")

options = Dict(
  "REFERENCE" => "RHF",
  "BASIS" => "PCSeg-1",
  "SCF_TYPE"  => "DIRECT",
  "GUESS" => "SAD",
  "SAD_SCF_TYPE" => "DIRECT",
  "DF_SCF_GUESS" => false,
  "PUREAM" => false,
  "E_CONVERGENCE" => 1e-10,
  "D_CONVERGENCE" => 1e-10,
  "PRINT" => 4
)
psi4.set_options(options)

name = "scf/"*options["BASIS"]
scf_e, scf_wfn = psi4.energy(name, molecule=geom, return_wfn=true)

S_psi4 = scf_wfn.S().to_array()

#== standalone setup ==#
debug = h5open("debug.h5", "r")

tmp = debug["RHF"]["Iteration-None"]["S"][:]
S_standalone = reshape(tmp,(Int64(sqrt(length(tmp))),
  Int64(sqrt(length(tmp)))))
#display(S_standalone)

close(debug)

#== comparison ==#
S_diff = S_standalone .- S_psi4
diff_array = []

println("Comparing Psi4 overlap to standalone overlap...")
for i in 1:size(S_diff)[1], j in 1:size(S_diff)[2]
  index = size(S_diff)[1]*(j-1) + i
  @assert isapprox(S_diff[i,j], S_diff[index])

  if S_diff[i,j] > THRESHOLD
    println((i,j), " => ", index, ": ", S_standalone[i,j], ", ",
      S_psi4[i,j], ", ", S_diff[i,j])
    push!(diff_array,(i,j))
  end
end
println()

#== check for symmetry issues ==#
println("Checking for symmetry issues in S_diff...")
diff_symm_array = []
for tuple in diff_array
  if S_diff[tuple...] - S_diff[reverse(tuple)...] > THRESHOLD
    push!(diff_symm_array,tuple)
    println(tuple, ": ", S_diff[tuple...], ", ",
      S_diff[reverse(tuple)...])
  end
end
println()

#== check diagonal elements ==#
println("Printing values of diagonal elements...")

#leng = Int64(sqrt(length(S_diff)))
#for i in 0:(leng-1)
#  index = (leng+1)*i + 1
#  println(index, " ", S_standalone[index], " ", S_psi4[index], " ", S_diff[index])
#end

for i in 1:size(S_diff)[1]
  index = size(S_diff)[1]*(i-1) + i
  println((i,i), " => ", index,  ": ", S_standalone[i, i], ", ",
    S_psi4[i, i], ", ", S_diff[i, i])
end
println()


text:
module CompatIterators
# Defining this module with the name `CompatIterators` and define
# `const Iterators = CompatIterators` at the end of this file so that
#
#     julia> using Compat.Iterators
#
#     julia> Iterators
#     Base.Iterators
#
#     julia> takewhile
#     takewhile (generic function with 1 method)
#
# works without conflict in identifiers in all Julia versions.  Users
# still need to explicitly import `Compat.Iterators` via `using
# Compat: Iterators` to use, e.g., `Iterators.map` before Julia 1.6.
# This is the case with and without the `CompatIterators` hack.

using Base: SizeUnknown

# We normally use `Base.f(...) = ...` style for overloading.  However,
# we use `import Base: f` style here to make synchronizing the code
# with `Base` easier.
import Base: IteratorEltype, IteratorSize, eltype, iterate

# Import exported APIs
for n in names(Base.Iterators)
    n === :Iterators && continue
    @eval begin
        using Base.Iterators: $n
        export $n
    end
end

# Import unexported public APIs
using Base.Iterators: filter

# https://github.com/JuliaLang/julia/pull/33437
if VERSION < v"1.4.0-DEV.291"  # 5f013d82f92026f7dfbe4234f283658beb1f8a2a
    export takewhile, dropwhile

    # takewhile
    struct TakeWhile{I,P<:Function}
        pred::P
        xs::I
    end

    takewhile(pred,xs) = TakeWhile(pred,xs)

    function iterate(ibl::TakeWhile, itr...)
        y = iterate(ibl.xs,itr...)
        y === nothing && return nothing
        ibl.pred(y[1]) || return nothing
        y
    end

    IteratorSize(::Type{<:TakeWhile}) = SizeUnknown()
    eltype(::Type{TakeWhile{I,P}}) where {I,P} = eltype(I)
    IteratorEltype(::Type{TakeWhile{I,P}}) where {I,P} = IteratorEltype(I)

    # dropwhile
    struct DropWhile{I,P<:Function}
        pred::P
        xs::I
    end

    dropwhile(pred,itr) = DropWhile(pred,itr)

    iterate(ibl::DropWhile,itr) = iterate(ibl.xs, itr)
    function iterate(ibl::DropWhile)
        y = iterate(ibl.xs)
        while y !== nothing
            ibl.pred(y[1]) || break
            y = iterate(ibl.xs,y[2])
        end
        y
    end

    IteratorSize(::Type{<:DropWhile}) = SizeUnknown()
    eltype(::Type{DropWhile{I,P}}) where {I,P} = eltype(I)
    IteratorEltype(::Type{DropWhile{I,P}}) where {I,P} = IteratorEltype(I)
end

# https://github.com/JuliaLang/julia/pull/34352
if VERSION < v"1.6.0-DEV.258"  # 1f8b44204fafb1caabc6a1cd6ca39458a550e2fc
    map(f, args...) = Base.Generator(f, args...)
else
    using Base.Iterators: map
end

end  # module

const Iterators = CompatIterators