text:
#using Pkg
# pkg"add Oceananigans CairoMakie JLD2"
ENV["GKSwstype"] = "100"
# pushfirst!(LOAD_PATH, @__DIR__)
# pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", "..")) # add Oceananigans

using Printf
using Statistics
using JLD2

using Oceananigans
using Oceananigans.Units
using Oceananigans.OutputReaders: FieldTimeSeries
using Oceananigans.Grids: xnode, ynode, znode
using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: CATKEVerticalDiffusivity

filename = "eddying_channel"

# Architecture
architecture = GPU()

# Domain
const Lx = 1000kilometers # zonal domain length [m]
const Ly = 2000kilometers # meridional domain length [m]
const Lz = 2kilometers    # depth [m]

# number of grid points
Nx = 64
Ny = 128
Nz = 40

save_fields_interval = 7days
stop_time = 5years
Δt₀ = 5minutes

#=
# # Stretched grid
# we implement here a linearly streched grid in which the top grid cell has Δzₜₒₚ
# and every other cell is bigger by a factor σ, e.g.,
# Δzₜₒₚ, Δzₜₒₚ * σ, Δzₜₒₚ * σ², ..., Δzₜₒₚ * σᴺᶻ⁻¹,
# so that the sum of all cell heights is Lz

# Given Lz and stretching factor σ > 1 the top cell height is Δzₜₒₚ = Lz * (σ - 1) / σ^(Nz - 1)

σ = 1.04 # linear stretching factor
Δz_center_linear(k) = Lz * (σ - 1) * σ^(Nz - k) / (σ^Nz - 1) # k=1 is the bottom-most cell, k=Nz is the top cell
linearly_spaced_faces(k) = k==1 ? -Lz : - Lz + sum(Δz_center_linear.(1:k-1))

grid = RectilinearGrid(architecture
                       topology = (Periodic, Bounded, Bounded),
                       size = (Nx, Ny, Nz),
                       halo = (3, 3, 3),
                       x = (0, Lx),
                       y = (0, Ly),
                       z = linearly_spaced_faces)

# The vertical spacing versus depth for the prescribed grid

plot(grid.Δzᵃᵃᶜ[1:Nz], grid.zᵃᵃᶜ[1:Nz],
     axis=(xlabel = "Vertical spacing (m)",
           ylabel = "Depth (m)"))

=#

# We choose a regular grid though because of numerical issues that yet need to be resolved
grid = RectilinearGrid(architecture = architecture,
                       topology = (Periodic, Bounded, Bounded),
                       size = (Nx, Ny, Nz),
                       halo = (3, 3, 3),
                       x = (0, Lx),
                       y = (0, Ly),
                       z = (-Lz, 0))

@info "Built a grid: $grid."

#####
##### Boundary conditions
#####

α  = 2e-4     # [K⁻¹] thermal expansion coefficient
g  = 9.8061   # [m s⁻²] gravitational constant
cᵖ = 3994.0   # [J K⁻¹] heat capacity
ρ  = 1024.0   # [kg m⁻³] reference density

parameters = (Ly = Ly,
              Lz = Lz,
              Qᵇ = 10 / (ρ * cᵖ) * α * g,          # buoyancy flux magnitude [m² s⁻³]
              y_shutoff = 5/6 * Ly,                # shutoff location for buoyancy flux [m]
              τ = 0.2/ρ,                           # surface kinematic wind stress [m² s⁻²]
              μ = 1 / 30days,                      # bottom drag damping time-scale [s⁻¹]
              ΔB = 8 * α * g,                      # surface vertical buoyancy gradient [s⁻²]
              H = Lz,                              # domain depth [m]
              h = 1000.0,                          # exponential decay scale of stable stratification [m]
              y_sponge = 19/20 * Ly,               # southern boundary of sponge layer [m]
              λt = 7.0days                         # relaxation time scale [s]
)

@inline function buoyancy_flux(i, j, grid, clock, model_fields, p)
    y = ynode(Center(), j, grid)
    return ifelse(y < p.y_shutoff, p.Qᵇ * cos(3π * y / p.Ly), 0.0)
end

buoyancy_flux_bc = FluxBoundaryCondition(buoyancy_flux, discrete_form=true, parameters=parameters)

@inline function u_stress(i, j, grid, clock, model_fields, p)
    y = ynode(Center(), j, grid)
    return - p.τ * sin(π * y / p.Ly)
end

u_stress_bc = FluxBoundaryCondition(u_stress, discrete_form=true, parameters=parameters)

@inline u_drag(i, j, grid, clock, model_fields, p) = @inbounds - p.μ * p.Lz * model_fields.u[i, j, 1]
@inline v_drag(i, j, grid, clock, model_fields, p) = @inbounds - p.μ * p.Lz * model_fields.v[i, j, 1]

u_drag_bc = FluxBoundaryCondition(u_drag, discrete_form=true, parameters=parameters)
v_drag_bc = FluxBoundaryCondition(v_drag, discrete_form=true, parameters=parameters)

b_bcs = FieldBoundaryConditions(top = buoyancy_flux_bc)

u_bcs = FieldBoundaryConditions(top = u_stress_bc, bottom = u_drag_bc)
v_bcs = FieldBoundaryConditions(bottom = v_drag_bc)

#####
##### Coriolis
#####

const f = -1e-4     # [s⁻¹]
const β =  1e-11    # [m⁻¹ s⁻¹]
coriolis = BetaPlane(f₀ = f, β = β)

#####
##### Forcing and initial condition
#####

@inline initial_buoyancy(z, p) = p.ΔB * (exp(z / p.h) - exp(-p.Lz / p.h)) / (1 - exp(-p.Lz / p.h))
@inline mask(y, p) = max(0.0, y - p.y_sponge) / (Ly - p.y_sponge)

@inline function buoyancy_relaxation(i, j, k, grid, clock, model_fields, p)
    timescale = p.λt
    y = ynode(Center(), j, grid)
    z = znode(Center(), k, grid)
    target_b = initial_buoyancy(z, p)
    b = @inbounds model_fields.b[i, j, k]
    return - 1 / timescale  * mask(y, p) * (b - target_b)
end

Fb = Forcing(buoyancy_relaxation, discrete_form = true, parameters = parameters)

# Turbulence closures

κh = 0.5e-5 # [m²/s] horizontal diffusivity
νh = 30.0   # [m²/s] horizontal viscocity
κz = 0.5e-5 # [m²/s] vertical diffusivity
νz = 3e-4   # [m²/s] vertical viscocity

horizontal_diffusivity = AnisotropicDiffusivity(νh=νh, νz=νz, κh=κh, κz=κz)

convective_adjustment = ConvectiveAdjustmentVerticalDiffusivity(convective_κz = 1.0,
                                                                convective_νz = 0.0)

catke = CATKEVerticalDiffusivity()

#####
##### Model building
#####

@info "Building a model..."

model = HydrostaticFreeSurfaceModel(grid = grid,
                                    free_surface = ImplicitFreeSurface(),
                                    momentum_advection = WENO5(),
                                    tracer_advection = WENO5(),
                                    buoyancy = BuoyancyTracer(),
                                    coriolis = coriolis,
                                    closure = (horizontal_diffusivity, catke),
                                    tracers = (:b, :e, :c),
                                    boundary_conditions = (b=b_bcs, u=u_bcs, v=v_bcs),
                                    forcing = (; b=Fb,)
                                    )

@info "Built $model."

#####
##### Initial conditions
#####

# resting initial condition
ε(σ) = σ * randn()
bᵢ(x, y, z) = parameters.ΔB * ( exp(z/parameters.h) - exp(-Lz / parameters.h) ) / (1 - exp(-Lz / parameters.h)) + ε(1e-8)
uᵢ(x, y, z) = ε(1e-8)
vᵢ(x, y, z) = ε(1e-8)
wᵢ(x, y, z) = ε(1e-8)

Δy = 100kilometers
Δz = 100
Δc = 2Δy
cᵢ(x, y, z) = exp(-y^2 / 2Δc^2) * exp(-(z + Lz/4)^2 / 2Δz^2)

set!(model, b=bᵢ, u=uᵢ, v=vᵢ, w=wᵢ, c=cᵢ)

#####
##### Simulation building

wizard = TimeStepWizard(cfl=0.1, max_change=1.1, max_Δt=20minutes)

wall_clock = [time_ns()]

function print_progress(sim)
    @printf("[%05.2f%%] i: %d, t: %s, wall time: %s, max(u): (%6.3e, %6.3e, %6.3e) m/s, next Δt: %s\n",
            100 * (sim.model.clock.time / sim.stop_time),
            sim.model.clock.iteration,
            prettytime(sim.model.clock.time),
            prettytime(1e-9 * (time_ns() - wall_clock[1])),
            maximum(abs, sim.model.velocities.u),
            maximum(abs, sim.model.velocities.v),
            maximum(abs, sim.model.velocities.w),
            prettytime(sim.Δt))

    wall_clock[1] = time_ns()

    return nothing
end

simulation = Simulation(model, Δt=Δt₀, stop_time=stop_time)

# add timestep wizardrogress callback
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(20))

# add progress callback
simulation.callbacks[:print_progress] = Callback(print_progress, IterationInterval(20))

#####
##### Diagnostics
#####

u, v, w = model.velocities
b, c = model.tracers.b, model.tracers.c

ζ = ComputedField(∂x(v) - ∂y(u))

B = AveragedField(b, dims=1)
U = AveragedField(u, dims=1)
V = AveragedField(v, dims=1)
W = AveragedField(w, dims=1)

b′ = b - B
v′ = v - V
w′ = w - W

v′b′ = AveragedField(v′ * b′, dims=1)
w′b′ = AveragedField(w′ * b′, dims=1)

outputs = (; b, ζ, u)

averaged_outputs = (; v′b′, w′b′, B, U)

#####
##### Build checkpointer and output writer
#####

simulation.output_writers[:checkpointer] = Checkpointer(model,
                                                        schedule = TimeInterval(1years),
                                                        prefix = filename,
                                                        force = true)

slicers = (west = FieldSlicer(i=1),
           east = FieldSlicer(i=grid.Nx),
           south = FieldSlicer(j=1),
           north = FieldSlicer(j=grid.Ny),
           bottom = FieldSlicer(k=1),
           top = FieldSlicer(k=grid.Nz))

for side in keys(slicers)
    field_slicer = slicers[side]

    simulation.output_writers[side] = JLD2OutputWriter(model, outputs,
                                                       schedule = TimeInterval(save_fields_interval),
                                                       field_slicer = field_slicer,
                                                       prefix = filename * "_$(side)_slice",
                                                       force = true)
end

simulation.output_writers[:zonal] = JLD2OutputWriter(model, (b=B, u=U),#, v=V, w=W, vb=v′b′, wb=w′b′),
                                                     schedule = TimeInterval(save_fields_interval),
                                                     prefix = filename * "_zonal_average",
                                                     force = true)
#=
simulation.output_writers[:averages] = JLD2OutputWriter(model, averaged_outputs,
                                                        schedule = AveragedTimeInterval(1days, window=1days, stride=1),
                                                        prefix = filename * "_averages",
                                                        verbose = true,
                                                        force = true)
=#

@info "Running the simulation..."

try
    run!(simulation, pickup=false)
catch err
    @info "run! threw an error! The error message is"
    showerror(stdout, err)
end

#####
##### Visualize
#####

using CairoMakie

fig = Figure(resolution = (2000, 1000))
ax_b = fig[1:5, 1] = LScene(fig)
ax_ζ = fig[1:5, 2] = LScene(fig)

# Extract surfaces on all 6 boundaries

iter = Node(0)
sides = keys(slicers)

zonal_file = jldopen(filename * "_zonal_average.jld2")
slice_files = NamedTuple(side => jldopen(filename * "_$(side)_slice.jld2") for side in sides)

grid = RectilinearGrid(architecture = CPU(),
                       topology = (Periodic, Bounded, Bounded),
                       size = (grid.Nx, grid.Ny, grid.Nz),
                       halo = (3, 3, 3),
                       x = (0, grid.Lx),
                       y = (0, grid.Ly),
                       z = (-grid.Lz, 0))

# Build coordinates, rescaling the vertical coordinate

xζ, yζ, zζ = nodes((Face, Face, Center), grid)
xu, yu, zu = nodes((Face, Center, Center), grid)
xv, yv, zv = nodes((Center, Face, Center), grid)
xw, yw, zw = nodes((Center, Center, Face), grid)
xb, yb, zb = nodes((Center, Center, Center), grid)

zscale = 300
zu = zu .* zscale
zb = zb .* zscale
zζ = zζ .* zscale

zonal_slice_displacement = 1.35

b_slices = (
      west = @lift(Array(slice_files.west["timeseries/b/"   * string($iter)][1, :, :])),
      east = @lift(Array(slice_files.east["timeseries/b/"   * string($iter)][1, :, :])),
     south = @lift(Array(slice_files.south["timeseries/b/"  * string($iter)][:, 1, :])),
     north = @lift(Array(slice_files.north["timeseries/b/"  * string($iter)][:, 1, :])),
    bottom = @lift(Array(slice_files.bottom["timeseries/b/" * string($iter)][:, :, 1])),
       top = @lift(Array(slice_files.top["timeseries/b/"    * string($iter)][:, :, 1]))
)

clims_b = @lift extrema(slice_files.top["timeseries/b/" * string($iter)][:])
kwargs_b = (colorrange=clims_b, colormap=:deep, show_axis=false)

surface!(ax_b, yb, zb, b_slices.west;   transformation = (:yz, xb[1]),   kwargs_b...)
surface!(ax_b, yb, zb, b_slices.east;   transformation = (:yz, xb[end]), kwargs_b...)
surface!(ax_b, xb, zb, b_slices.south;  transformation = (:xz, yb[1]),   kwargs_b...)
surface!(ax_b, xb, zb, b_slices.north;  transformation = (:xz, yb[end]), kwargs_b...)
surface!(ax_b, xb, yb, b_slices.bottom; transformation = (:xy, zb[1]),   kwargs_b...)
surface!(ax_b, xb, yb, b_slices.top;    transformation = (:xy, zb[end]), kwargs_b...)

# b_avg = @lift zonal_file["timeseries/b/" * string($iter)][1, :, :]
# u_avg = @lift zonal_file["timeseries/u/" * string($iter)][1, :, :]

clims_u = @lift extrema(zonal_file["timeseries/u/" * string($iter)][1, :, :])
clims_u = (-0.4, 0.4)

# contour!(ax_b, yb, zb, b_avg; levels = 25, color = :black, linewidth = 2, transformation = (:yz, zonal_slice_displacement * xb[end]), show_axis=false)
# surface!(ax_b, yu, zu, u_avg; transformation = (:yz, zonal_slice_displacement * xu[end]), colorrange=clims_u, colormap=:balance)

rotate_cam!(ax_b.scene, (π/24, -π/6, 0))

ζ_slices = (
      west = @lift(Array(slice_files.west["timeseries/ζ/"   * string($iter)][1, :, :])),
      east = @lift(Array(slice_files.east["timeseries/ζ/"   * string($iter)][1, :, :])),
     south = @lift(Array(slice_files.south["timeseries/ζ/"  * string($iter)][:, 1, :])),
     north = @lift(Array(slice_files.north["timeseries/ζ/"  * string($iter)][:, 1, :])),
    bottom = @lift(Array(slice_files.bottom["timeseries/ζ/" * string($iter)][:, :, 1])),
       top = @lift(Array(slice_files.top["timeseries/ζ/"    * string($iter)][:, :, 1]))
)

u_slices = (
      west = @lift(Array(slice_files.west["timeseries/u/"   * string($iter)][1, :, :])),
      east = @lift(Array(slice_files.east["timeseries/u/"   * string($iter)][1, :, :])),
     south = @lift(Array(slice_files.south["timeseries/u/"  * string($iter)][:, 1, :])),
     north = @lift(Array(slice_files.north["timeseries/u/"  * string($iter)][:, 1, :])),
    bottom = @lift(Array(slice_files.bottom["timeseries/u/" * string($iter)][:, :, 1])),
       top = @lift(Array(slice_files.top["timeseries/u/"    * string($iter)][:, :, 1]))
)

clims_ζ = @lift extrema(slice_files.top["timeseries/ζ/" * string($iter)][:])
clims_ζ = (-1f-4, 1f-4)
kwargs_ζ = (colormap=:curl, show_axis=false, colorrange=clims_ζ)
kwargs_east = (colormap=:curl, show_axis=false, colorrange=(-5f-5, 5f-5))
clims_u = @lift extrema(slice_files.top["timeseries/u/" * string($iter)][:])
kwargs_u = (colormap=:balance, show_axis=false)

surface!(ax_ζ, yζ, zζ, ζ_slices.west;   transformation = (:yz, xζ[1]),   kwargs_ζ...)
surface!(ax_ζ, yζ, zζ, ζ_slices.east;   transformation = (:yz, xζ[end]), kwargs_east...)
surface!(ax_ζ, xζ, zζ, ζ_slices.south;  transformation = (:xz, yζ[1]),   kwargs_ζ...)
surface!(ax_ζ, xζ, zζ, ζ_slices.north;  transformation = (:xz, yζ[end]), kwargs_ζ...)
surface!(ax_ζ, xζ, yζ, ζ_slices.bottom; transformation = (:xy, zζ[1]),   kwargs_ζ...)
surface!(ax_ζ, xζ, yζ, ζ_slices.top;    transformation = (:xy, zζ[end]), kwargs_ζ...)

# b_avg = @lift zonal_file["timeseries/b/" * string($iter)][1, :, :]
# u_avg = @lift zonal_file["timeseries/u/" * string($iter)][1, :, :]

clims_u = @lift extrema(zonal_file["timeseries/u/" * string($iter)][1, :, :])
clims_u = (-0.4, 0.4)

# contour!(ax_ζ, yb, zb, b_avg; levels = 25, color = :black, linewidth = 2, transformation = (:yz, zonal_slice_displacement * xb[end]), show_axis=false)
# surface!(ax_ζ, yu, zu, u_avg; transformation = (:yz, zonal_slice_displacement * xu[end]), colorrange=clims_u, colormap=:balance)

rotate_cam!(ax_ζ.scene, (π/24, -π/6, 0))

title = @lift(string("Buoyancy, vertical vorticity, and zonally-averaged u at t = ",
                     prettytime(zonal_file["timeseries/t/" * string($iter)])))

fig[0, :] = Label(fig, title, textsize=30)

display(fig)

iterations = parse.(Int, keys(zonal_file["timeseries/t"]))

record(fig, filename * ".mp4", iterations, framerate=12) do i
    @info "Plotting iteration $i of $(iterations[end])..."
    iter[] = i
end

display(fig)

for file in slice_files
    close(file)
end

close(zonal_file)


text:
using ITensors
using ITensorUnicodePlots
using Graphs

g = grid((5,))
tn = itensornetwork(g; linkspaces=10, sitespaces=2)
@visualize fig tn

fig

ext: jl


text:
export qrfactUnblocked, updateqrfactUnblocked!

# Add the non in-place version of qrfactUnblocked! defined in
# https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/qr.jl
function qrfactUnblocked(A::AbstractMatrix{T}, arg...; kwargs...) where T
    LinearAlgebra.require_one_based_indexing(A)
    AA = similar(A, LinearAlgebra._qreltype(T), size(A))
    copyto!(AA, A)
    return LinearAlgebra.qrfactUnblocked!(AA, arg...; kwargs...)
end

function updateqrfactUnblocked!(S::QR{T,Array{T,2}}, a::AbstractVector{T}) where {T}
    m, n = size(S.factors)
    # Add one entry to F.τ
    # Add one column to F.factors
    # factors = hcat(S.factors, S.Q'*a)

    factors = zeros(m, n+1)
    view(factors,:,1:n) .= S.factors
    mul!(view(factors,:,n+1), S.Q', a)
    # factors = hcat(S.factors, S.Q'*a)
    x = view(factors, n+1:m, n+1)
    τk = LinearAlgebra.reflector!(x)
    push!(S.τ, τk)

    return QR(factors, S.τ)
end

function updateqrfactUnblocked!(S::QR{T,Array{T,2}}, A::AbstractMatrix{T}) where {T}
    m, n = size(S.factors)
    mA, nA = size(A)
    @assert mA == m "Length of the columns do not match"
    # Add one entry to F.τ

    # Add one column to F.factors
    factors = zeros(m, n+nA)
    factors .= hcat(S.factors, S.Q'*A)

    @inbounds for k = n+1:min(m - 1 + !(T<:Real), n+nA)
        x = view(factors, k:m, k)
        τk = LinearAlgebra.reflector!(x)
        push!(S.τ, τk)
        LinearAlgebra.reflectorApply!(x, τk, view(factors, k:m, k + 1:n + nA))
    end
    return QR(factors, S.τ)
end

# in-place version of qrfactUnblocked! defined in
# https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/qr.jl
# function qrfactUnblocked!(A::AbstractMatrix{T}) where {T}
#     LinearAlgebra.require_one_based_indexing(A)
#     m, n = size(A)
#     τ = zeros(T, min(m,n))
#     for k = 1:min(m - 1 + !(T<:Real), n)
#         x = view(A, k:m, k)
#         τk = LinearAlgebra.reflector!(x)
#         τ[k] = τk
#         LinearAlgebra.reflectorApply!(x, τk, view(A, k:m, k + 1:n))
#     end
#         QR(A, τ)
# end

# struct UpdatableQR{T} <: Factorization{T}
#     factors::Array{Array{T,1},1}
#     τ::Array{Float64,1}
#     # cache::Array{Float64,1}
# end
#

# in-place version of qrfactUnblocked! defined in
# https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/qr.jl
# function qrfactUnblocked!(A::AbstractMatrix{T}) where {T}
#     LinearAlgebra.require_one_based_indexing(A)
#     m, n = size(A)
#     τ = zeros(T, min(m,n))
#     for k = 1:min(m - 1 + !(T<:Real), n)
#         x = view(A, k:m, k)
#         τk = LinearAlgebra.reflector!(x)
#         τ[k] = τk
#         LinearAlgebra.reflectorApply!(x, τk, view(A, k:m, k + 1:n))
#     end
#         QR(A, τ)
# end

# function updateqrfactUnblocked!(S::UpdatableQR{T}, a::AbstractVector{T}) where {T}
#     m, n = size(S.factors)
#     # Add one entry to F.τ
#
#     # Add one column to F.factors
#     # factors = hcat(S.factors, S.Q'*a)
#
#     # @inbounds for k = n+1:min(m - 1 + !(T<:Real), n+1)
#         mul!(S.cache, , a)
#         push!(S.factors, S.Q'*a)
#         # view(factors,:,1:n) .= S.factors
#         # mul!(view(factors,:,n+1), , a)
#         # factors = hcat(S.factors, S.Q'*a)
#         x = view(S.factors, n+1:m, n+1)
#         τk = LinearAlgebra.reflector!(x)
#         push!(S.τ, τk)
#         # LinearAlgebra.reflectorApply!(x, τk, view(factors, k:m, k + 1:n + 1))
#     # end
#     # return QR(factors, S.τ)
# end
# function UpdatableQR(A::AbstractMatrix{T}) where T
#     return UpdatableQR{T}
# end
#
# function updateqrfactUnblocked!(S::UpdatableQR{T}, a::AbstractVector{T}) where {T}
#     m, n = size(S.factors)
#     # Add one entry to F.τ
#     push!(S.τ, 0.0)
#
#     # Add one column to F.factors
#     push!(S.factors, S.Q'*a)
#
#     @inbounds for k = n+1:min(m - 1 + !(T<:Real), n+1)
#         x = view(factors, k:m, k)
#         τk = LinearAlgebra.reflector!(x)
#         S.τ[k] = τk
#         LinearAlgebra.reflectorApply!(x, τk, view(factors, k:m, k + 1:n + 1))
#     end
#     S = QR(factors, S.τ)
#     return S
# end

max_forks_repo_name: mleprovost/QROMP


text:
using ClinicalTrialSampleSize
using Base.Test

# Run tests

println("Running tests:")
tic()
println("Test 1:")
@time @test include("OneSampleTest.jl")
println("Test 2:")
@time @test include("TwoSampleTest.jl")
println("Test 3:")
@time @test include("CheckArgsTest.jl")
toc()


text:
using Documenter
using PencilFFTs

# This is to make sure that doctests in docstrings are executed correctly.
DocMeta.setdocmeta!(
    PencilFFTs, :DocTestSetup, :(using PencilFFTs); recursive=true)

function main()
    @time makedocs(
        modules = [PencilFFTs],
        authors = "Juan Ignacio Polanco <jipolanc@gmail.com> and contributors",
        repo = "https://github.com/jipolanco/PencilFFTs.jl/blob/{commit}{path}#L{line}",
        sitename = "PencilFFTs.jl",
        format = Documenter.HTML(
            prettyurls=true,
            canonical = "https://jipolanco.github.io/PencilFFTs.jl",
            # load assets in <head>
            assets=[
                "assets/custom.css",
                "assets/tomate.js",  # matomo code
            ],
        ),
        pages=[
            "Home" => "index.md",
            "tutorial.md",
            "More examples" => [
                "examples/in-place.md",
                "examples/gradient.md",
            ],
            "Library" => [
                "PencilFFTs.md",
                "Transforms.md",
                "PencilFFTs_timers.md",
                "Internals" => ["GlobalFFTParams.md"],
            ],
            "benchmarks.md",
        ],
        doctest=true,
        linkcheck=true,
        linkcheck_ignore=[
            # This URL is correct, but gets incorrectly flagged by linkcheck.
            "https://jipolanco.github.io/PencilArrays.jl/dev/PencilArrays/#PencilArrays.Pencils.range_local-Tuple{Union{PencilArray,%20Union{Tuple{Vararg{A,N}%20where%20N},%20AbstractArray{A,N}%20where%20N}%20where%20A%3C:PencilArray}}",
        ],
        checkdocs=:all,
    )

    nothing
end

main()

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/jipolanco/PencilFFTs.jl",
    forcepush = true,
)


text:
using Images, SIUnits.ShortUnits
using Base.Test

if VERSION.minor == 2
    macro mytest_throws(extype, arg)
        :(@test_throws $(esc(arg)))
    end
else
    macro mytest_throws(extype, arg)
        :(@test_throws $(esc(extype)) $(esc(arg)))
    end
end

B = rand(1:20,3,5)
cmap = uint8(repmat(linspace(12,255,20),1,3))
img = ImageCmap(copy(B),cmap,["colorspace" => "RGB", "pixelspacing" => [2.0, 3.0], "spatialorder" => Images.yx])
imgd = convert(Image, img)
imgd["pixelspacing"] = [2.0mm, 3.0mm]
@test eltype(img) == Uint8
@test eltype(imgd) == Uint8

# basic information
@test size(img) == (3,5)
@test size(imgd) == (3,5,3)
@test ndims(img) == 2
@test ndims(imgd) == 3

# copy/similar
A = randn(3,5,3)
imgc = copy(img)
@test imgc.data == img.data
imgc = copy(imgd, A)
@test imgc.data == A
img2 = similar(img)
@test isa(img2, ImageCmap)
@test img2.data != img.data
img2 = similar(imgd)
@test isa(img2, Image)

# ref/assign
prev = img[4]
@test prev == B[4]
img[4] = prev+1
@test img.data[4] == prev+1
@test img[4] == prev+1

# properties
@test colorspace(img) == "RGB"
@test colordim(img) == 0
@test colordim(imgd) == 3
@test timedim(img) == 0
@test limits(img) == (12,255)
@test limits(imgd) == (0,255)
@test pixelspacing(img) == [2.0, 3.0]
@test pixelspacing(imgd) == [2.0mm, 3.0mm]
@test spacedirections(img) == Vector{Float64}[[2.0, 0], [0, 3.0]]
@test spacedirections(imgd) == Vector{SIUnits.SIQuantity{Float64,1,0,0,0,0,0,0}}[[2.0mm, 0.0mm], [0.0mm, 3.0mm]]

@test sdims(img) == sdims(imgd)
@test coords_spatial(img) == coords_spatial(imgd)
@test size_spatial(img) == size_spatial(imgd)

@test storageorder(img) == Images.yx
@test storageorder(imgd) == [Images.yx, "color"]

A = rand(4,4,3)
@test colordim(A) == 3
Aimg = permutedims(convert(Image, A), [3,1,2])
@test colordim(Aimg) == 1

# sub/slice
s = sub(img, 2, 1:4)
@test ndims(s) == 2
@test sdims(s) == 2
@test size(s) == (1,4)
s = subim(img, 2, 1:4)
@test ndims(s) == 2
@test sdims(s) == 2
@test size(s) == (1,4)
s = sliceim(img, 2, 1:4)
@test ndims(s) == 1
@test sdims(s) == 1
@test size(s) == (4,)
s = sliceim(imgd, 2, 1:4, 1:3)
@test ndims(s) == 2
@test sdims(s) == 1
@test colordim(s) == 2
@test spatialorder(s) == ["x"]
s = sliceim(imgd, 2:2, 1:4, 1:3)
@test ndims(s) == 3
@test sdims(s) == 2
@test colordim(s) == 3
@test spatialorder(s) == ["y","x"]

# reslicing
D = randn(3,5,4)
sd = SliceData(D, 2)
C = slice(D, sd, 2)
@test C == reshape(D[1:end, 2, 1:end], size(C))
reslice!(C, sd, 3)
@test C == reshape(D[1:end, 3, 1:end], size(C))
sd = SliceData(D, 3)
C = slice(D, sd, 2)
@test C == reshape(D[1:end, 1:end, 2], size(C))

sd = SliceData(imgd, 2)
s = sliceim(imgd, sd, 2)
@test colordim(s) == 2
@test colorspace(s) == "RGB"
@test spatialorder(s) == ["y"]
@test s.data == reshape(imgd[:,2,:], size(s))
sd = SliceData(imgd, 3)
s = sliceim(imgd, sd, 2)
@test colordim(s) == 0
@test colorspace(s) == "Unknown"
@test spatialorder(s) == Images.yx
@test s.data == imgd[:,:,2]
reslice!(s, sd, 3)
@test s.data == imgd[:,:,3]

# named indexing
@test dimindex(imgd, "color") == 3
@test dimindex(imgd, "y") == 1
@test dimindex(imgd, "z") == 0
imgdp = permutedims(imgd, [3,1,2])
@test dimindex(imgdp, "y") == 2
@test img["y", 2, "x", 4] == B[2,4]
@test img["x", 4, "y", 2] == B[2,4]
chan = imgd["color", 2]
Blookup = reshape(cmap[B[:],1], size(B))
@test chan == Blookup

sd = SliceData(imgd, "x")
s = sliceim(imgd, sd, 2)
@test spatialorder(s) == ["y"]
@test s.data == reshape(imgd[:,2,:], size(s))
sd = SliceData(imgd, "y")
s = sliceim(imgd, sd, 2)
@test spatialorder(s) == ["x"]
@test s.data == reshape(imgd[2,:,:], size(s))
sd = SliceData(imgd, "x", "y")
s = sliceim(imgd, sd, 2, 1)
@test s.data == reshape(imgd[1,2,:], 3)

# spatial order, width/height, and permutations
@test spatialpermutation(Images.yx, imgd) == [1,2]
@test widthheight(imgd) == (5,3)
C = convert(Array, imgd)
@test C == imgd.data
imgd["spatialorder"] = ["x", "y"]
@test spatialpermutation(Images.xy, imgd) == [1,2]
@test widthheight(imgd) == (3,5)
C = convert(Array, imgd)
@test C == permutedims(imgd.data, [2,1,3])
imgd.properties["spatialorder"] = ["y", "x"]
@test spatialpermutation(Images.xy, imgd) == [2,1]
imgd.properties["spatialorder"] = ["x", "L"]
@test spatialpermutation(Images.xy, imgd) == [1,2]
imgd.properties["spatialorder"] = ["L", "x"]
@test spatialpermutation(Images.xy, imgd) == [2,1]
@test spatialpermutation(Images.xy, A) == [2,1]
@test spatialpermutation(Images.yx, A) == [1,2]

imgd.properties["spatialorder"] = Images.yx
imgp = permutedims(imgd, ["x", "y", "color"])
@test imgp.data == permutedims(imgd.data, [2,1,3])
imgp = permutedims(imgd, ("color", "x", "y"))
@test imgp.data == permutedims(imgd.data, [3,2,1])
@test pixelspacing(imgp) == [3.0mm, 2.0mm]
imgc = copy(imgd)
imgc["spacedirections"] = spacedirections(imgc)
delete!(imgc, "pixelspacing")
imgp = permutedims(imgc, ["x", "y", "color"])
@test spacedirections(imgp) == Vector{SIUnits.SIQuantity{Float64,1,0,0,0,0,0,0}}[[0.0mm, 3.0mm],[2.0mm, 0.0mm]]
@test pixelspacing(imgp) == [3.0mm, 2.0mm]

## Convenience constructors
# colorim
@test colordim(colorim(rand(Uint8, 3, 5, 5))) == 1
@test colordim(colorim(rand(Uint8, 5, 5, 3))) == 3
@test spatialorder(colorim(rand(Uint8, 3, 5, 5))) == ["x", "y"]
@test spatialorder(colorim(rand(Uint8, 5, 5, 3))) == ["y", "x"]
@test colordim(colorim(rand(Uint8, 4, 5, 5), "RGBA")) == 1
@test colordim(colorim(rand(Uint8, 4, 5, 5), "ARGB")) == 1
@test colordim(colorim(rand(Uint8, 5, 5, 4), "RGBA")) == 3
@test colordim(colorim(rand(Uint8, 5, 5, 4), "ARGB")) == 3
@test spatialorder(colorim(rand(Uint8, 4, 5, 5), "ARGB")) == ["x", "y"]
@test spatialorder(colorim(rand(Uint8, 5, 5, 4), "ARGB")) == ["y", "x"]
@mytest_throws ErrorException colorim(rand(Uint8, 3, 5, 3))
@mytest_throws ErrorException colorim(rand(Uint8, 4, 5, 5))
@mytest_throws ErrorException colorim(rand(Uint8, 5, 5, 4))
@mytest_throws ErrorException colorim(rand(Uint8, 4, 5, 4), "ARGB")
@mytest_throws ErrorException colorim(rand(Uint8, 5, 5, 5), "foo")
@mytest_throws ErrorException colorim(rand(Uint8, 2, 2, 2), "bar")


text:
module COBREXA

using Dates
using Distributed
using DistributedData
using JSON
using JuMP
using LinearAlgebra
using MAT
using MacroTools
using OrderedCollections
using Random
using SparseArrays
using Statistics

import Base: findfirst, getindex, show
import Pkg
import SBML # conflict with Reaction struct name

include("banner.jl")
_print_banner()

# autoloading
const _inc(path...) = include(joinpath(path...))
const _inc_all(dir) = _inc.(joinpath.(dir, filter(fn -> endswith(fn, ".jl"), readdir(dir))))

_inc_all.(
    joinpath.(
        @__DIR__,
        [
            joinpath("base", "types", "abstract"),
            joinpath("base", "logging"),
            joinpath("base", "types"),
            "base",
            "io",
            joinpath("io", "show"),
            "reconstruction",
            joinpath("analysis", "modifications"),
            "analysis",
            joinpath("analysis", "sampling"),
            joinpath("base", "utils"),
        ],
    ),
)

# export everything that isn't prefixed with _ (inspired by JuMP.jl, thanks!)
for sym in names(@__MODULE__, all = true)
    if sym in [Symbol(@__MODULE__), :eval, :include] || startswith(string(sym), ['_', '#'])
        continue
    end
    @eval export $sym
end

end # module


text:
mutable struct UpperBoundNonStochastic{S,A,O}
    upper_bound_act::Vector{A}
    upper_bound_memo::Vector{Float64}

    # Constructor
    function UpperBoundNonStochastic{S,A,O}(pomdp::POMDP{S,A,O}) where {S,A,O}
        this = new()
        # this executes just once per problem run
        this.upper_bound_act = Array{A}(n_states(pomdp))    # upper_bound_act
        this.upper_bound_memo = Array{Float64}(n_states(pomdp)) # upper_bound_memo
        return this
    end
end

fringe_upper_bound{S,A,O}(pomdp::POMDP{S,A,O}, state::S) =
    error("$(typeof(pomdp)) does not implement fringe_upper_bound")

function init_bound{S,A,O}(ub::UpperBoundNonStochastic{S,A,O},
                    pomdp::POMDP{S,A,O},
                    config::DESPOTConfig)

# function DESPOT.init_upper_bound{S,A,O}(ub::UpperBoundNonStochastic{S,A,O},
#                     pomdp::POMDP{S,A,O},
#                     config::DESPOTConfig)

    current_level_ub_memo = Array{Float64}(n_states(pomdp))
    next_level_ub_memo = Array{Float64}(n_states(pomdp))

    next_state = S()
    r::Float64 = 0.0
    trans_distribution = create_transition_distribution(pomdp)
    rng = DESPOTRandomNumber(0) # dummy RNG

    fill!(current_level_ub_memo, -Inf)

    for s in iterator(states(pomdp))
        next_level_ub_memo[state_index(pomdp,s)+1] = fringe_upper_bound(pomdp,s) # 1-based indexing
    end

    for i in 1:config.search_depth # length of horizon
        for s in iterator(states(pomdp))
            for a in iterator(actions(pomdp))
                trans_distribution.state = s
                trans_distribution.action = a
                next_state = POMDPs.rand(rng, trans_distribution, next_state)
                r = reward(pomdp, s, a)
                possibly_improved_value =
                    r + pomdp.discount * next_level_ub_memo[state_index(pomdp,next_state)+1]
                if (possibly_improved_value > current_level_ub_memo[state_index(pomdp,s)+1])
                    current_level_ub_memo[state_index(pomdp,s)+1] = possibly_improved_value
                    if i == config.search_depth
                        # Set best actions when last level is being computed
                        ub.upper_bound_act[state_index(pomdp,s)+1] = a
                    end
                end
            end # for a
        end # for s

        # swap array references
        tmp = current_level_ub_memo
        current_level_ub_memo = next_level_ub_memo
        next_level_ub_memo = tmp

        fill!(current_level_ub_memo,-Inf)
    end

    #TODO: this can probably be done more optimally (by referencing upper_bound_memo to start with),
    # however, this only runs once per problem and is probably not a big deal. Leave it as is for now.
    copy!(ub.upper_bound_memo, next_level_ub_memo)

    return nothing
end

function upper_bound{S,A,O}(ub::UpperBoundNonStochastic{S,A,O},
                     pomdp::POMDP{S,A,O},
                     particles::Vector{DESPOTParticle{S}},
                     config::DESPOTConfig)
  weight_sum = 0.
  total_cost = 0.

  for p in particles
    weight_sum += p.weight
    total_cost += p.weight * ub.upper_bound_memo[state_index(pomdp,p.state)+1]
  end
  return total_cost / weight_sum
end


text:
using Multilane
using MCTS
using POMDPToolbox
using POMDPs
# using POMCP
using Missings
using DataFrames
using CSV

@everywhere using Missings
@everywhere using Multilane
@everywhere using POMDPToolbox

using Gallium

dpws = DPWSolver(depth=20,
                 n_iterations=1_000_000,
                 max_time=1.0,
                 exploration_constant=50.0,
                 k_state=4.0,
                 alpha_state=1/8,
                 enable_action_pw=false,
                 estimate_value=RolloutEstimator(SimpleSolver()))

solver = SingleBehaviorSolver(dpws, Multilane.NORMAL)

@show N = 500
alldata = DataFrame()

for lambda in 2.0.^(-2:5)
# for lambda in 2.^1
# for lambda in [1.0, 10.0, 100.0, 1000.0]

    @show lambda

    behaviors = standard_uniform(correlation=0.75)
    pp = PhysicalParam(4, lane_length=100.0)
    dmodel = NoCrashIDMMOBILModel(10, pp,
                                  behaviors=behaviors,
                                  p_appear=1.0,
                                  lane_terminate=true,
                                  brake_terminate_thresh=4.0,
                                  max_dist=1000.0
                                 )
    rmodel = SuccessReward(lambda=lambda)
    pomdp = NoCrashPOMDP{typeof(rmodel)}(dmodel, rmodel, 0.95, false)
    mdp = NoCrashMDP{typeof(rmodel)}(dmodel, rmodel, 0.95, false)
    planner = solve(solver, mdp)

    sim_pomdp = deepcopy(pomdp)
    sim_pomdp.throw=true

    sims = []

    for i in 1:N
        rng_seed = i+40000
        rng = MersenneTwister(rng_seed)
        agg_up = AggressivenessUpdater(sim_pomdp, 500, 0.1, 0.1, WeightUpdateParams(smoothing=0.0, wrong_lane_factor=0.5), MersenneTwister(rng_seed+50000))
        metadata = Dict(:rng_seed=>rng_seed,
                        :lambda=>lambda,
                        :solver=>"baseline",
                        :dt=>pp.dt
                   )

        is = initial_state(pomdp, MersenneTwister(i+30000))
        ips = MLPhysicalState(is)

        hr = HistoryRecorder(max_steps=100, rng=rng, capture_exception=true)
        push!(sims, Sim(sim_pomdp, planner, agg_up, ips, is,
                        simulator=hr,
                        metadata=metadata
                       ))
        @assert problem(last(sims)).throw
    end

    # data = run(sims) do sim, hist
    data = run_parallel(sims) do sim, hist

        if isnull(exception(hist))
            p = problem(sim)
            steps_in_lane = 0
            steps_to_lane = missing
            nb_brakes = 0
            crashed = false
            min_speed = Inf
            min_ego_speed = Inf
            for (k,(s,sp)) in enumerate(eachstep(hist, "s,sp"))

                nb_brakes += detect_braking(p, s, sp)

                if sp.cars[1].y == p.rmodel.target_lane
                    steps_in_lane += 1
                end
                if sp.cars[1].y == p.rmodel.target_lane
                    if ismissing(steps_to_lane)
                        steps_to_lane = k
                    end
                end

                if is_crash(p, s, sp)
                    crashed = true
                end

                min_speed = min(minimum(c.vel for c in sp.cars), min_speed)
                min_ego_speed = min(min_ego_speed, sp.cars[1].vel)
            end
            time_to_lane = steps_to_lane*p.dmodel.phys_param.dt
            distance = last(state_hist(hist)).x

            return [:n_steps=>n_steps(hist),
                    :mean_iterations=>mean(ai[:tree_queries] for ai in eachstep(hist, "ai")),
                    :reward=>discounted_reward(hist),
                    :crashed=>crashed,
                    :steps_to_lane=>steps_to_lane,
                    :steps_in_lane=>steps_in_lane,
                    :nb_brakes=>nb_brakes,
                    :exception=>false,
                    :distance=>distance,
                    :mean_ego_speed=>distance/(n_steps(hist)*p.dmodel.phys_param.dt),
                    :min_speed=>min_speed,
                    :min_ego_speed=>min_ego_speed,
                    :terminal=>string(get(last(state_hist(hist)).terminal, missing))
                   ]
        else
            warn("Error in Simulation")
            showerror(STDERR, get(exception(hist)))
            return [:exception=>true,
                    :ex_type=>string(typeof(get(exception(hist))))
                   ]
        end
    end

    success = 100.0*sum(data[:terminal].=="lane")/N
    brakes = 100.0*sum(data[:nb_brakes].>=1)/N
    @printf("%% reaching:%5.1f; %% braking:%5.1f\n", success, brakes)

    @show extrema(data[:distance])
    @show mean(data[:mean_ego_speed])
    @show minimum(data[:min_speed])

    if isempty(alldata)
        alldata = data
    else
        alldata = vcat(alldata, data)
    end
end

# @show alldata

datestring = Dates.format(now(), "E_d_u_HH_MM")
filename = Pkg.dir("Multilane", "data", "baseline_curve_"*datestring*".csv")
println("Writing data to $filename")
CSV.write(filename, alldata)


text:
@testset "Instruction Tests" begin
    model = ALifeBenchmark.TierraModel(
        ALifeBenchmark.TierrianInstruction[
        ALifeBenchmark.NOP_0()
        ])

    ancestor = model.organisms[model.organism_keys[1]]

    # NOP
    ALifeBenchmark.apply!(ALifeBenchmark.NOP_0(), ancestor, model)
    @test sum(model.memory) == 0
    ALifeBenchmark.apply!(ALifeBenchmark.NOP_1(), ancestor, model)
    @test sum(model.memory) == 0

    # OR
    ALifeBenchmark.apply!(ALifeBenchmark.OR_1(), ancestor, model)
    @test ancestor.c == 1
    ALifeBenchmark.apply!(ALifeBenchmark.OR_1(), ancestor, model)
    @test ancestor.c == 1

    # Shift
    ALifeBenchmark.apply!(ALifeBenchmark.SH_L(), ancestor, model)
    @test ancestor.c == 2
    ALifeBenchmark.apply!(ALifeBenchmark.SH_L(), ancestor, model)
    @test ancestor.c == 4
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
function coverage_robot(inputs::Array{Float64})
    # TODO: visualize why this doesn't work
    # wall or covered spot in front
    # green = RoboGrid.color_to_float("green")
    # if inputs[3] == 1.0 || inputs[4] == green
    #     return rand([3, 4])
    # end
    # # paint green in self spot
    # if inputs[10] != "green"
    #     return 6
    # end
    # 1
    rand(1:4)
end

no_terminate(e::Episode) = false

@testset "Exit reward" begin
    e = Episode(Grid(), reward=RoboGrid.exit_reward)
    total_reward = run!(e, coverage_robot; terminate=no_terminate)
    println("Found exit in ", e.meta["step"], " steps, ", e.meta["time"], " s, ",
            e.meta["memory"], " bytes")
    @test total_reward == 1.0
end

@testset "Food reward" begin
    g = Grid()
    foodind = findfirst([RoboGrid.object(c) == "food" for c in g.cells])
    g.exits[1] = LinearIndices(g.cells)[foodind]
    e = Episode(g, reward=RoboGrid.food_reward)
    total_reward = run!(e, coverage_robot; terminate=no_terminate)
    println("Found food in ", e.meta["step"], " steps, ", e.meta["time"], " s, ",
            e.meta["memory"], " bytes")
    @test total_reward == 1.0
    @test sum(map(c->RoboGrid.object(c) == "food", e.grid.cells)) == 0
end

@testset "Coverage reward" begin
    e = Episode(Grid(), reward=RoboGrid.coverage_reward)
    total_reward = run!(e, coverage_robot; terminate=no_terminate)
    println("Found exit in ", e.meta["step"], " steps, ", e.meta["time"], " s, ",
            e.meta["memory"], " bytes")
    visited = sum(e.meta["visited"])
    println("Covered $visited cells for $total_reward")
    @test total_reward >= 1
    @test visited == total_reward
end


text:
Base.show(io::IO, j::Object) = _show(io, j)
# use the Base fallback AbstractArray show method instead
# Base.show(io::IO, j::Array) = _show(io, j)
_show(io::IO, x, indent=0, offset=0) = show(io, x)

function _show(io::IO, obj::Object, indent=0, offset=0)
    if isempty(obj)
        print(io, "{}")
        return
    end
    println(io, "{")
    indent += 1
    keyvals = collect(obj)
    keys = map(x->x[1], keyvals)
    vals = map(x->x[2], keyvals)
    maxlen = maximum(map(sizeof, keys)) + 5
    # @show maxlen
    for i = 1:length(keys)
        Base.write(io, "  "^indent)
        Base.write(io, lpad("\"$(keys[i])\"" * ": ", maxlen + offset, ' '))
        _show(io, vals[i], indent, maxlen + offset)
        if i == length(keys)
            indent -= 1
            Base.write(io, "\n" * ("  "^indent * " "^offset) * "}")
        else
            Base.write(io, ",\n")
        end
    end
    return
end

function _show(io::IO, arr::Array, indent=0, offset=0)
    if isempty(arr)
        print(io, "[]")
        return
    end
    println(io, "[")
    indent += 1
    vals = collect(arr)
    for (i, val) in enumerate(vals)
        Base.write(io, "  "^indent * " "^offset)
        _show(io, vals[i], indent, offset)
        if i == length(vals)
            indent -= 1
            Base.write(io, "\n" * ("  "^indent * " "^offset) * "]")
        else
            Base.write(io, ",\n")
        end
    end
    return
end


text:
using DataFrames, CSV, McCormick, EAGO, JuMP

include("nitrification_surrogate.jl")

# alternatively McCormick can be used instead EAGO in this example
include("nitrification_surrogate.jl")

function Ammonia_oxidation_rate(x)
    V = 1000.0
    m_in = V / 240 # L/s
    m_out = V / 240 # L/s
    X_AO = 0.505 # g VSS/L
    M_N = 14.0067 #
    c_NH = x[1]
    c_O = x[2]
    r_AOmax = 0.67 # mg N-NH+ /g VSS_AO min
    K_OAO = 0.3 # mg O2/L
    K_SAO = 0.34 # mg N-NH4+/L
    K_IAO = 6200 # mg N-NH4+/L
    r_AO = r_AOmax*(c_NH/(K_SAO + c_NH + c_NH^2/K_IAO))*(c_O/(K_OAO + c_O)); # mg N-NH+ /g VSS_AO min
    return r_AO
end

## Kinetics Model
function DataModel(t, y, p, u)

    ## Parameters
    K_SAO = 0.24 # mg N-NH/L
    K_IAO = 6200.0 # mg N-NH/L
    K_OAO = 0.3 # mg/L
    K_SNO = 1.6 # mg N-NI/L
    K_INO = 13000.0 # mg N-NI/L
    K_ONO = 1.5 # mg/L
    X_AO = 0.505 # g VSS/L
    X_NO = 0.151 # g VSS/L
    rNO_MAX = 1.07/60 # mg N-NI/(g-VSS_NO*s)

    # Volumetric flow
    V = 1000.0 # L
    m_in = V/240 # L/s
    m_out = V/240 # L/s
    Cs = 31.0 # mg N-NH/L

    # Oxygen flow
    SOTE = 0.1
    cO_s = 9.1 # mg/L
    Phi_OAO = 2.5 # mg O2/ mg N-NH
    Phi_ONO = 0.32 # mg O2/ mg N-NI

    cNH, cNI, cNA, cO = y

    if t >= 3080.0 && t <= 3100.0
        Cin = p[1]
    else
        Cin = Cs
    end

    # ANN
    #@show cNH, cO
    rAO = ANN1([cNH, cO])/60.0 # mg N-NH/(g-VSS_AO*s)
    #rAO = Ammonia_oxidation_rate([cNH, cO])
    t1 =(cNI/positive(K_SNO + cNI + cNI^2 /K_INO))
    t2 =(cO/positive(K_ONO + cO))
    rNO = rNO_MAX*t1*t2

    # Reaction Rate
    R_NH = -rAO*X_AO  # mg N-NH/(L*s)
    R_NI = rAO*X_AO - rNO*X_NO  # mg N-NI/(L*s)
    R_NA = rNO*X_NO  # mg N-NA/(L*s)

    W = 0.2967*u[1] # mg/s
    SOTR = SOTE*W # mg/s
    kla = SOTR/cO_s/V # s^(-1)
    R_O = kla*(cO_s - cO) - rAO*Phi_OAO*X_AO - rNO*Phi_ONO*X_NO # mg O2/(L*s)

    # Mass Balacne
    h1 = 1 /V*(m_in*Cin - m_out*cNH) + R_NH
    h2 = R_NI
    h3 = R_NA
    h4 = R_O

    return h1, h2, h3, h4
end

t0 = 0.0
h = 20.0
tspan_0 = 3080.0
n_out_0 = convert(Int, round(tspan_0 / h))

function EE_Form_Initial(u0)
    cNH0 = 30.0
    cNI0 = 0.0
    cNA0 = 0.0
    cO0 = 0.0
    y0 = cNH0, cNI0, cNA0, cO0
    p0 = 40.0
    y = y0
    y1, y2, y3, y4 = y0
    t = t0
    for j = 1:n_out_0
        t = t + h
        dy1, dy2, dy3, dy4 = DataModel(t, y, p0, u0)    # Perform Step
        #y1 = bnd(y1 + h*dy1, 0.0, 40.0)  # Enforce physical bounds on quantity at each step
        #y2 = bnd(y2 + h*dy2, 0.0, 40.0)  # Enforce physical bounds on quantity at each step
        #y3 = bnd(y3 + h*dy3, 0.0, 40.0)  # Enforce physical bounds on quantity at each step
        #y4 = bnd(y4 + h*dy4, 0.0, 9.1)
        y1 = y1 + h*dy1
        y2 = y2 + h*dy2
        y3 = y3 + h*dy3
        y4 = y4 + h*dy4
        y = y1, y2, y3, y4              # Repackage as tuple
    end
    return t, y
end

u0 = 440.0
tc, yc = EE_Form_Initial(u0)
# 500, 50 bad....
#h = 0.1
tspan = 500.0
n_out = convert(Int, round(tspan / h))
function EE_Form(p, u)
    y1, y2, y3, y4 = yc
    y = y1, y2, y3, y4
    t = tc
    for j = 1:n_out
            t = t + h
            #@show j
            dy1, dy2, dy3, dy4 = DataModel(t, y, p, u)
            #y1 = bnd(y1 + h*dy1, 0.0, 40.0)  # Enforce physical bounds on quantity at each step
            #y2 = bnd(y2 + h*dy2, 0.0, 40.0)  # Enforce physical bounds on quantity at each step
            #y3 = bnd(y3 + h*dy3, 0.0, 40.0)  # Enforce physical bounds on quantity at each step
            #y4 = bnd(y4 + h*dy4, 0.0, 9.1)
            #y1 = positive(y1 + h*dy1)
            #y2 = positive(y2 + h*dy2)
            #y3 = positive(y3 + h*dy3)
            #y4 = positive(y4 + h*dy4)
            y1 = y1 + h*dy1
            y2 = y2 + h*dy2
            y3 = y3 + h*dy3
            y4 = positive(y4 + h*dy4)
            y = y1, y2, y3, y4              # Repackage as tuple
            #@show y
    end
    return t, y
end

#=
function gsip(p, u)
    @show p, u
    t, y = EE_Form(p, u)
    cNH_f, cNI_f, cNA_f, cO_f = y
    @show cNH_f - SP
    return cNH_f - SP
end
=#
function gsip(u,p)
    SP = 30.0
    y1, y2, y3, y4 = yc
    y = y1, y2, y3, y4
    t = tc
    for j = 1:n_out
        t = t + h
        dy1, dy2, dy3, dy4 = DataModel(t, y, p, u)
        y1 = y1 + h*dy1
        y2 = y2 + h*dy2
        y3 = y3 + h*dy3
        y4 = positive(y4 + h*dy4)
        y = y1, y2, y3, y4              # Repackage as tuple
    end
    cNH_f, cNI_f, cNA_f, cO_f = y
    return cNH_f - SP
end

f(u) = tspan * u[1]

function solve_nitrification_sip()

    sip_x_lo = Float64[440.0]
    sip_x_hi = Float64[2000.0]
    sip_p_lo = Float64[31.0]
    sip_p_hi = Float64[40.0]

    opt = EAGO.Optimizer
    sip_result = explicit_sip_solve(
        sip_x_lo,
        sip_x_hi,
        sip_p_lo,
        sip_p_hi,
        f,
        [gsip],
        sip_absolute_tolerance = 3000.00, # a little less than 1E-2 rel tolerance (but not significant)
        sip_optimizer = opt,
        upper_bounding_depth = 12,
        verbosity = 0
    )
    return sip_result
end

polymer_sip_result = solve_nitrification_sip()

max_stars_repo_name: PSORLab/RobustHybridModels

max_issues_repo_name: PSORLab/RobustHybridModels

max_forks_repo_name: PSORLab/RobustHybridModels


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

struct NetworkWatchersApi <: SwaggerApi
    client::Swagger.Client
end

"""

Verifies the possibility of establishing a direct TCP connection from a virtual machine to a given endpoint including another VM or an arbitrary remote server.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: parameters::ConnectivityParameters (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ConnectivityInformation
"""
function networkWatchersCheckConnectivity(_api::NetworkWatchersApi, resourceGroupName::String, networkWatcherName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", ConnectivityInformation, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/connectivityCheck", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Creates or updates a network watcher in the specified resource group.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: parameters::NetworkWatcher (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NetworkWatcher
"""
function networkWatchersCreateOrUpdate(_api::NetworkWatchersApi, resourceGroupName::String, networkWatcherName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", NetworkWatcher, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Deletes the specified network watcher resource.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function networkWatchersDelete(_api::NetworkWatchersApi, resourceGroupName::String, networkWatcherName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets the specified network watcher by resource group.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NetworkWatcher
"""
function networkWatchersGet(_api::NetworkWatchersApi, resourceGroupName::String, networkWatcherName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", NetworkWatcher, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Queries status of flow log on a specified resource.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: parameters::FlowLogStatusParameters (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: FlowLogInformation
"""
function networkWatchersGetFlowLogStatus(_api::NetworkWatchersApi, resourceGroupName::String, networkWatcherName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", FlowLogInformation, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/queryFlowLogStatus", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets the next hop from the specified VM.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: parameters::NextHopParameters (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NextHopResult
"""
function networkWatchersGetNextHop(_api::NetworkWatchersApi, resourceGroupName::String, networkWatcherName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", NextHopResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/nextHop", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets the current network topology by resource group.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: parameters::TopologyParameters (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Topology
"""
function networkWatchersGetTopology(_api::NetworkWatchersApi, resourceGroupName::String, networkWatcherName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Topology, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/topology", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Initiate troubleshooting on a specified resource
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: parameters::TroubleshootingParameters (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: TroubleshootingResult
"""
function networkWatchersGetTroubleshooting(_api::NetworkWatchersApi, resourceGroupName::String, networkWatcherName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", TroubleshootingResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/troubleshoot", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Get the last completed troubleshooting result on a specified resource
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: parameters::QueryTroubleshootingParameters (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: TroubleshootingResult
"""
function networkWatchersGetTroubleshootingResult(_api::NetworkWatchersApi, resourceGroupName::String, networkWatcherName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", TroubleshootingResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/queryTroubleshootResult", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets the configured and effective security group rules on the specified VM.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: parameters::SecurityGroupViewParameters (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: SecurityGroupViewResult
"""
function networkWatchersGetVMSecurityRules(_api::NetworkWatchersApi, resourceGroupName::String, networkWatcherName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", SecurityGroupViewResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/securityGroupView", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets all network watchers by resource group.
Param: resourceGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NetworkWatcherListResult
"""
function networkWatchersList(_api::NetworkWatchersApi, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", NetworkWatcherListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets all network watchers by subscription.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NetworkWatcherListResult
"""
function networkWatchersListAll(_api::NetworkWatchersApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", NetworkWatcherListResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Network/networkWatchers", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Configures flow log on a specified resource.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: parameters::FlowLogInformation (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: FlowLogInformation
"""
function networkWatchersSetFlowLogConfiguration(_api::NetworkWatchersApi, resourceGroupName::String, networkWatcherName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", FlowLogInformation, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/configureFlowLog", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Verify IP flow from the specified VM to a location given the currently configured NSG rules.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: parameters::VerificationIPFlowParameters (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: VerificationIPFlowResult
"""
function networkWatchersVerifyIPFlow(_api::NetworkWatchersApi, resourceGroupName::String, networkWatcherName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", VerificationIPFlowResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/ipFlowVerify", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export networkWatchersCheckConnectivity, networkWatchersCreateOrUpdate, networkWatchersDelete, networkWatchersGet, networkWatchersGetFlowLogStatus, networkWatchersGetNextHop, networkWatchersGetTopology, networkWatchersGetTroubleshooting, networkWatchersGetTroubleshootingResult, networkWatchersGetVMSecurityRules, networkWatchersList, networkWatchersListAll, networkWatchersSetFlowLogConfiguration, networkWatchersVerifyIPFlow


text:
# ---
# title: New columns on the fly
# cover: assets/new_columns_on_the_fly.png
# description: Transforming data before generating the plot
# author: "[Pietro Vertechi](https://github.com/piever)"
# ---

using AlgebraOfGraphics, CairoMakie
set_aog_theme!() #src

# Use a `Tuple` to pass combine several columns into a unique operation.

df = (x=rand(100), y=rand(100), z=rand(100), c=rand(["a", "b"], 100))
layers = linear() + mapping(color=:z)
plt = data(df) * layers * mapping(:x, (:x, :y, :z) => (+) => "x + y + z", layout=:c)
fg = draw(plt)

# save cover image #src
mkpath("assets") #src
save("assets/new_columns_on_the_fly.png", fg) #src


text:
module Constants

export JNI_VERSION_1_1, JNI_VERSION_1_2, JNI_VERSION_1_4, JNI_VERSION_1_6, JNI_VERSION_1_8,
    JNI_TRUE, JNI_FALSE,
    JNI_OK, JNI_ERR, JNI_EDETACHED, JNI_EVERSION, JNI_ENOMEM, JNI_EEXIST, JNI_EINV

const JNI_VERSION_1_1 = convert(Cint, 0x00010001)
const JNI_VERSION_1_2 = convert(Cint, 0x00010002)
const JNI_VERSION_1_4 = convert(Cint, 0x00010004)
const JNI_VERSION_1_6 = convert(Cint, 0x00010006)
const JNI_VERSION_1_8 = convert(Cint, 0x00010008)
const JNI_VERSION_9   = convert(Cint, 0x00090000)
const JNI_VERSION_10  = convert(Cint, 0x000a0000)

const JNI_TRUE = convert(Cuchar, 1)
const JNI_FALSE = convert(Cuchar, 0)

# Return Values
const JNI_OK           = convert(Cint, 0)               #/* success */
const JNI_ERR          = convert(Cint, -1)              #/* unknown error */
const JNI_EDETACHED    = convert(Cint, -2)              #/* thread detached from the VM */
const JNI_EVERSION     = convert(Cint, -3)              #/* JNI version error */
const JNI_ENOMEM       = convert(Cint, -4)              #/* not enough memory */
const JNI_EEXIST       = convert(Cint, -5)              #/* VM already created */
const JNI_EINVAL       = convert(Cint, -6)              #/* invalid arguments */

end


text:
using Dates
using Downloads: download
using TimeZones: DEPS_DIR

const LATEST_FILE = joinpath(DEPS_DIR, "latest")
const LATEST_FORMAT = Dates.DateFormat("yyyy-mm-ddTHH:MM:SS")
const LATEST_DELAY = Hour(1)  # In 1996 a correction to a release was made an hour later

function read_latest(io::IO)
    version = readline(io)
    retrieved_utc = DateTime(readline(io), LATEST_FORMAT)
    return version, retrieved_utc
end

function read_latest(filename::AbstractString)
    open(filename, "r") do io
        read_latest(io)
    end
end

function write_latest(io::IO, version::AbstractString, retrieved_utc::DateTime=now(Dates.UTC))
    write(io, version)
    write(io, "\n")
    write(io, Dates.format(retrieved_utc, LATEST_FORMAT))
end

const LATEST = let T = Tuple{AbstractString, DateTime}
    isfile(LATEST_FILE) ? Ref{T}(read_latest(LATEST_FILE)) : Ref{T}()
end

function set_latest_cached(version::AbstractString, retrieved_utc::DateTime=now(Dates.UTC))
    LATEST[] = version, retrieved_utc
    open(LATEST_FILE, "w") do io
        write_latest(io, version, retrieved_utc)
    end
end

function latest_cached(now_utc::DateTime=now(Dates.UTC))
    if isassigned(LATEST)
        latest_version, latest_retrieved_utc = LATEST[]

        if now_utc - latest_retrieved_utc < LATEST_DELAY
            return latest_version
        end
    end

    return nothing
end

"""
    tzdata_versions() -> Vector{String}

Retrieves all of the currently available tzdata versions from IANA. The version list is
ordered from earliest to latest.

# Examples
```julia
julia> last(tzdata_versions())  # Current latest available tzdata version
"2020a"
```
"""
function tzdata_versions()
    releases_file = download("https://data.iana.org/time-zones/releases/")

    html = try
        read(releases_file, String)
    finally
        rm(releases_file)
    end

    versions = [
        m[:version]
        for m in eachmatch(r"href=\"tzdata(?<version>(?:\d{2}){1,2}[a-z]?).tar.gz\"", html)
    ]

    # Correctly order releases which include 2-digit years (e.g. "96")
    sort!(versions, by=v -> length(v) < 5 ? "19$v" : v)

    return versions
end

"""
    tzdata_url(version="latest") -> AbstractString

Generates a HTTPS URL for the specified tzdata version. Typical version strings are
formatted as 4-digit year followed by a lowercase ASCII letter. Available versions start
with "tzdata" and are listed on "https://data.iana.org/time-zones/releases/" or
"ftp://ftp.iana.org/tz/releases/".

# Examples
```julia
julia> tzdata_url("2017a")
"https://data.iana.org/time-zones/releases/tzdata2017a.tar.gz"
```
"""
function tzdata_url(version::AbstractString="latest")
    # Note: We could also support FTP but the IANA server is unreliable and likely
    # to break if working from behind a firewall.
    if version == "latest"
        "https://data.iana.org/time-zones/tzdata-latest.tar.gz"
    else
        "https://data.iana.org/time-zones/releases/tzdata$version.tar.gz"
    end
end

"""
    tzdata_download(version="latest", dir=tempdir()) -> AbstractString

Downloads a tzdata archive from IANA using the specified `version` to the specified
directory. See `tzdata_url` for details on tzdata version strings.
"""
function tzdata_download(version::AbstractString="latest", dir::AbstractString=tempdir())
    now_utc = now(Dates.UTC)
    if version == "latest"
        latest_version = latest_cached(now_utc)
        if latest_version !== nothing
            archive = joinpath(dir, "tzdata$(latest_version).tar.gz")
            isfile(archive) && return archive
        end
    end

    url = tzdata_url(version)
    archive = download(url, joinpath(dir, basename(url)))  # Overwrites the local file if any

    # Note: An "HTTP 404 Not Found" may result in the 404 page being downloaded. Also,
    # catches issues with corrupt archives
    if !isarchive(archive)
        rm(archive)
        error("Unable to download $version tzdata")
    end

    # Rename the file to have an explicit version
    if version == "latest"
        version = tzdata_version_archive(archive)

        archive_versioned = joinpath(dir, "tzdata$version.tar.gz")
        mv(archive, archive_versioned, force=true)
        archive = archive_versioned

        set_latest_cached(version, now_utc)
    end

    return archive
end


text:
"""
`renamefile(path::AbstractString, newname)` rename the file-name part in the path only; return the new path of new name.
"""
function renamefile(path::AbstractString, newname)
    dir = dirname(path);
    #  oldname = basename(path);
    return newpath = joinpath(dir, newname)
end

max_stars_repo_name: CGRG-lab/HypertextTools

max_issues_repo_name: CGRG-lab/HypertextTools

max_forks_repo_name: CGRG-lab/HypertextTools


text:
using Base: @propagate_inbounds

mutable struct OffsetVec{T,VV<:AbstractVector{T}} <: AbstractVector{T}
    parent::VV
    offset::Int
end
offset(V::OffsetVec) = V.offset
Base.parent(V::OffsetVec) = V.parent
@propagate_inbounds Base.getindex(V::OffsetVec, i::Int) = parent(V)[i-offset(V)]
@propagate_inbounds function Base.setindex!(V::OffsetVec, val, i::Int)
    parent(V)[i-offset(V)] = val
    V
end

mutable struct OffsetMat{T,MM<:AbstractMatrix{T}} <: AbstractMatrix{T}
    parent::MM
    offset1::Int
    offset2::Int
end
offset1(M::OffsetMat) = M.offset1
offset2(M::OffsetMat) = M.offset2
Base.parent(M::OffsetMat) = M.parent
@propagate_inbounds Base.getindex(M::OffsetMat, i::Int, j::Int) = parent(M)[i-offset1(M), j-offset2(M)]
#@propagate_inbounds Base.getindex(M::OffsetMat, ::Colon, j::Int) = parent(M)[:, j-offset2(M)]


text:
module Frames
import Base

"""
"""
abstract type AbstractCANFrame end

"""
"""
mutable struct CANFrame <: AbstractCANFrame
    frame_id::UInt32
    data::Array{UInt8,1}
    is_extended::Bool

    function CANFrame(frame_id::UInt32, data::AbstractArray{UInt8};
                      is_extended::Bool=false)
        if length(data) > 8
            throw(DomainError(data, "CANFrame allows a maximum of 8 bytes"))
        end
        return new(frame_id, data, is_extended)
    end
end

"""
"""
function CANFrame(frame_id::Integer, data::A; is_extended::Bool=false) where
    {A <: AbstractArray{<:Integer}}
    return CANFrame(convert(UInt32, frame_id), UInt8[data...]; is_extended=is_extended)
end

"""
"""
function CANFrame(frame_id::Integer, data::Integer...; is_extended::Bool=false)
    return CANFrame(convert(UInt32, frame_id), UInt8[data...]; is_extended=is_extended)
end

"""
"""
function CANFrame(frame_id::Integer; is_extended=false)
    return CANFrame(convert(UInt32, frame_id), UInt8[]; is_extended=is_extended)
end

"""
"""
mutable struct CANFdFrame <: AbstractCANFrame
    frame_id::UInt32
    data::Array{UInt8,1}
    is_extended::Bool

    """
    """
    function CANFdFrame(frame_id::UInt32, data::A; is_extended::Bool=false) where
        {A <: AbstractArray{UInt8}}
        if length(data) > 64
            throw(DomainError(data, "CANFdFrame allows a maximum of 64 bytes"))
        end
        return new(frame_id, data, is_extended)
    end
end

"""
"""
function CANFdFrame(frame_id::Integer, data::A; is_extended::Bool=false) where
    {A <: AbstractArray{<:Integer}}
    return CANFdFrame(convert(UInt32, frame_id), UInt8[data...]; is_extended)
end

"""
"""
function CANFdFrame(frame_id::Integer, data::Integer...; is_extended::Bool=false)
    return CANFdFrame(convert(UInt32, frame_id), UInt8[data...]; is_extended)
end

"""
"""
function Base.:(==)(lhs::AbstractCANFrame, rhs::AbstractCANFrame)::Bool
    return false
end

"""
"""
function Base.:(==)(lhs::CANFrame, rhs::CANFrame)::Bool
    if frame_id(lhs) != frame_id(rhs)
        return false
    end

    if data(lhs) != data(rhs)
        return false
    end

    if is_extended(lhs) != is_extended(rhs)
        return false
    end

    return true
end

"""
"""
function frame_id(frame::AbstractCANFrame)::UInt32
    if is_extended(frame)
        return frame.frame_id & 0x1F_FF_FF_FF
    else
        return frame.frame_id & 0x7_FF
    end
end

"""
"""
function data(frame::AbstractCANFrame)::Array{UInt8,1}
    return frame.data
end

"""
"""
function dlc(frame::AbstractCANFrame)::UInt8
    return length(frame.data)
end

"""
"""
function is_standard(frame::AbstractCANFrame)::Bool
    return !frame.is_extended
end

"""
"""
function is_extended(frame::AbstractCANFrame)::Bool
    return frame.is_extended
end

"""
"""
function max_size(::Type{AbstractCANFrame})::UInt8
    return 8
end

"""
"""
function max_size(::Type{CANFrame})::UInt8
    return 8
end

"""
"""
function max_size(::Type{CANFdFrame})::UInt8
    return 64
end

export CANFdFrame, CANFrame
export frame_id, data, dlc, is_extended, is_standard, max_size
end


text:
# Julia wrapper for header: /opt/local/include/libbladeRF.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function bladerf_get_device_list(devices::Ptr{Ptr{bladerf_devinfo}})
    ccall((:bladerf_get_device_list,libbladerf),Cint,(Ptr{Ptr{bladerf_devinfo}},),devices)
end

function bladerf_free_device_list(devices::Ptr{bladerf_devinfo})
    ccall((:bladerf_free_device_list,libbladerf),Void,(Ptr{bladerf_devinfo},),devices)
end

function bladerf_open_with_devinfo(device::Ptr{Ptr{bladerf}},devinfo::Ptr{bladerf_devinfo})
    ccall((:bladerf_open_with_devinfo,libbladerf),Cint,(Ptr{Ptr{bladerf}},Ptr{bladerf_devinfo}),device,devinfo)
end

function bladerf_open(device::Ptr{Ptr{bladerf}},device_identifier::Ptr{Uint8})
    ccall((:bladerf_open,libbladerf),Cint,(Ptr{Ptr{bladerf}},Ptr{Uint8}),device,device_identifier)
end

function bladerf_close(device::Ptr{bladerf})
    ccall((:bladerf_close,libbladerf),Void,(Ptr{bladerf},),device)
end

function bladerf_init_devinfo(info::Ptr{bladerf_devinfo})
    ccall((:bladerf_init_devinfo,libbladerf),Void,(Ptr{bladerf_devinfo},),info)
end

function bladerf_get_devinfo(dev::Ptr{bladerf},info::Ptr{bladerf_devinfo})
    ccall((:bladerf_get_devinfo,libbladerf),Cint,(Ptr{bladerf},Ptr{bladerf_devinfo}),dev,info)
end

function bladerf_get_devinfo_from_str(devstr::Ptr{Uint8},info::Ptr{bladerf_devinfo})
    ccall((:bladerf_get_devinfo_from_str,libbladerf),Cint,(Ptr{Uint8},Ptr{bladerf_devinfo}),devstr,info)
end

function bladerf_devinfo_matches(a::Ptr{bladerf_devinfo},b::Ptr{bladerf_devinfo})
    ccall((:bladerf_devinfo_matches,libbladerf),Bool,(Ptr{bladerf_devinfo},Ptr{bladerf_devinfo}),a,b)
end

function bladerf_devstr_matches(dev_str::Ptr{Uint8},info::Ptr{bladerf_devinfo})
    ccall((:bladerf_devstr_matches,libbladerf),Bool,(Ptr{Uint8},Ptr{bladerf_devinfo}),dev_str,info)
end

function bladerf_backend_str(backend::bladerf_backend)
    ccall((:bladerf_backend_str,libbladerf),Ptr{Uint8},(bladerf_backend,),backend)
end

function bladerf_enable_module(dev::Ptr{bladerf},m::bladerf_module,enable::Bool)
    ccall((:bladerf_enable_module,libbladerf),Cint,(Ptr{bladerf},bladerf_module,Bool),dev,m,enable)
end

function bladerf_set_loopback(dev::Ptr{bladerf},l::bladerf_loopback)
    ccall((:bladerf_set_loopback,libbladerf),Cint,(Ptr{bladerf},bladerf_loopback),dev,l)
end

function bladerf_get_loopback(dev::Ptr{bladerf},l::Ptr{bladerf_loopback})
    ccall((:bladerf_get_loopback,libbladerf),Cint,(Ptr{bladerf},Ptr{bladerf_loopback}),dev,l)
end

function bladerf_set_sample_rate(dev::Ptr{bladerf},_module::bladerf_module,rate::Uint32,actual::Ptr{Uint32})
    ccall((:bladerf_set_sample_rate,libbladerf),Cint,(Ptr{bladerf},bladerf_module,Uint32,Ptr{Uint32}),dev,_module,rate,actual)
end

function bladerf_set_rational_sample_rate(dev::Ptr{bladerf},_module::bladerf_module,rate::Ptr{bladerf_rational_rate},actual::Ptr{bladerf_rational_rate})
    ccall((:bladerf_set_rational_sample_rate,libbladerf),Cint,(Ptr{bladerf},bladerf_module,Ptr{bladerf_rational_rate},Ptr{bladerf_rational_rate}),dev,_module,rate,actual)
end

function bladerf_set_sampling(dev::Ptr{bladerf},sampling::bladerf_sampling)
    ccall((:bladerf_set_sampling,libbladerf),Cint,(Ptr{bladerf},bladerf_sampling),dev,sampling)
end

function bladerf_get_sampling(dev::Ptr{bladerf},sampling::Ptr{bladerf_sampling})
    ccall((:bladerf_get_sampling,libbladerf),Cint,(Ptr{bladerf},Ptr{bladerf_sampling}),dev,sampling)
end

function bladerf_get_sample_rate(dev::Ptr{bladerf},_module::bladerf_module,rate::Ptr{Uint32})
    ccall((:bladerf_get_sample_rate,libbladerf),Cint,(Ptr{bladerf},bladerf_module,Ptr{Uint32}),dev,_module,rate)
end

function bladerf_get_rational_sample_rate(dev::Ptr{bladerf},_module::bladerf_module,rate::Ptr{bladerf_rational_rate})
    ccall((:bladerf_get_rational_sample_rate,libbladerf),Cint,(Ptr{bladerf},bladerf_module,Ptr{bladerf_rational_rate}),dev,_module,rate)
end

function bladerf_set_correction(dev::Ptr{bladerf},_module::bladerf_module,corr::bladerf_correction,value::Int16)
    ccall((:bladerf_set_correction,libbladerf),Cint,(Ptr{bladerf},bladerf_module,bladerf_correction,Int16),dev,_module,corr,value)
end

function bladerf_get_correction(dev::Ptr{bladerf},_module::bladerf_module,corr::bladerf_correction,value::Ptr{Int16})
    ccall((:bladerf_get_correction,libbladerf),Cint,(Ptr{bladerf},bladerf_module,bladerf_correction,Ptr{Int16}),dev,_module,corr,value)
end

function bladerf_set_txvga2(dev::Ptr{bladerf},gain::Cint)
    ccall((:bladerf_set_txvga2,libbladerf),Cint,(Ptr{bladerf},Cint),dev,gain)
end

function bladerf_get_txvga2(dev::Ptr{bladerf},gain::Ptr{Cint})
    ccall((:bladerf_get_txvga2,libbladerf),Cint,(Ptr{bladerf},Ptr{Cint}),dev,gain)
end

function bladerf_set_txvga1(dev::Ptr{bladerf},gain::Cint)
    ccall((:bladerf_set_txvga1,libbladerf),Cint,(Ptr{bladerf},Cint),dev,gain)
end

function bladerf_get_txvga1(dev::Ptr{bladerf},gain::Ptr{Cint})
    ccall((:bladerf_get_txvga1,libbladerf),Cint,(Ptr{bladerf},Ptr{Cint}),dev,gain)
end

function bladerf_set_tx_gain(dev::Ptr{bladerf},gain::Cint)
    ccall((:bladerf_set_tx_gain,libbladerf),Cint,(Ptr{bladerf},Cint),dev,gain)
end

function bladerf_set_lna_gain(dev::Ptr{bladerf},gain::bladerf_lna_gain)
    ccall((:bladerf_set_lna_gain,libbladerf),Cint,(Ptr{bladerf},bladerf_lna_gain),dev,gain)
end

function bladerf_get_lna_gain(dev::Ptr{bladerf},gain::Ptr{bladerf_lna_gain})
    ccall((:bladerf_get_lna_gain,libbladerf),Cint,(Ptr{bladerf},Ptr{bladerf_lna_gain}),dev,gain)
end

function bladerf_set_rxvga1(dev::Ptr{bladerf},gain::Cint)
    ccall((:bladerf_set_rxvga1,libbladerf),Cint,(Ptr{bladerf},Cint),dev,gain)
end

function bladerf_get_rxvga1(dev::Ptr{bladerf},gain::Ptr{Cint})
    ccall((:bladerf_get_rxvga1,libbladerf),Cint,(Ptr{bladerf},Ptr{Cint}),dev,gain)
end

function bladerf_set_rxvga2(dev::Ptr{bladerf},gain::Cint)
    ccall((:bladerf_set_rxvga2,libbladerf),Cint,(Ptr{bladerf},Cint),dev,gain)
end

function bladerf_get_rxvga2(dev::Ptr{bladerf},gain::Ptr{Cint})
    ccall((:bladerf_get_rxvga2,libbladerf),Cint,(Ptr{bladerf},Ptr{Cint}),dev,gain)
end

function bladerf_set_gain(dev::Ptr{bladerf},mod::bladerf_module,gain::Cint)
    ccall((:bladerf_set_gain,libbladerf),Cint,(Ptr{bladerf},bladerf_module,Cint),dev,mod,gain)
end

function bladerf_set_bandwidth(dev::Ptr{bladerf},_module::bladerf_module,bandwidth::Uint32,actual::Ptr{Uint32})
    ccall((:bladerf_set_bandwidth,libbladerf),Cint,(Ptr{bladerf},bladerf_module,Uint32,Ptr{Uint32}),dev,_module,bandwidth,actual)
end

function bladerf_get_bandwidth(dev::Ptr{bladerf},_module::bladerf_module,bandwidth::Ptr{Uint32})
    ccall((:bladerf_get_bandwidth,libbladerf),Cint,(Ptr{bladerf},bladerf_module,Ptr{Uint32}),dev,_module,bandwidth)
end

function bladerf_set_lpf_mode(dev::Ptr{bladerf},_module::bladerf_module,mode::bladerf_lpf_mode)
    ccall((:bladerf_set_lpf_mode,libbladerf),Cint,(Ptr{bladerf},bladerf_module,bladerf_lpf_mode),dev,_module,mode)
end

function bladerf_get_lpf_mode(dev::Ptr{bladerf},_module::bladerf_module,mode::Ptr{bladerf_lpf_mode})
    ccall((:bladerf_get_lpf_mode,libbladerf),Cint,(Ptr{bladerf},bladerf_module,Ptr{bladerf_lpf_mode}),dev,_module,mode)
end

function bladerf_select_band(dev::Ptr{bladerf},_module::bladerf_module,frequency::Uint32)
    ccall((:bladerf_select_band,libbladerf),Cint,(Ptr{bladerf},bladerf_module,Uint32),dev,_module,frequency)
end

function bladerf_set_frequency(dev::Ptr{bladerf},_module::bladerf_module,frequency::Uint32)
    ccall((:bladerf_set_frequency,libbladerf),Cint,(Ptr{bladerf},bladerf_module,Uint32),dev,_module,frequency)
end

function bladerf_get_frequency(dev::Ptr{bladerf},_module::bladerf_module,frequency::Ptr{Uint32})
    ccall((:bladerf_get_frequency,libbladerf),Cint,(Ptr{bladerf},bladerf_module,Ptr{Uint32}),dev,_module,frequency)
end

function bladerf_expansion_attach(dev::Ptr{bladerf},xb::bladerf_xb)
    ccall((:bladerf_expansion_attach,libbladerf),Cint,(Ptr{bladerf},bladerf_xb),dev,xb)
end

function bladerf_expansion_get_attached(dev::Ptr{bladerf},xb::Ptr{bladerf_xb})
    ccall((:bladerf_expansion_get_attached,libbladerf),Cint,(Ptr{bladerf},Ptr{bladerf_xb}),dev,xb)
end

function bladerf_xb200_set_filterbank(dev::Ptr{bladerf},mod::bladerf_module,filter::bladerf_xb200_filter)
    ccall((:bladerf_xb200_set_filterbank,libbladerf),Cint,(Ptr{bladerf},bladerf_module,bladerf_xb200_filter),dev,mod,filter)
end

function bladerf_xb200_get_filterbank(dev::Ptr{bladerf},_module::bladerf_module,filter::Ptr{bladerf_xb200_filter})
    ccall((:bladerf_xb200_get_filterbank,libbladerf),Cint,(Ptr{bladerf},bladerf_module,Ptr{bladerf_xb200_filter}),dev,_module,filter)
end

function bladerf_xb200_set_path(dev::Ptr{bladerf},_module::bladerf_module,path::bladerf_xb200_path)
    ccall((:bladerf_xb200_set_path,libbladerf),Cint,(Ptr{bladerf},bladerf_module,bladerf_xb200_path),dev,_module,path)
end

function bladerf_xb200_get_path(dev::Ptr{bladerf},_module::bladerf_module,path::Ptr{bladerf_xb200_path})
    ccall((:bladerf_xb200_get_path,libbladerf),Cint,(Ptr{bladerf},bladerf_module,Ptr{bladerf_xb200_path}),dev,_module,path)
end

function bladerf_init_stream(stream::Ptr{Ptr{bladerf_stream}},dev::Ptr{bladerf},callback::bladerf_stream_cb,buffers::Ptr{Ptr{Ptr{Void}}},num_buffers::Csize_t,format::bladerf_format,samples_per_buffer::Csize_t,num_transfers::Csize_t,user_data::Ptr{Void})
    ccall((:bladerf_init_stream,libbladerf),Cint,(Ptr{Ptr{bladerf_stream}},Ptr{bladerf},bladerf_stream_cb,Ptr{Ptr{Ptr{Void}}},Csize_t,bladerf_format,Csize_t,Csize_t,Ptr{Void}),stream,dev,callback,buffers,num_buffers,format,samples_per_buffer,num_transfers,user_data)
end

function bladerf_stream(stream::Ptr{bladerf_stream},_module::bladerf_module)
    ccall((:bladerf_stream,libbladerf),Cint,(Ptr{bladerf_stream},bladerf_module),stream,_module)
end

function bladerf_submit_stream_buffer(stream::Ptr{bladerf_stream},buffer::Ptr{Void},timeout_ms::Uint32)
    ccall((:bladerf_submit_stream_buffer,libbladerf),Cint,(Ptr{bladerf_stream},Ptr{Void},Uint32),stream,buffer,timeout_ms)
end

function bladerf_deinit_stream(stream::Ptr{bladerf_stream})
    ccall((:bladerf_deinit_stream,libbladerf),Void,(Ptr{bladerf_stream},),stream)
end

function bladerf_set_stream_timeout(dev::Ptr{bladerf},_module::bladerf_module,timeout::Uint32)
    ccall((:bladerf_set_stream_timeout,libbladerf),Cint,(Ptr{bladerf},bladerf_module,Uint32),dev,_module,timeout)
end

function bladerf_get_stream_timeout(dev::Ptr{bladerf},_module::bladerf_module,timeout::Ptr{Uint32})
    ccall((:bladerf_get_stream_timeout,libbladerf),Cint,(Ptr{bladerf},bladerf_module,Ptr{Uint32}),dev,_module,timeout)
end

function bladerf_sync_config(dev::Ptr{bladerf},_module::bladerf_module,format::bladerf_format,num_buffers::Uint32,buffer_size::Uint32,num_transfers::Uint32,stream_timeout::Uint32)
    ccall((:bladerf_sync_config,libbladerf),Cint,(Ptr{bladerf},bladerf_module,bladerf_format,Uint32,Uint32,Uint32,Uint32),dev,_module,format,num_buffers,buffer_size,num_transfers,stream_timeout)
end

function bladerf_sync_tx(dev::Ptr{bladerf},samples::Ptr{Void},num_samples::Uint32,metadata::Ptr{bladerf_metadata},timeout_ms::Uint32)
    ccall((:bladerf_sync_tx,libbladerf),Cint,(Ptr{bladerf},Ptr{Void},Uint32,Ptr{bladerf_metadata},Uint32),dev,samples,num_samples,metadata,timeout_ms)
end

function bladerf_sync_rx(dev::Ptr{bladerf},samples::Ptr{Void},num_samples::Uint32,metadata::Ptr{bladerf_metadata},timeout_ms::Uint32)
    ccall((:bladerf_sync_rx,libbladerf),Cint,(Ptr{bladerf},Ptr{Void},Uint32,Ptr{bladerf_metadata},Uint32),dev,samples,num_samples,metadata,timeout_ms)
end

function bladerf_get_serial(dev::Ptr{bladerf},serial::Ptr{Uint8})
    ccall((:bladerf_get_serial,libbladerf),Cint,(Ptr{bladerf},Ptr{Uint8}),dev,serial)
end

function bladerf_get_vctcxo_trim(dev::Ptr{bladerf},trim::Ptr{Uint16})
    ccall((:bladerf_get_vctcxo_trim,libbladerf),Cint,(Ptr{bladerf},Ptr{Uint16}),dev,trim)
end

function bladerf_get_fpga_size(dev::Ptr{bladerf},size::Ptr{bladerf_fpga_size})
    ccall((:bladerf_get_fpga_size,libbladerf),Cint,(Ptr{bladerf},Ptr{bladerf_fpga_size}),dev,size)
end

function bladerf_fw_version(dev::Ptr{bladerf},version::Ptr{bladerf_version})
    ccall((:bladerf_fw_version,libbladerf),Cint,(Ptr{bladerf},Ptr{bladerf_version}),dev,version)
end

function bladerf_is_fpga_configured(dev::Ptr{bladerf})
    ccall((:bladerf_is_fpga_configured,libbladerf),Cint,(Ptr{bladerf},),dev)
end

function bladerf_fpga_version(dev::Ptr{bladerf},version::Ptr{bladerf_version})
    ccall((:bladerf_fpga_version,libbladerf),Cint,(Ptr{bladerf},Ptr{bladerf_version}),dev,version)
end

function bladerf_device_speed(dev::Ptr{bladerf})
    ccall((:bladerf_device_speed,libbladerf),bladerf_dev_speed,(Ptr{bladerf},),dev)
end

function bladerf_flash_firmware(dev::Ptr{bladerf},firmware::Ptr{Uint8})
    ccall((:bladerf_flash_firmware,libbladerf),Cint,(Ptr{bladerf},Ptr{Uint8}),dev,firmware)
end

function bladerf_load_fpga(dev::Ptr{bladerf},fpga::Ptr{Uint8})
    ccall((:bladerf_load_fpga,libbladerf),Cint,(Ptr{bladerf},Ptr{Uint8}),dev,fpga)
end

function bladerf_flash_fpga(dev::Ptr{bladerf},fpga_image::Ptr{Uint8})
    ccall((:bladerf_flash_fpga,libbladerf),Cint,(Ptr{bladerf},Ptr{Uint8}),dev,fpga_image)
end

function bladerf_erase_stored_fpga(dev::Ptr{bladerf})
    ccall((:bladerf_erase_stored_fpga,libbladerf),Cint,(Ptr{bladerf},),dev)
end

function bladerf_device_reset(dev::Ptr{bladerf})
    ccall((:bladerf_device_reset,libbladerf),Cint,(Ptr{bladerf},),dev)
end

function bladerf_jump_to_bootloader(dev::Ptr{bladerf})
    ccall((:bladerf_jump_to_bootloader,libbladerf),Cint,(Ptr{bladerf},),dev)
end

function bladerf_strerror(error::Cint)
    ccall((:bladerf_strerror,libbladerf),Ptr{Uint8},(Cint,),error)
end

function bladerf_version(version::Ptr{bladerf_version})
    ccall((:bladerf_version,libbladerf),Void,(Ptr{bladerf_version},),version)
end

function bladerf_log_set_verbosity(level::bladerf_log_level)
    ccall((:bladerf_log_set_verbosity,libbladerf),Void,(bladerf_log_level,),level)
end

function bladerf_alloc_image(_type::bladerf_image_type,address::Uint32,length::Uint32)
    ccall((:bladerf_alloc_image,libbladerf),Ptr{bladerf_image},(bladerf_image_type,Uint32,Uint32),_type,address,length)
end

function bladerf_alloc_cal_image(fpga_size::bladerf_fpga_size,vctcxo_trim::Uint16)
    ccall((:bladerf_alloc_cal_image,libbladerf),Ptr{bladerf_image},(bladerf_fpga_size,Uint16),fpga_size,vctcxo_trim)
end

function bladerf_free_image(image::Ptr{bladerf_image})
    ccall((:bladerf_free_image,libbladerf),Void,(Ptr{bladerf_image},),image)
end

function bladerf_image_write(image::Ptr{bladerf_image},file::Ptr{Uint8})
    ccall((:bladerf_image_write,libbladerf),Cint,(Ptr{bladerf_image},Ptr{Uint8}),image,file)
end

function bladerf_image_read(image::Ptr{bladerf_image},file::Ptr{Uint8})
    ccall((:bladerf_image_read,libbladerf),Cint,(Ptr{bladerf_image},Ptr{Uint8}),image,file)
end

function bladerf_si5338_read(dev::Ptr{bladerf},address::Uint8,val::Ptr{Uint8})
    ccall((:bladerf_si5338_read,libbladerf),Cint,(Ptr{bladerf},Uint8,Ptr{Uint8}),dev,address,val)
end

function bladerf_si5338_write(dev::Ptr{bladerf},address::Uint8,val::Uint8)
    ccall((:bladerf_si5338_write,libbladerf),Cint,(Ptr{bladerf},Uint8,Uint8),dev,address,val)
end

function bladerf_si5338_set_tx_freq(dev::Ptr{bladerf},freq::Uint32)
    ccall((:bladerf_si5338_set_tx_freq,libbladerf),Cint,(Ptr{bladerf},Uint32),dev,freq)
end

function bladerf_si5338_set_rx_freq(dev::Ptr{bladerf},freq::Uint32)
    ccall((:bladerf_si5338_set_rx_freq,libbladerf),Cint,(Ptr{bladerf},Uint32),dev,freq)
end

function bladerf_lms_read(dev::Ptr{bladerf},address::Uint8,val::Ptr{Uint8})
    ccall((:bladerf_lms_read,libbladerf),Cint,(Ptr{bladerf},Uint8,Ptr{Uint8}),dev,address,val)
end

function bladerf_lms_write(dev::Ptr{bladerf},address::Uint8,val::Uint8)
    ccall((:bladerf_lms_write,libbladerf),Cint,(Ptr{bladerf},Uint8,Uint8),dev,address,val)
end

function bladerf_lms_set_dc_cals(dev::Ptr{bladerf},dc_cals::Ptr{bladerf_lms_dc_cals})
    ccall((:bladerf_lms_set_dc_cals,libbladerf),Cint,(Ptr{bladerf},Ptr{bladerf_lms_dc_cals}),dev,dc_cals)
end

function bladerf_lms_get_dc_cals(dev::Ptr{bladerf},dc_cals::Ptr{bladerf_lms_dc_cals})
    ccall((:bladerf_lms_get_dc_cals,libbladerf),Cint,(Ptr{bladerf},Ptr{bladerf_lms_dc_cals}),dev,dc_cals)
end

function bladerf_config_gpio_read(dev::Ptr{bladerf},val::Ptr{Uint32})
    ccall((:bladerf_config_gpio_read,libbladerf),Cint,(Ptr{bladerf},Ptr{Uint32}),dev,val)
end

function bladerf_config_gpio_write(dev::Ptr{bladerf},val::Uint32)
    ccall((:bladerf_config_gpio_write,libbladerf),Cint,(Ptr{bladerf},Uint32),dev,val)
end

function bladerf_expansion_gpio_read(dev::Ptr{bladerf},val::Ptr{Uint32})
    ccall((:bladerf_expansion_gpio_read,libbladerf),Cint,(Ptr{bladerf},Ptr{Uint32}),dev,val)
end

function bladerf_expansion_gpio_write(dev::Ptr{bladerf},val::Uint32)
    ccall((:bladerf_expansion_gpio_write,libbladerf),Cint,(Ptr{bladerf},Uint32),dev,val)
end

function bladerf_expansion_gpio_dir_read(dev::Ptr{bladerf},val::Ptr{Uint32})
    ccall((:bladerf_expansion_gpio_dir_read,libbladerf),Cint,(Ptr{bladerf},Ptr{Uint32}),dev,val)
end

function bladerf_expansion_gpio_dir_write(dev::Ptr{bladerf},val::Uint32)
    ccall((:bladerf_expansion_gpio_dir_write,libbladerf),Cint,(Ptr{bladerf},Uint32),dev,val)
end

function bladerf_get_timestamp(dev::Ptr{bladerf},mod::bladerf_module,value::Ptr{Uint64})
    ccall((:bladerf_get_timestamp,libbladerf),Cint,(Ptr{bladerf},bladerf_module,Ptr{Uint64}),dev,mod,value)
end

function bladerf_dac_write(dev::Ptr{bladerf},val::Uint16)
    ccall((:bladerf_dac_write,libbladerf),Cint,(Ptr{bladerf},Uint16),dev,val)
end

function bladerf_xb_spi_write(dev::Ptr{bladerf},val::Uint32)
    ccall((:bladerf_xb_spi_write,libbladerf),Cint,(Ptr{bladerf},Uint32),dev,val)
end

function bladerf_calibrate_dc(dev::Ptr{bladerf},_module::bladerf_cal_module)
    ccall((:bladerf_calibrate_dc,libbladerf),Cint,(Ptr{bladerf},bladerf_cal_module),dev,_module)
end

function bladerf_erase_flash(dev::Ptr{bladerf},erase_block::Uint32,count::Uint32)
    ccall((:bladerf_erase_flash,libbladerf),Cint,(Ptr{bladerf},Uint32,Uint32),dev,erase_block,count)
end

function bladerf_read_flash(dev::Ptr{bladerf},buf::Ptr{Uint8},page::Uint32,count::Uint32)
    ccall((:bladerf_read_flash,libbladerf),Cint,(Ptr{bladerf},Ptr{Uint8},Uint32,Uint32),dev,buf,page,count)
end

function bladerf_write_flash(dev::Ptr{bladerf},buf::Ptr{Uint8},page::Uint32,count::Uint32)
    ccall((:bladerf_write_flash,libbladerf),Cint,(Ptr{bladerf},Ptr{Uint8},Uint32,Uint32),dev,buf,page,count)
end


text:
"""get all groups"""
function getgrp(node_memory::Vector{Dict{Int,Int}})
    groups = Dict{Int,Vector{Int}}()
    for i=1:length(node_memory)
        for k in keys(node_memory[i])
            if haskey(groups, k)
                push!(groups[k], i)
            else
                groups[k] = [i]
            end
        end
    end
    collect(values(groups))
end

"""
calculate the entropy of a probility distribution
"""
function entropy(P::Vector{Float64})
    H = 0.0
    for i=1:length(P)
        # treat 0*log(0) as being equal to zero
        H += P[i] > 0.0 ? P[i]*log(P[i]) : 0.0
    end
    -H > 0 ? -H : 0.0
end

"""read groups from file"""
function readgrp(filename)
    groups = Vector{Int}[]
    f = open(filename, "r")
    for ll in eachline(f)
        push!(groups, [parse(Int, i) for i in split(chomp(ll))])
    end
    groups
end

"""write groups to file"""
function writegrp(filename, groups::Vector{Vector{Int}})
    f = open(filename, "w")
    for i=1:length(groups)
        for j=1:length(groups[i])-1
            print(f, groups[i][j],' ')
        end
        print(f, groups[i][end],'\n')
    end
    close(f)
end

"""read membership from file"""
function readmsp(filename)
    membership = Dict{Int, Vector{Int}}()
    f = open(filename, "r")
    for ll in eachline(f)
        entries = [parse(Int, i) for i in split(chomp(ll))]
        membership[entries[1]] = entries[2:end]
    end
    membership
end

"""write membership to file"""
function writemsp(filename, membership::Dict{Int, Vector{Int}})
    f = open(filename, "w")
    for k in sort(collect(keys(membership)))
        print(f, k, '\t')
        for j=1:length(membership[k])-1
            print(f, membership[k][j], ' ')
        end
        print(f, membership[k][end], '\n')
    end
    close(f)
end

"""transform membership to groups"""
function msp2grp(membership::Dict{Int, Vector{Int}})
    groups = Dict{Int, Vector{Int}}()
    for (k,v) in membership
        for i in v
            if haskey(groups, i)
                push!(groups[i], k)
            else
                groups[i] = [k]
            end
        end
    end
    collect(values(groups))
end

"""transform groups to membership"""
function grp2msp(groups::Vector{Vector{Int}})
    membership = Dict{Int, Vector{Int}}()
    for i=1:length(groups)
        for j=1:length(groups[i])
            if haskey(membership, groups[i][j])
                push!(membership[groups[i][j]], i)
            else
                membership[groups[i][j]] = [i]
            end
        end
    end
    membership
end